偏微分方程数值解上机实习.DOC

1、偏微分方程数值解上机实习数值求解二维扩散方程的初边值问题22(,0)sini(1,),utxyyttu (0,10),()01,0xytytx古典显式格式：1,1,1,1,22nnnnjljljljljljljljluuhh将原格式化为： 1 2,1,1,1,1,(4)/nnnnnjljljljljl jlu h中附源程序：%-运用古典显式差分格式求解二维扩散方程的初边值问题;function gdxs(ti,h,t)%-ti:时间步长;%-h:空间步长;k=t/ti;m=1/h+1;r=ti/h2; %- r 为网格比;w=ones(m,m);u=ones(m,m);for i=2:m-1f

2、or j=2:m-1u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1);endendticfor l=1:kfor i=2:m-1for j=2:m-1w(i,j)=r*u(i-1,j)+r*u(i,j-1)+r*u(i+1,j)+r*u(i,j+1)+(1-4*r)*u(i,j);endendu=w;endtoct=tocumesh(u)交替方向隐式格式（P-R 格式）：111222,12111122,2nnnnjljljljljljljljlnnnnjljljljljljljljluuuhh将原差分格式化为： 111222, ,1,121222,1,1, , /

3、nnnnjl jljljl jljljl jljljl jljluuuuh 中+-中 中1, ,12 2,1 ,2,22 222n nl jl njmnmlu uuu AA AAAA 121 ,2 ,1 1, 2,22nljn lnjm nmluu AA代入边界条件，转化为三对角矩阵 12, ,233,1 ,12,2 202 n nl jl njmnmlu uuu AA AAA 21 ,2 13 3,1 1, 2,20nljn lnjm nmluu AA AAA附追赶法源程序：%-追赶法求解三对角方程组;function x=zg(a,b,c,d)%-a:方程组系数矩阵 A 的下对角元素;%-

4、b:方程组系数矩阵 A 的主对角元素;%-c:方程组系数矩阵 A 的上对角元素;%-d:追赶法所求方程的右端向量;%-l:系数矩阵 A 所分解成的下三角阵 L 中的下对角元素了 l(i);%-u:系数矩阵 A 所分解成的下三角阵 U 中的主对角元素了 u(i);n=length(b);u(1)=b(1);y(1)=d(1);for i=1:n-1 %-追赶法求解之追过程 求解 Ly=d;l(i)=a(i)/u(i);u(i+1)=b(i+1)-l(i)*c(i);y(i+1)=d(i+1)-l(i)*y(i);endx(n)=y(n)/u(n); %-追赶法求解之赶过程 求解 Uz=y;for

5、 j=n-1:-1:1if u(j)=0 break;elsex(j)=(y(j)-c(j)*x(j+1)/u(j);endend%-运用 P-R 差分格式求解二维扩散方程的初边值问题;function pr(ti,h,t) %-ti:时间步长 h:空间步长;k=t/ti+1;m=1/h+1;r=ti/h2; %- r 为网格比;w=ones(m,m);u=ones(m,m); %-输入初始值v=ones(m,m);for i=2:m-1for j=2:m-1u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1);endend%-输入用 P-R 差分格式求解的三对角矩阵

6、b=ones(1,m-2)*(2+2*r);a=-r*ones(1,m-3);c=-r*ones(1,m-3);A=zeros(m-2,m-2);for i=1:m-2A(i,i)=2-2*r; endfor i=1:m-3A(i,i+1)=r;A(i+1,i)=r;endp=zeros(m-2,1);p(1)=2*r;p(m-2)=2*r;ticfor l=1:kfor i=2:m-1 d1=A*u(i,2:m-1)+p;d1=d1;w(2:m-1,i)=zg(a,b,c,d1); %-调用追赶法求解d2=A*w(2:m-1,i)+p;v(i,2:m-1)=zg(a,b,c,d2); %-调

7、用追赶法求解end u=v;endtoc t=tocumesh(0:0.1:1,0:0.1:1,u)局部一维格式：111222,1,1,211111222,2( )( )nnnnjljljljljljljljlnnnnjljljljljljljljluuuuhh将原格式化为： 111222, ,1, ,1, 2111222, , ,()()/nnnnjljljljljljljl jljljl jljluuuuh 中12, 1,2,1 ,2,2 2222n nl ll nmlnmlu uuu A AAAAA 121 ,2 ,1 1, 2,22jjn njnjm njmu u AA代入边界条件，转

8、化为三对角矩阵 12, 1,23,1 ,2,2 202n nl ll nmlnml uu A AAA 121 ,3 ,1 1, 2,202 2jjn njnjm njmuu u AA AA附源程序：%-运用局部一维格式求解二维扩散方程的初边值问题;function god(ti,hi,t) %-ti 为时间步长 ， hi 为空间步长;m=1/hi;n=t/ti;g=ti/(hi2); %- g 为网格比u=ones(m+1,m+1); %-输入初始值for i=2:mfor j=2:mu(i,j)=sin(pi*(i-1)*hi)*sin(2*pi*(j-1)*hi);endenda(1:m-

9、2)=-0.5*g;b(1:m-1)=1+g;c(1:m-2)=-0.5*g; %-输入用局部一维差分格式求解的三对角矩阵B=zeros(m-1,m+1);for i=1:m-1B(i,i)=0.5*g;B(i,i+1)=1-g; B(i,i+2)=0.5*g;endf=zeros(m-1,1);f(1,1)=0.5*g;f(m-1,1)=0.5*g;w=ones(m+1,m+1);for i=1:nfor j=2:m d=B*u(:,j)+f;%-调用追赶法求解x=zg(a,b,c,d); w(2:m,j)=x;endfor j=2:me=B*w(j,:)+f;x=zg(a,b,c,e);

