资源描述：
《偏微分方程数值解上机》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、偏微分方程数值解上机实习数值求解二维扩散方程的初边值问题古典显式格式：将原格式化为：附源程序：%-------------------------------------------运用古典显式差分格式求解二维扩散方程的初边值问题;functiongdxs(ti,h,t)%-------------------------------------------ti:时间步长;%-------------------------------------------h:空间步长;k=t/ti;m=1/h+1;r=ti/h^2;%--------
2、----------------------r为网格比;w=ones(m,m);u=ones(m,m);fori=2:m-1forj=2:m-1u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1));endendticforl=1:kfori=2:m-1forj=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)8交替方向隐式格式（P-R格式）：将原差分格式化为：代
3、入边界条件，转化为三对角矩阵8附追赶法源程序：%-------------------------------------------追赶法求解三对角方程组;functionx=zg(a,b,c,d)%--------------------------------------------a:方程组系数矩阵A的下对角元素;%--------------------------------------------b:方程组系数矩阵A的主对角元素;%--------------------------------------------c:方
4、程组系数矩阵A的上对角元素;%--------------------------------------------d:追赶法所求方程的右端向量;%--------------------------------------------l:系数矩阵A所分解成的下三角阵L中的下对角元素了l(i);%--------------------------------------------u:系数矩阵A所分解成的下三角阵U中的主对角元素了u(i);n=length(b);u(1)=b(1);y(1)=d(1);fori=1:n-1%-----
5、---------------------追赶法求解之追过程求解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;forj=n-1:-1:1ifu(j)==0break;elsex(j)=(y(j)-c(j)*x(j+1))/u(j);endend%---------------------------------------------
6、--运用P-R差分格式求解二维扩散方程的初边值问题;functionpr(ti,h,t)%-------------------------------------------ti:时间步长h:空间步长;k=t/ti+1;m=1/h+1;r=ti/h^2;%------------------------------r为网格比;w=ones(m,m);u=ones(m,m);%------------------------输入初始值v=ones(m,m);fori=2:m-1forj=2:m-1u(i,j)=sin(pi*(i-1)*h
7、)*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);8A=zeros(m-2,m-2);fori=1:m-2A(i,i)=2-2*r;endfori=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;ticforl=1:kfori=2:m-1d1=A*u(i,2:m
8、-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
