《数值求解二维扩散方程的初边值问题.doc》由会员分享，可在线阅读，更多相关《数值求解二维扩散方程的初边值问题.doc（6页珍藏版）》请在金锄头文库上搜索。
1、 数值求解二维扩散方程的初边值问题古典显式格式:将原格式化为:附源程序:%-运用古典显式差分格式求解二维扩散方程的初边值问题;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-1for 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
2、)+r*u(i,j+1)+(1-4*r)*u(i,j);endendu=w;endtoct=tocumesh(u)交替方向隐式格式(P-R格式):将原差分格式化为:代入边界条件,转化为三对角矩阵附追赶法源程序:%-追赶法求解三对角方程组;function x=zg(a,b,c,d)%-a:方程组系数矩阵A的下对角元素;%-b:方程组系数矩阵A的主对角元素;%-c:方程组系数矩阵A的上对角元素;%-d:追赶法所求方程的右端向量;%-l:系数矩阵A所分解成的下三角阵L中的下对角元素了l(i);%-u:系数矩阵A所分解成的下三角阵U中的主对角元素了u(i);n=length(b);u(1)=b(1)
3、;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 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; %
4、- 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差分格式求解的三对角矩阵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
5、;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); %-调用追赶法求解end u=v;endtoc t=tocumesh(0:0.1:1,0:0.1:1,u)局部一维格式:将原格式化为:代入边界条件,转化为三对角矩阵附源程序:%-运用局部一维格式求解二维扩散方程的初边值问题;function god(ti,hi,t) %-ti为时间步长 , hi为空间步长;m=1/hi;n=
6、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-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
7、=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); %-调用追赶法求解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.00000
8、0000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999707990.999999999445260.999999999235500.999999999102410.999999999055020.999999999102400.999999999235500.999999999445260.999999999707991.000000000000001.000000000000000.999999999445260.999999998943480.9
9、99999998547660.999999998290520.999999998204810.999999998290520.999999998547660.999999998943480.999999999445261.000000000000001.000000000000000.999999999235500.999999998547660.999999997998510.999999997650070.999999997526020.999999997650070.999999997998510.999999998547660.999999999235501.0000000000000
10、01.000000000000000.999999999102400.999999998290520.999999997650070.999999997234010.999999997095320.999999997234010.999999997650070.999999998290520.999999999102401.000000000000001.000000000000000.999999999055020.999999998204810.999999997526020.999999997095320.999999996941990.999999997095320.999999997526020.999999998204810.999999999055021.000000000000001.000000000000000.999999999102400.999999998290520.99999999765007 温馨提示：如果这篇偏微分方程数值解上机文档不完整，请从下面的
