一维抛物线方程混合边界问题matlab求解(一维扩散方程)..docx

在差分方法中分为前差、后差、中心差以及显式、隐式、中心式。这些概念分别对应着对空间和时间的离散。针对以上抛物线方程的求解方法，分别采用向前、向后、对称六点和三行式进行计算。一维抛物线的一般方程为此题为混合边界问题1、前差的一般格式为将成为网比，记做r，则差分格式可以写成将上标为k+1的放在一边，k的放在一边，这样就可以写成矩阵形式。根据已知条件便可以求解。（下面是前差的matlab代码）r=10;x=0:0.1:1;t=0:0.1:1;U=[];U(:,1)=0; U(:,11)=0;U(1,:)=sin(pi.*x);for i=2:11 %%行for j=2:10 %% 列 U(i,j)=r*U(i-1,j-1)+(1-2*r)*U(i-1,j)+r*U(i-1,j+1);endendfigure;surf(x,t,U);例题中并没有给出t的具体值，这里取了1，如同正方形laplace方程一样，纯粹为了好求。上面两个for循环代替了矩阵的作用。如果想改成矩阵，很简单，自行解决。2、后差的一般格式为这里请注意前差后差的不同点在于右式中上标的变化。附上matlab代码：h=0.1;k=0.1;N=1/h;M=1/0.1;r=k/h^2;for i=1:N-1B(i)=sin(i/N*pi);endU(:,1)=B;A=diag(ones(1,N-1)*(1+2*r))-diag(ones(1,N-2)*r,1)-diag(ones(1,N-2)*r,-1);for i=2:M+1U(:,i)=A\U(:,i-1);endZ=[zeros(11,1),U,zeros(11,1)];[X,Y]=meshgrid(linspace(1,0,11),linspace(1,0,11));surf(X,Y,Z);对称六点一般格式整理上式可以得到类似AX=BY形式的矩阵，h=0.1;k=0.1;N=1/h;M=1/0.1;r=k/h^2;for i=1:N-1B(i)=sin(i/N*pi)endU(:,1)=BA=diag(ones(1,N-1)*(1+r))-diag(ones(1,N-2)*r/2,1)-diag(ones(1,N-2)*r/2,-1)C=diag(ones(1,N-1)*(1-r))+diag(ones(1,N-2)*r/2,1)+diag(ones(1,N-2)*r/2,-1)for i=2:M+1U(:,i)=A\(C*U(:,i-1))endZ=[zeros(11,1),U,zeros(11,1)];[X,Y]=meshgrid(linspace(0,1,11),linspace(0,1,11));surf(X,Y,Z);三行式的稳定性很差，做出来的结果也令人吃惊，不知是对是错，大家自己好好检查一下。还有一点就是三行式三行式，肯定用到了三行，而目前已知条件是无法用三行式来求解的，我这里选择了用后差法获得了前两行，用前两行来求解第三行，以此类推。当然也可以选择前差或者对称六点，在这里我只是想到了在学龙格库塔方法或者米莱尔—哈明预测校正方法时说隐式比较稳定采用的，具体也不懂。代码如下：h=0.1;k=0.1;N=1/0.1;M=1/0.1;r=k/h^2;U=[];for i=1:N-1B(i)=sin(i/N*pi)endA=diag(ones(1,N-1)*(1+2*r))-diag(ones(1,N-2)*r,1)-diag(ones(1,N-2)*r,-1);U(:,1)=B;U(:,2)=A\U(:,1);I=diag(ones(1,N-1)*(-4*r))+diag(ones(1,N-2)*2*r,1)+diag(ones(1,N-2)*2*r,-1)for i=3:M+1 U(:,i)=I*U(:,i-1)+eye(N-1)*U(:,i-2)endZ=[zeros(11,1),U,zeros(11,1)][X,Y]=meshgrid(linspace(0,1,11),linspace(0,1,11));surf(X,Y,Z);结果不难看出，确实稳定性好差，不过对与不对，大家自己去检查吧，这里只提供最简单的思路。