二维抛物线方程数值解法(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 [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-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);
for l=2:m2-2
d(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的矩阵元素不输入%一定要注意输入元素的正确性
end
d(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);
for l=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);
end
u(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%
for l=m2-2:-1:1
u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);
end
end
end
for k=1:n+1
for i=1:m2+1
for j=1:m1+1
p(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为误差end
end
end
2.function [u p e x y t]=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的情况% 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));
%编程时-t0(k)写成了+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))];
end
end
r=h2/(h1*h1);r4=1+r;r5=r/2;
for k=1:n
for i=2:m2
u0(i,[1 m1+1],k)=r4*u(i,[1 m1+1],k+1)-r5*(u(i-1,[1 m1+1],... k+1)+u(i+1,[1 m1+1],k+1));
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))];
end
end
r1=r-r*r;r2=2*(r-1)*(r-1);r3=r*r/2;
for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2
a=-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);
for l=2:m1-2
d(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的矩阵元素不输入%一定要注意输入元素的正确性
end
d(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);
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-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);
for l=2:m2-2
d(l)=2*u0(l+1,j,k);
%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性
end
d(m2-1)=2*u0(m2,j,k)+r*u(m2+1,j,k+1);
for l=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);
end
u(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%
for l=m2-2:-1:1
u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);
end
end
end
for k=1:n+1
for i=1:m2+1
for j=1:m1+1
p(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为误差end
end
end
3.function [u u0 p e x y t]=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;
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))]; u1(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;
end
end
r=h2/(h1*h1);r1=2*(1-r);r2=r/4;r3=2*(1+r);
for k=1:n
for i=2:m2
u0(i,[1 m1+1],k)=u1(i,[1 m1+1],k)-r2*(u(i-1,[1 m1+1],k+1)-...
2*u(i,[1 m1+1],k+1)+u(i+1,[1 m1+1],k+1)-u(i-1,[1 m1+1],k)+...
2*u(i,[1 m1+1],k)-u(i+1,[1 m1+1],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))];
end
end
for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组%
for i=2:m2
a=-r*ones(1,m1-1);
c=a;a(1)=0;c(m1-1)=0;
b=r3*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-1);
c=a;a(1)=0;c(m2-1)=0;
b=r3*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);
for l=2:m2-2
d(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的矩阵元素不输入%一定要注意输入元素的正确性
end
d(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);
for l=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);
end
u(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%
for l=m2-2:-1:1
u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);
end
end
end
for k=1:n+1
for i=1:m2+1
for j=1:m1+1
p(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为误差end
end
end
[up e x y t]=ADI2(0.01,0.001,100,100,1000);surf(x,y,e(:,:,1001)) t=1的误差曲面
下面是三种方法的误差比较：
1.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2）代替）（t=1时的误差）
2.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)
(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)
3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) surf(x,y,e(:,:,11))(表示t=1时的误差）
下面是相关数据：
1: [u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)
截断中间过渡层用u(i,j,k+1/2)代替）
e(:,:,11) =
Columns 1 through 6
0 0 0 0 0 0
0 0.00040947 0.00025182 0.00019077 0.00017112 0.00017604
0 0.00057359 0.00042971 0.00035402 0.00032565 0.00033628
0 0.00066236 0.00054689 0.00047408 0.00044596 0.00046267
0 0.00072152 0.00062001 0.00055081 0.00052442 0.00054553
0 0.00076164 0.0006576 0.00058522 0.00055732 0.00057984
0 0.00078336 0.00065993 0.00057557 0.00054161 0.00056209
0 0.00078161 0.00061872 0.00051646 0.00047429 0.00048964
0 0.00073621 0.0005148 0.00039979 0.00035439 0.00036313
0 0.00056964 0.00031688 0.00022051 0.0001884 0.00019192
0 0 0 0 0 0
2.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)
(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)
e(:,:,11) =
Columns 1 through 6
0 0 0 0 0 0
0 0.00027006 0.00016305 0.00012104 0.0001071 0.00010995
0 0.00037754 0.00027817 0.0002253 0.00020483 0.00021116
0 0.00043539 0.00035386 0.00030207 0.00028124 0.0002914
0 0.00047398 0.00040104 0.00035113 0.00033111 0.00034405
0 0.0005003 0.00042535 0.00037309 0.0003519 0.00036571
0 0.00051479 0.00042699 0.00036681 0.00034164 0.0003541
0 0.00051415 0.00040056 0.00032887 0.0002985 0.00030764
0 0.00048504 0.0003335 0.00025411 0.0002221 0.00022706
0 0.00037609 0.00020532 0.00013956 0.00011718 0.00011902
0 0 0 0 0 0
3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式)
e(:,:,11) =
Columns 1 through 6
0 0 0 0 0 0
0 8.6469e-006 1.4412e-005 1.8364e-005 2.091e-005 2.2174e-005
0 1.4412e-005 2.4777e-005 3.2047e-005 3.6716e-005 3.8961e-005
0 1.8364e-005 3.2047e-005 4.1789e-005 4.8054e-005 5.1008e-005
0 2.091e-005 3.6716e-005 4.8054e-005 5.5353e-005 5.8764e-005
0 2.2174e-005 3.8961e-005 5.1008e-005 5.8764e-005 6.2389e-005
0 2.2118e-005 3.8698e-005 5.0523e-005 5.8126e-005 6.171e-005
0 2.055e-005 3.5581e-005 4.6157e-005 5.2942e-005 5.6197e-005
0 1.707e-005 2.8951e-005 3.7128e-005 4.2365e-005 4.4952e-005
0 1.0851e-005 1.7698e-005 2.2265e-005 2.5203e-005 2.672e-005
0 0 0 0 0 0
1.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)
截断中间过渡层用u(i,j,k+1/2)代替）
Columns 7 through 11
0 0 0 0 0
0.00020348 0.00026228 0.00038338 0.00066008 0
0.00038607 0.00048321 0.00064717 0.00091668 0
0.00052635 0.00064203 0.00081637 0.0010517 0
0.0006174 0.00074272 0.00092111 0.0011417 0
0.00065651 0.00078964 0.00097724 0.0012051 0
0.00064051 0.00078116 0.00098594 0.0012433 0
0.00056474 0.00070822 0.00093332 0.0012478 0
0.00042547 0.00055526 0.00078616 0.0011844 0
0.00022735 0.00030946 0.00049004 0.00092402 0
0 0 0 0 0
2.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)
(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)
Columns 7 through 11
0 0 0 0 0
0.00012826 0.00016798 0.00024986 0.00043637 0
0.00024444 0.00031023 0.00042173 0.00060513 0
0.00033401 0.00041257 0.00053179 0.00069358 0
0.00039216 0.00047742 0.00059986 0.00075263 0
0.00041704 0.00050761 0.00063642 0.00079439 0
0.00040657 0.0005021 0.00064226 0.00081984 0
0.00035784 0.000455 0.00060828 0.00082334 0
0.00026866 0.00035628 0.00051263 0.0007824 0
0.00014262 0.00019789 0.00031956 0.0006113 0
0 0 0 0 0 3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) Columns 7 through 11
0 0 0 0 0
2.2118e-005 2.055e-005 1.707e-005 1.0851e-005 0
3.8698e-005 3.5581e-005 2.8951e-005 1.7698e-005 0
5.0523e-005 4.6157e-005 3.7128e-005 2.2265e-005 0
5.8126e-005 5.2942e-005 4.2365e-005 2.5203e-005 0
6.171e-005 5.6197e-005 4.4952e-005 2.672e-005 0
6.1116e-005 5.5803e-005 4.4817e-005 2.6785e-005 0
5.5803e-005 5.1239e-005 4.1529e-005 2.5153e-005 0
4.4817e-005 4.1529e-005 3.42e-005 2.126e-005 0
2.6785e-005 2.5153e-005 2.126e-005 1.3869e-005 0
0 0 0 0 0。