偏微分方程数值解法题解

1、偏微分方程数值解法（带程序）例1 求解初边值问题要求采用树脂格式，完成下列计算：（1） 取分别计算时刻的数值解。（2） 取分别计算时刻的数值解。（3） 取分别计算时刻的数值解。并与解析解进行比较。解：程序function A=zhongxinchafen(x,y,la)U=zeros(length(x),length(y);for i=1:size(x,2) if x(i)>0&x(i)<=0.5 U(i,1)=2*x(i); elseif x(i)>0.5&x(i)<1 U(i,1)=2*(1-x(i); endendfor j=1:length(y)

2、-1 for i=1:length(x)-2 U(i+1,j+1)=U(i+1,j)+la*(U(i+2,j)-2*U(i+1,j)+U(i,j); endendA=U(:,size(U,2)function u=jiexijie1(x,t)for i=1:size(x,2) k=3; a1=(1/(12)*sin(1*pi/2)*sin(1*pi*x(i)*exp(-12*pi2*t); a2=a1+(1/(22)*sin(2*pi/2)*sin(2*pi*x(i)*exp(-22*pi2*t); while abs(a2-a1)>0.00001 a1=a2; a2=a1+(1/(k2

3、)*sin(k*pi/2)*sin(k*pi*x(i)*exp(-k2*pi2*t); k=k+1; end u(i)=8/(pi2)*a2;endclc; %第1题第1问clear;t1=0.01;t2=0.02;t3=0.1;x=0:0.1:1;y1=0:0.001:t1;y2=0:0.001:t2;y3=0:0.001:t3;la=0.1;subplot(131)A1=zhongxinchafen(x,y1,la);u1=jiexijie1(x,t1)line(x,A1,'color','r','linestyle',':'

4、,'linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1);A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x,t2)line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,u2,'color','b','linewidth',1);A3=z

5、hongxinchafen(x,y3,la);u3=jiexijie1(x,t3)line(x,A3,'color','r','linestyle',':','linewidth',1.5);line(x,u3,'color','b','linewidth',1);title('例1(1)');subplot(132);line(x,u1,'color','b','linewidth',1);line(x

6、,u2,'color','b','linewidth',1);line(x,u3,'color','b','linewidth',1);title('解析解');subplot(133);line(x,A1,'color','r','linestyle',':','linewidth',1.5);line(x,A2,'color','r','linestyle

7、9;,':','linewidth',1.5);line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解');clc; %第1题第2问clear;t1=0.01;t2=0.02;t3=0.1;x=0:0.1:1;y1=0:0.005:t1;y2=0:0.005:t2;y3=0:0.005:t3;la=0.5;subplot(131);A1=zhongxinchafen(x,y1,la);u1=j

8、iexijie1(x,t1)line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1);A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x,t2)line(x,A2,'color','r','linestyle',':','li

9、newidth',1.5);line(x,u2,'color','b','linewidth',1); A3=zhongxinchafen(x,y3,la);u3=jiexijie1(x,t3)line(x,A3,'color','r','linestyle',':','linewidth',1.5);line(x,u3,'color','b','linewidth',1);title('例1(2)

10、9;);subplot(132);line(x,u1,'color','b','linewidth',1);line(x,u2,'color','b','linewidth',1);line(x,u3,'color','b','linewidth',1);title('解析解');subplot(133);line(x,A1,'color','r','linestyle',':&#

11、39;,'linewidth',1.5);line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解');clc; %第1题第3问clear;t1=0.01;t2=0.02;t3=0.1;x=0:0.1:1;y1=0:0.01

12、:t1;y2=0:0.01:t2;y3=0:0.01:t3;la=1.0;subplot(131);A1=zhongxinchafen(x,y1,la);u1=jiexijie1(x,t1)line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1);A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x

13、,t2)line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,u2,'color','b','linewidth',1);A3=zhongxinchafen(x,y3,la);u3=jiexijie1(x,t3)line(x,A3,'color','r','linestyle',':','linewidth',1.5);

14、line(x,u3,'color','b','linewidth',1);title('例1(3)');subplot(132);line(x,u1,'color','b','linewidth',1);line(x,u2,'color','b','linewidth',1);line(x,u3,'color','b','linewidth',1);title('解析解')

15、;subplot(133);line(x,A1,'color','r','linestyle',':','linewidth',1.5);line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,A3,'color','r','linestyle',':','linewidth',1.5);t

16、itle('数值解');运行结果：表1：取时刻的解析解与数值解时间(t)解析解数值解时间(t)解析解数值解时间(t)解析解数值解0000000.22690.19960.20560.19390.09340.09440.43170.39680.39110.37810.17760.17960.59410.58220.53830.53730.24440.24720.69840.72810.63280.64870.28730.29070.73440.78670.66540.68910.30210.30560.69840.72810.63280.64870.28730.29070.5941

17、0.58220.53830.53730.24440.24720.43170.39680.39110.37810.17760.17960.22690.19960.20560.19390.09340.09440.000000.000000.00000表2：取时刻的解析解与数值解时间(t)解析解数值解时间(t)解析解数值解时间(t)解析解数值解0000000.22690.20000.20560.20000.09340.09490.43170.40000.39110.37500.17760.17170.59410.60000.53830.55000.24440.24840.69840.70000.63

18、280.62500.28730.27780.73440.80000.66540.70000.30210.30710.69840.70000.63280.62500.28730.27780.59410.60000.53830.55000.24440.24840.43170.40000.39110.37500.17760.17170.22690.20000.20560.20000.09340.09490.000000.000000.00000表3：取时刻的解析解与数值解时间(t)解析解数值解时间(t)解析解数值解时间(t)解析解数值解0000000.22690.20000.20560.20000.

19、0934220.2000.43170.40000.39110.40000.1776-453.6000.59410.60000.53830.60000.2444684.6000.69840.80000.63280.40000.2873-861.2000.73440.60000.66541.00000.3021929.0000.69840.80000.63280.40000.2873-861.2000.59410.60000.53830.60000.2444684.6000.43170.40000.39110.40000.1776-453.6000.22690.20000.20560.20000.

20、0934220.2000.000000.000000.00000图1：取时刻的解析解与数值解图2：取时刻的解析解与数值解图3：取时刻的解析解与数值解例2 用Crank-Nicolson格式完成例1的所有任务。解：格式将第层和第层放到等式两端，得，（在程序中用表示），为时间步长，为空间步长（在程序中用来表示）化简：左右两边同时乘以得：程序：function A=granknicolson(x,y,la,a)U=zeros(size(x,2),size(y,2);for i=1:size(x,2) if x(i)>0&x(i)<=0.5 U(i,1)=2*x(i); elsei

21、f x(i)>0.5&x(i)<1 U(i,1)=2*(1-x(i); endendx=0:0.1:1;V1=zeros(size(x,2)-2);V2=zeros(size(x,2)-2);for i=1:size(x,2)-2 for j=1:size(x,2)-2 if i=j V1(i,j)=2+2*a*la; V2(i,j)=2-2*a*la; elseif i=j+1|j=i+1 V1(i,j)=-a*la; V2(i,j)=a*la; end endendfor j=2:size(y,2) U(2:size(x,2)-1,j)=inv(V1)*V2*U(2:s

22、ize(x,2)-1,j-1);endA=U(:,size(U,2);function u=jiexijie1(x,t)for i=1:size(x,2) k=3; a1=(1/(12)*sin(1*pi/2)*sin(1*pi*x(i)*exp(-12*pi2*t); a2=a1+(1/(22)*sin(2*pi/2)*sin(2*pi*x(i)*exp(-22*pi2*t); while abs(a2-a1)>0.00001 a1=a2; a2=a1+(1/(k2)*sin(k*pi/2)*sin(k*pi*x(i)*exp(-k2*pi2*t); k=k+1; end u(i)=8

23、/(pi2)*a2;endclc; %第2题第1问clear;t1=0.01;t2=0.02;t3=0.1;x=0:0.1:1;y1=0:0.001:t1;y2=0:0.001:t2;y3=0:0.001:t3;la=0.1;a=1;subplot(131);A1=granknicolson(x,y1,la,a);u1=jiexijie1(x,t1)line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color&

24、#39;,'b','linewidth',1);A2=granknicolson(x,y2,la,a);u2=jiexijie1(x,t2)line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,u2,'color','b','linewidth',1);A3=granknicolson(x,y3,la,a);u3=jiexijie1(x,t3)line(x,A3,&#

25、39;color','r','linestyle',':','linewidth',1.5);line(x,u3,'color','b','linewidth',1);title('例2(1)');subplot(132);line(x,u1,'color','b','linewidth',1);line(x,u2,'color','b','linewidth',1

26、);line(x,u3,'color','b','linewidth',1);title('解析解');subplot(133);line(x,A1,'color','r','linestyle',':','linewidth',1.5);line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,A3,

27、9;color','r','linestyle',':','linewidth',1.5);title('数值解');clc; %第2题第2问clear;t1=0.01;t2=0.02;t3=0.1;x=0:0.1:1;y1=0:0.005:t1;y2=0:0.005:t2;y3=0:0.005:t3;la=0.5;a=1;subplot(131);A1=granknicolson(x,y1,la,a);u1=jiexijie1(x,t1)line(x,A1,'color','r&#

28、39;,'linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1);A2=granknicolson(x,y2,la,a);u2=jiexijie1(x,t2)line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,u2,'color',

29、9;b','linewidth',1);A3=granknicolson(x,y3,la,a);u3=jiexijie1(x,t3)line(x,A3,'color','r','linestyle',':','linewidth',1.5);line(x,u3,'color','b','linewidth',1);title('例2(2)');subplot(132);line(x,u1,'color','

30、;b','linewidth',1);line(x,u2,'color','b','linewidth',1);line(x,u3,'color','b','linewidth',1);title('解析解');subplot(133);line(x,A1,'color','r','linestyle',':','linewidth',1.5);line(x,A2,'col

31、or','r','linestyle',':','linewidth',1.5);line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解');clc; %第2题第3问clear;t1=0.01;t2=0.02;t3=0.1;x=0:0.1:1;y1=0:0.01:t1;y2=0:0.01:t2;y3=0:0.01:t3;la=1.0;a=1;subplo

32、t(131);A1=granknicolson(x,y1,la,a);u1=jiexijie1(x,t1)line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1);A2=granknicolson(x,y2,la,a);u2=jiexijie1(x,t2)line(x,A2,'color','r'

33、,'linestyle',':','linewidth',1.5);line(x,u2,'color','b','linewidth',1);A3=granknicolson(x,y3,la,a);u3=jiexijie1(x,t3)line(x,A3,'color','r','linestyle',':','linewidth',1.5);line(x,u3,'color','b',&#

34、39;linewidth',1);title('例2(3)');subplot(132);line(x,u1,'color','b','linewidth',1);line(x,u2,'color','b','linewidth',1);line(x,u3,'color','b','linewidth',1);title('解析解');subplot(133);line(x,A1,'color',

35、'r','linestyle',':','linewidth',1.5);line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解');运行结果：表4：取时刻的解析解与Cran

36、k-Nicolson格式数值解时间(t)解析解数值解时间(t)解析解数值解时间(t)解析解数值解0000000.22690.19930.20560.19330.09340.09490.43170.39590.39110.37730.17760.18050.59410.58100.53830.53710.24440.24840.69840.72880.63280.65000.28730.29210.73440.79040.66540.69140.30210.30710.69840.72880.63280.65000.28730.29210.59410.58100.53830.53710.2444

37、0.24840.43170.39590.39110.37730.17760.18050.22690.19930.20560.19330.09340.09490.000000.000000.00000表5：取时刻的解析解与Crank-Nicolson格式数值解时间(t)解析解数值解时间(t)解析解数值解时间(t)解析解数值解0000000.22690.19920.20560.19340.09340.09490.43170.39590.39110.37760.17760.18040.59410.58200.53830.53750.24440.24840.69840.72930.63280.6497

38、0.28730.29200.73440.78790.66540.69060.30210.30700.69840.72930.63280.64970.28730.29200.59410.58200.53830.53750.24440.24840.43170.39590.39110.37760.17760.18040.22690.19920.20560.19340.09340.09490.000000.000000.00000表6：取时刻的解析解与Crank-Nicolson格式数值解时间(t)解析解数值解时间(t)解析解数值解时间(t)解析解数值解0000000.22690.19890.2056

39、0.19360.09340.09480.43170.39560.39110.37890.17760.18030.59410.58340.53830.53970.24440.24820.69840.73810.63280.64610.28730.29180.73440.76910.66540.69210.30210.30690.69840.73810.63280.64610.28730.29180.59410.58340.53830.53970.24440.24820.43170.39560.39110.37890.17760.18030.22690.19890.20560.19360.0934

40、0.09480.000000.000000.00000图4：取时刻的解析解与Crank-Nicolson格式数值解图5：取时刻的解析解与Crank-Nicolson格式数值解图6：取时刻的解析解与Crank-Nicolson格式数值解例3 数值求解初边值问题要求边界条件离散采用（1） 向前（后）一阶差商代替微商。（2） 中心差商（二阶）代替微商。数值格式仍采用，这里取计算注：本例的解析解为其中是的正根。解：程序：function A=yijiechashang(x,y,la)U=zeros(length(x),length(y);for i=2:length(x)-1 U(i,1)=1;end

41、U(1,1)=U(2,1)/1.1;U(length(x),1)=U(length(x)-1,1)/1.1;for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=U(i+1,j)+la*(U(i+2,j)-2*U(i+1,j)+U(i,j); U(1,j+1)=U(2,j+1)/1.1; U(length(x),j+1)=U(length(x)-1,j+1)/1.1; endendA=U(:,size(U,2);function u=jiexijie2(x,t,d)r=zhenggen(d);for i=1:length(x) k=3; a1=e

42、xp(-4*r(1)2*t)*cos(2*r(1)*(x(i)-0.5)/(3+4*r(1)2)*cos(r(1); a2=a1+exp(-4*r(2)2*t)*cos(2*r(2)*(x(i)-0.5)/(3+4*r(2)2)*cos(r(2); while abs(a2-a1)>0.00001 a1=a2; a2=a1+exp(-4*r(k)2*t)*cos(2*r(k)*(x(i)-0.5)/(3+4*r(k)2)*cos(r(k); k=k+1; end u(i)=4*a2;endclc; %第3题第1问clear;t1=0.005;t2=0.01;t3=0.015;t4=0.0

43、2;t5=0.10;t6=0.20;t7=0.25;t8=0.5;t9=1.0;la=0.25;x=0:0.1:1;y1=0:0.0025:t1;y2=0:0.0025:t2;y3=0:0.0025:t3;y4=0:0.0025:t4;y5=0:0.0025:t5;y6=0:0.0025:t6;y7=0:0.0025:t7;y8=0:0.0025:t8;y9=0:0.0025:t9;d=0:0.000005:100;A1=yijiechashang(x,y1,la);u1=jiexijie2(x,t1,d);A2=yijiechashang(x,y2,la);u2=jiexijie2(x,t2

44、,d);A3=yijiechashang(x,y3,la);u3=jiexijie2(x,t3,d);A4=yijiechashang(x,y4,la);u4=jiexijie2(x,t4,d);A5=yijiechashang(x,y5,la);u5=jiexijie2(x,t5,d);A6=yijiechashang(x,y6,la);u6=jiexijie2(x,t6,d);A7=yijiechashang(x,y7,la);u7=jiexijie2(x,t7,d);A8=yijiechashang(x,y8,la);u8=jiexijie2(x,t8,d);A9=yijiechasha

45、ng(x,y9,la);u9=jiexijie2(x,t9,d);subplot(121);line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1);line(x,A2,'color','r','linestyle',':','linewidth',

46、1.5);line(x,u2,'color','b','linewidth',1); line(x,A3,'color','r','linestyle',':','linewidth',1.5);line(x,u3,'color','b','linewidth',1);line(x,A4,'color','r','linestyle',':','li

47、newidth',1.5);line(x,u4,'color','b','linewidth',1);title('例3(1)');legend('数值解','解析解');subplot(122);line(x,A5,'color','r','linestyle',':','linewidth',1.5);line(x,u5,'color','b','linewidth&#

48、39;,1);line(x,A6,'color','r','linestyle',':','linewidth',1.5);line(x,u6,'color','b','linewidth',1);line(x,A7,'color','r','linestyle',':','linewidth',1.5);line(x,u7,'color','b','

49、;linewidth',1);line(x,A8,'color','r','linestyle',':','linewidth',1.5);line(x,u8,'color','b','linewidth',1);line(x,A9,'color','r','linestyle',':','linewidth',1.5);line(x,u9,'color','

50、;b','linewidth',1);legend('数值解','解析解');clc; %第3题第2问clear;t1=0.005;t2=0.01;t3=0.015;t4=0.02;t5=0.10;t6=0.20;t7=0.25;t8=0.5;t9=1.0;la=0.25;x=0:0.1:1;y1=0:0.0025:t1;y2=0:0.0025:t2;y3=0:0.0025:t3;y4=0:0.0025:t4;y5=0:0.0025:t5;y6=0:0.0025:t6;y7=0:0.0025:t7;y8=0:0.0025:t8;y9=0:

51、0.0025:t9;d=0:0.000005:100;A1=zhongxinchafen2(x,y1,la);u1=jiexijie2(x,t1,d);A2=zhongxinchafen2(x,y2,la);u2=jiexijie2(x,t2,d);A3=zhongxinchafen2(x,y3,la);u3=jiexijie2(x,t3,d);A4=zhongxinchafen2(x,y4,la);u4=jiexijie2(x,t4,d);A5=zhongxinchafen2(x,y5,la);u5=jiexijie2(x,t5,d);A6=zhongxinchafen2(x,y6,la);

52、u6=jiexijie2(x,t6,d);A7=zhongxinchafen2(x,y7,la);u7=jiexijie2(x,t7,d);A8=zhongxinchafen2(x,y8,la);u8=jiexijie2(x,t8,d);A9=zhongxinchafen2(x,y9,la);u9=jiexijie2(x,t9,d);subplot(121);line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'c

53、olor','b','linewidth',1); line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,u2,'color','b','linewidth',1);line(x,A3,'color','r','linestyle',':','linewidth',1.5);lin

54、e(x,u3,'color','b','linewidth',1);line(x,A4,'color','r','linestyle',':','linewidth',1.5);line(x,u4,'color','b','linewidth',1);title('例3(2)');legend('数值解','解析解');subplot(122);line(x,A5,'

55、;color','r','linestyle',':','linewidth',1.5);line(x,u5,'color','b','linewidth',1);line(x,A6,'color','r','linestyle',':','linewidth',1.5);line(x,u6,'color','b','linewidth',1);li

56、ne(x,A7,'color','r','linestyle',':','linewidth',1.5);line(x,u7,'color','b','linewidth',1);line(x,A8,'color','r','linestyle',':','linewidth',1.5);line(x,u8,'color','b','linewidt

57、h',1);line(x,A9,'color','r','linestyle',':','linewidth',1.5);line(x,u9,'color','b','linewidth',1);legend('数值解','解析解'); 运行结果：表7：取用向前（后）一阶差商代替微商计算时刻的解析解与数值解时间(t)解析解数值解时间解析解数值解时间(t)解析解数值解0.0050.92500.87340.010.89650.850

58、70.0150.87550.83310.98410.96070.96270.93580.94440.91640.99840.99430.99050.98010.98020.96620.99991.00000.99840.99610.99450.98951.00001.00000.99980.99960.99880.99761.00001.00001.00001.00000.99960.99931.00001.00000.99980.99960.99880.99760.99991.00000.99840.99610.99450.98950.99840.99430.99050.98010.9802

59、0.96620.98410.96070.96270.93580.94440.91640.92500.87340.89650.85070.87550.83310.020.85850.81840.100.71760.68690.200.60400.57090.92860.90030.78280.75560.65910.62800.96950.95320.83420.81020.70290.67350.98910.98170.87130.84980.73480.70670.99670.99410.89360.87380.75410.72680.99850.99730.90110.88180.7606

60、0.73360.99670.99410.89360.87380.75410.72680.98910.98170.87130.84980.73480.70670.96950.95320.83420.81020.70290.67350.92860.90030.78280.75560.65910.62800.85850.81840.71760.68690.60400.57090.250.55460.52060.50.36190.32831.000.15420.13050.60520.57270.39490.36110.16820.14350.64540.61420.42120.38730.17940.15400.67470.64440.44030.40630.18750.16150.69240.66280.45190.41790.19250.16610.69840.66890.45580.42180.19410.16770.69240.66280.45190.41790.19250.16610.67470.64440.44030.40630.18750