2维抛物线方程数值解法(ADI隐式交替法)方法.docVIP

ADI隐式交替法三种解法及误差分析（一般的教材上只说第一种） 理论部分参看孙志忠：偏微分方程数值解法 注意： 最好不要直接看程序，中间很多公式很烦人的（一定要小心），我写了两天，终于写对了。 中间：例如r*(u(i-1,m1,k)+u(i+1,m1,k))形式写成分形式：r*u(i-1,m1,k)+r*u(i+1,m1,k)后面会出错，我也不是很清楚为什么，可能由于舍入误差，或者大数吃掉小数的影响。 下面有三个程序 具体理论看书，先仔细看书（孙志忠：偏微分方程数值解法）或者网上搜一些理论。 Matlab程序： 1.function [u u0 p e x y t]=ADI1(h1,h2,m1,m2,n) %ADI解二维抛物线型偏微分方程(P-R交替隐式，截断) %此程序用的是追赶法解线性方程组 %h1为空间步长，h2为时间步长 %m1,m2分别为x方向，y方向网格数，n为时间网格数 %p为精确解，u为数值解,e为误差 %定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数 x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况% y=(0:m2)*h1+0; t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2; for k=1:n+1 for i=1:m2+1 for j=1:m1+1 f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k)); end end end for i=1:m2+1 for j=1:m1+1 u(i,j,1)=exp(0.5*(x(j)+y(i))); end end for k=1:n+1 for i=1:m2+1 u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u0(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ; end end for k=1:n+1 for j=1:m1+1 u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))]; u0([1 m2+1],j,k)=[exp(0.5*x(j)-t0(k)) exp(0.5*(1+x(j))-t0(k))]; end end r=h2/(h1*h1);r1=2*(1-r);r2=2*(1+r); for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2 a=-r*ones(1,m1-1); c=a;a(1)=0;c(m1-1)=0; b=r2*ones(1,m1-1); d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+... h2*f(i,2,k); for l=2:m1-2 d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+... h2*f(i,l+1,k); %输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性 end d(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))... +r1*u(i,m1,k)+h2*f(i,m1,k); for l=1:m1-2 %开始解线性方程组 消元过程 a(l+1)=-a(l+1)/b(l); b(l+1)=b(l+1)+a(l+1)*c(l); d(l+1)=d(l+1)+a(l+1)*d(l); end u0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程% for l=m1-2:-1:1 u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l); end end for j=2:m1 a=-r*ones(1,m2-