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

偏微分方程数值解上机实习 数值求解二维扩散方程的初边值问题 古典显式格式： 将原格式化为： 附源程序： %-------------------------------------------运用古典显式差分格式求解二维扩散方程的初边值问题; function gdxs(ti,h,t) %-------------------------------------------ti:时间步长; %-------------------------------------------h:空间步长; k=t/ti; m=1/h+1; r=ti/h^2; %------------------------------ r为网格比; w=ones(m,m); u=ones(m,m); for i=2:m-1 for j=2:m-1 u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1)); end end tic for l=1:k for i=2:m-1 for j=2:m-1 w(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); end end u=w; end toc t=toc u mesh(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); 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); end x(n)=y(n)/u(n); %------------------------追赶法求解之赶过程 求解Uz=y; for j=n-1:-1:1 if u(j)==0 break; else x(j)=(y(j)-c(j)*x(j+1))/u(j); end end %-----------------------------------------------运用P-R差分格式求解二维扩散方程的初边值问题; function pr(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); for i=2:m-1 for j=2:m-1 u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1)); end end %------------------------输入用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); fo