微分方程數(shù)值解法編程作業(yè)二

1、微分方程數(shù)值解法編程作業(yè)二 微分方程數(shù)值解法編程作業(yè)二考慮擴散方程的初邊值問題 考慮擴散方程的初邊值問題? ?? ? ? ?22 ,0 1, 0;,0 sin ,0 1;0, 1, 0, 0.u u x t t xu x x xu t u t t?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?四種差分格式為： 四種差分格式為：（1）向前差分格式 ）向前差分格式? ?11 1 1 2 n n n

2、 nj j j j u u u u ? ? ? ?? ? ? ? ? ?（2）向后差分格式 ）向后差分格式? ?1 1 11 1 1 2 n n n nj j j j u u u u ? ? ? ? ? ?? ? ? ? ? ? ?（3）Crank-Nicolson Crank-Nicolson 格式 格式? ? ? ?1 1 11 1 1 1 1 1 2 2 2 2n n n n n nj j j j j j u u u u u u

3、? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?（4）Du Du Fort-Frankel Fort-Frankel 格式 格式? ? ? ?1 11 1 1 2 1 2 2 2 n n n nj j j j u u u u ? ? ? ? ? ?? ? ? ? ? ? ?其中， 其中， （1）為二步顯式格式， ）為二步顯式格式， （2） 、 （3）為二步隱式格式， ）為二步隱式格式， （4）為三步顯式格

4、式。 ）為三步顯式格式。 （2） 、（3）格式需要求解線性方程組。用四種方法計算數(shù)值解，并與初值問題的解析解相比較， ）格式需要求解線性方程組。用四種方法計算數(shù)值解，并與初值問題的解析解相比較，得出你的結(jié)論并簡要說明。 得出你的結(jié)論并簡要說明。解： 解：（1）function kuosan1 %向前差分h=0.1;lambna=[0.1 0.5];x=1;t=0.1;J=round(x/h);x=0:h:

5、1;u=sin(pi*x');for k=1:2N=round(t/(lambna(k)*h^2));for i=1:Nfor j=2:Ju(j,i+1)=u(j,i)+lambna(k)*(u(j+1,i)-2*u(j,i)+u(j-1,i));endX=0:0.05:1;U=exp(-pi^2*0.1)*sin(pi*X);x=0:0.1:1;[u1]=CN(0.1,100);[u2]=CN(0.5,20); plot(x,

6、u1,'r:',x,u2,'g-.',X,U,'b','LineWidth',2);xlabel('x');ylabel('u');legend('\lambda=0.1','\lambda=0.5','準確解');title('C-N格式');function [u]=CN(A

7、,T) x=0:0.1:1;m=length(x);u=sin(pi*x);C1=sparse(1:m-2,1:m-2,1+A,m-2,m-2); C2=sparse(2:m-2,1:m-3,-0.5*A,m-2,m-2);C3=sparse(1:m-3,2:m-2,-0.5*A,m-2,m-2);C=C1+C2+C3; %系數(shù)矩陣for n=1:Tb=(1-A)*u(2:m-1)+(0.5*

8、A)*u(3:m)+(0.5*A)*u(1:m-2);u(2:m-1)=C\b';end（4）function kuosan4 %Du Fort-Frankel 格式 格式X=0:0.05:1;U=exp(-pi^2*0.1)*sin(pi*X);x=0:0.1:1;[u1]=DF(0.1,100);[u2]=DF(0.5,20);plot(x,u1,'r:',x,u2,'

9、;g-.',X,U,'b','LineWidth',2);xlabel('x');ylabel('u');legend('\lambda=0.1','\lambda=0.5','準確解 準確解');title('Du Fort-Frankel 格式 格式');function [u]=DF(A,T)