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

ADI隐式交替法三种解法及误差分析（一般的教材上只说第一种）理论部分参看孙志忠：偏微分方程数值解法注意：1.最好不要直接看程序，中间很多公式很烦人的（一定要小心），我写了两天，终于写对了。2.中间：例如r*(u(i-1,m1,k)+u(i+1,m1,k))形式写成分形式：r*u(i-1,m1,k)+r*u(i+1,m1,k)后面会出错，我也不是很清楚为什么，可能由于舍入误差，或者大数吃掉小数的影响。3.下面有三个程序4.具体理论看书，先仔细看书（孙志忠：偏微分方程数值解法）或者网上搜一些理论。Matlab程序：1.function[uu0pexyt]=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;fork=1:n+1fori=1:m2+1forj=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfori=1:m2+1forj=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfork=1:n+1fori=1:m2+1u(i,[1m1+1],k)=[exp(0.5*y(i)-t(k))exp(0.5*(1+y(i))-t(k))];u0(i,[1m1+1],k)=[exp(0.5*y(i)-t0(k))exp(0.5*(1+y(i))-t0(k))];endendfork=1:n+1forj=1:m1+1u([1m2+1],j,k)=[exp(0.5*x(j)-t(k))exp(0.5*(1+x(j))-t(k))];u0([1m2+1],j,k)=[exp(0.5*x(j)-t0(k))exp(0.5*(1+x(j))-t0(k))];endendr=h2/(h1*h1);r1=2*(1-r);r2=2*(1+r);fork=1:n%外循环，先固定每一时间层，每一时间层上解一线性方程组%fori=2:m2a=-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);forl=2:m1-2d(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的矩阵元素不输入%一定要注意输入元素的正确性endd(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);forl=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);endu0(i,m1,k)=d(m1-1)/b(m1-1);%回代过程%forl=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendforj=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r2*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);forl=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+...h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);forl=1:m2-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);endu(m2,j,k+1)=d(m2-1)/b(m2-1);%回代过程%forl=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfork=1:n+1fori=1:m2+1forj=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k));%p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k));%e为误差endendend2.function[upexyt]=ADI2(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(D'Yakonov交替方向隐格式)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u'(i,j,k)(引入的过渡层),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;%定义t0是为了f(x,y,t)~=0的情况%fork=1:n+1fori=1:m2+1forj=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));%编程时-t0(k)写成了+t0(k),导致错误;endendend%初始条件fori=1:m2+1forj=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endend%边界条件fork=1:n+1fori=1:m2+1u(i,[1m1+1],k)=[exp(0.5*y(i)-t(k))exp(0.5*(1+y(i))-t(k))];endendr=h2/(h1*h1);r4=1+r;r5=r/2;fork=1:nfori=2:m2u0(i,[1m1+1],k)=r4*u(i,[1m1+1],k+1)-r5*(u(i-1,[1m1+1],...k+1)+u(i+1,[1m1+1],k+1));endendfork=1:n+1forj=1:m1+1u([1m2+1],j,k)=[exp(0.5*x(j)-t(k))exp(0.5*(1+x(j))-t(k))];endendr1=r-r*r;r2=2*(r-1)*(r-1);r3=r*r/2;fork=1:n%外循环，先固定每一时间层，每一时间层上解一线性方程组%fori=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=2*r4*ones(1,m1-1);d(1)=r*u0(i,1,k)+r1*(u(i-1,2,k)+u(i,1,k)+u(i+1,2,k)+...u(i,3,k))+r2*u(i,2,k)+r3*(u(i-1,1,k)+...u(i+1,1,k)+u(i-1,3,k)+u(i+1,3,k))+2*h2*f(i,2,k);forl=2:m1-2d(l)=r1*(u(i-1,l+1,k)+u(i,l,k)+u(i+1,l+1,k)+...u(i,l+2,k))+r2*u(i,l+1,k)+r3*(u(i-1,l,k)+...u(i+1,l,k)+u(i-1,l+2,k)+u(i+1,l+2,k))+2*h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r1*(u(i-1,m1,k)+u(i,m1-1,k)+...u(i+1,m1,k)+u(i,m1+1,k))+r2*u(i,m1,k)+...r3*(u(i-1,m1-1,k)+...u(i+1,m1-1,k)+u(i-1,m1+1,k)+u(i+1,m1+1,k))+2*h2*f(i,m1,k);forl=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);forl=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendforj=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=2*r4*ones(1,m2-1);d(1)=r*u(1,j,k+1)+2*u0(2,j,k);forl=2:m2-2d(l)=2*u0(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=2*u0(m2,j,k)+r*u(m2+1,j,k+1);forl=1:m2-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);endu(m2,j,k+1)=d(m2-1)/b(m2-1);%回代过程%forl=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfork=1:n+1fori=1:m2+1forj=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k));%p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k));%e为误差endendend3.function[uu0pexyt]=ADI5(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;fork=1:n+1fori=1:m2+1forj=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfori=1:m2+1forj=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfork=1:n+1fori=1:m2+1u(i,[1m1+1],k)=[exp(0.5*y(i)-t(k))exp(0.5*(1+y(i))-t(k))];u1(i,[1m1+1],k)=[exp(0.5*y(i)-t0(k))exp(0.5*(1+y(i))-t0(k))];endendr=h2/(h1*h1);r1=2*(1-r);r2=r/4;r3=2*(1+r);fork=1:nfori=2:m2u0(i,[1m1+1],k)=u1(i,[1m1+1],k)-r2*(u(i-1,[1m