10、%-调用追赶法求解u(j,2:m)=x;end endumesh(u)古典显式在 t=1 时运行结果： gdxs（0.0025，0.1，1）所用时间 t=0 1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999707990.99999999944

11、5260.999999999235500.999999999102410.999999999055020.999999999102400.999999999235500.999999999445260.999999999707991.000000000000001.000000000000000.999999999445260.999999998943480.999999998547660.999999998290520.999999998204810.999999998290520.999999998547660.999999998943480.999999999445261.0000000

12、00000001.000000000000000.999999999235500.999999998547660.999999997998510.999999997650070.999999997526020.999999997650070.999999997998510.999999998547660.999999999235501.000000000000001.000000000000000.999999999102400.999999998290520.999999997650070.999999997234010.999999997095320.999999997234010.999

13、999997650070.999999998290520.999999999102401.000000000000001.000000000000000.999999999055020.999999998204810.999999997526020.999999997095320.999999996941990.999999997095320.999999997526020.999999998204810.999999999055021.000000000000001.000000000000000.999999999102400.999999998290520.999999997650070

14、.999999997234010.999999997095320.999999997234010.999999997650070.999999998290520.999999999102401.000000000000001.000000000000000.999999999235500.999999998547660.999999997998510.999999997650070.999999997526020.999999997650070.999999997998510.999999998547660.999999999235501.000000000000001.00000000000

15、0000.999999999445260.999999998943480.999999998547660.999999998290520.999999998204810.999999998290520.999999998547660.999999998943480.999999999445261.000000000000001.000000000000000.999999999707990.999999999445260.999999999235500.999999999102400.999999999055020.999999999102400.999999999235500.9999999

16、99445260.999999999707991.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000000P-R 格式 t=1 时运行结果：pr（0.0025，0.1，1）所用时间 t=0.360000000000001.000000000000001.00000000000

17、0001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999544420.999999999133430.999999998807270.999999998597870.999999998525710.999999998597870.999999998807270.999999999133430.9999999

18、99544421.000000000000001.000000000000000.999999999133430.999999998351690.999999997731300.999999997332980.999999997195730.999999997332980.999999997731300.999999998351690.999999999133431.000000000000001.000000000000000.999999998807270.999999997731300.999999996877400.999999996329170.999999996140260.999

19、999996329170.999999996877400.999999997731300.999999998807271.000000000000001.000000000000000.999999998597870.999999997332980.999999996329170.999999995684680.999999995462600.999999995684680.999999996329170.999999997332980.999999998597871.000000000000001.000000000000000.999999998525710.999999997195730

20、.999999996140260.999999995462600.999999995229100.999999995462600.999999996140260.999999997195730.999999998525711.000000000000001.000000000000000.999999998597870.999999997332980.999999996329170.999999995684680.999999995462600.999999995684680.999999996329170.999999997332980.999999998597871.00000000000

21、0001.000000000000000.999999998807270.999999997731300.999999996877400.999999996329170.999999996140260.999999996329170.999999996877400.999999997731300.999999998807271.000000000000001.000000000000.999999999130.999999998350.999999997730.999999997330.999999997190.999999997330.999999997730.999999998350.99

22、9999999131.00000000000000 343 169 130 298 573 298 130 169 343 0001.000000000000000.999999999544420.999999999133430.999999998807270.999999998597870.999999998525710.999999998597870.999999998807270.999999999133430.999999999544421.000000000000001.000000000000001.000000000000001.000000000000001.000000000

23、000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000000局部一维格式 t=1 时的运行结果： god（0.0025，0.1，1）所用时间 t= 0.390000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000

24、000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999521570.999999999089960.999999998747440.999999998527530.999999998451750.999999998527530.999999998747440.999999999089960.999999999521571.000000000000001.000000000000000.999999999089960.999999998269010.999999997617490.99999

25、9997199200.999999997055060.999999997199200.999999997617490.999999998269010.999999999089961.000000000000001.000000000000000.999999998747440.999999997617490.999999996720760.999999996145030.999999995946640.999999996145030.999999996720760.999999997617500.999999998747441.000000000000001.000000000000000.9

26、99999998527530.999999997199200.999999996145030.999999995468210.999999995234990.999999995468210.999999996145030.999999997199200.999999998527531.000000000000001.000000000000000.999999998451750.999999997055060.999999995946640.999999995234990.999999994989770.999999995234990.999999995946640.9999999970550

27、60.999999998451751.000000000000001.000000000000000.999999998527530.999999997199200.999999996145030.999999995468210.999999995234990.999999995468210.999999996145030.999999997199200.999999998527531.000000000000001.000000000000000.999999998747440.999999997617490.999999996720760.999999996145030.999999995

28、946640.999999996145030.999999996720760.999999997617500.999999998747441.000000000000001.000000000000000.999999999089960.999999998269010.999999997617490.999999997199200.999999997055060.999999997199200.999999997617500.999999998269010.999999999089961.000000000000001.000000000000000.999999999521570.99999

29、9999089960.999999998747440.999999998527530.999999998451750.999999998527530.999999998747440.999999999089960.999999999521571.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000000结论：由上面的表格数据可知：古典显式格式的计算速度最快，且当 时，才是稳定的，局部一维1/4格式的计算速度最慢，但是它是无条件稳定的，P-R 格式的速度居中，也是无条件稳定的。