ADI(交替方向隐格式)求解二维抛物方程(含matlab程序)

ADI 法求解二维抛物方程

学校：中国石油大学（华东） 学院：理学院 姓名：张道德 时间：2013.4.27

1、ADI 法介绍

作为模型，考虑二维热传导方程的边值问题：

（3.6.1）,0,,0(,,0)(,)(0,,)(,,)(,0,)(,,)0t xx yy u u u x y l t u x y x y u y t u l y t u x t u x l t ?=+<<>??

=??====?

取空间步长1

h M

=，时间步长0t >，作两族平行于坐标轴的网线：

,,,0,1,

,,j k x x jh y y kh j k M =====将区域0,x y l ≤≤分割成2M 个小矩形。第

一个ADI 算法（交替方向隐格式）是Peaceman 和Rachford （1955）提出的。 方法：

由第n 层到第n+1层计算分为两步：

（1） 第一步： 12,1

2

n j k xx yy u +从n->n+，求u 对向后差分,u 向前差分，构造出差分格

式为：

1

(3.6.1)11112222,,1,,1,,1,,1

2

2

1

222,,2-22=2

1()n n n n n n n n j k

j k

j k

j k j k

j k j k j k n n x j k y j k h

h

h

τ

δδ+

+

+

+

+-+-+-+-+=+u

u u

u

u

u u u （

+

）

u u

（2） 第二步：12,1

2

n j k xx yy u +从n+->n+1，求u 对向前差分,u 向后差分，构造出差分格

式为：

2

(3.6.1)1111111222,,1,,1,,1,,1

2

2

1

2212,,2-22=2

1()n n n n n n n n j k

j k

j k

j k j k

j k j k j k n n x j k y j k h

h h

τ

δδ+

+

+

+++++-+-++-+-+=+u

u u

u

u

u u u （

+

）

u u

其中12

11

,1,,1,0,1,2,

,()22n j k M n n n τ+=-=+=+上表表示在t=t 取值。

假定第n 层的,n j k

u

已求得，则由1(3.6.1)求出12

,n j k u +，这只需按行

(1,

,1)j M =-解一些具有三对角系数矩阵的方程组；再由2(3.6.1)求出

1,n j k u +，这只需按列(1,

,1)k M =-解一些具有三对角系数矩阵的方程组，

所以计算时容易实现的。

2、数值例子

（1）问题

用ADI 法求解二维抛物方程的初边值问题：

21

(),(,)(0,1)(0,1),0,4(0,,)(1,,)0,01,0,(,0,)(,1,)0,01,0(,,0)sin cos .

xx yy y y u u u x y G t t u y t u y t y t u x t u x t x t u x y x y ππ??=+∈=?>????==<<>??==<<>?=?? 已知（精确解为：2

(,,)sin cos exp()8

u x y t x y t πππ=-）

设(0,1,

,),(0,1,

,),(0,1,

,)j k n x jh j J y kh k K t n n N τ======差分解为,n j k u ，

则边值条件为：0,,,0,1,1,0,0,1,,,,0,1,

,n n

k J k n

n n n

j j j K j K u u k K

u u u u j J

-?===??===??

初值条件为：0,sin cos j k j k u x y ππ=

取空间步长12140h h h ===，时间步长11600τ=网比21r h τ==。用ADI 法分别计算到时间层1t =。

（2）计算过程 从n 到n+1时，

根据边值条件：0,,0,0,1

,,n n

k J k u u k K ===，

已经知道第0列和第K 列数值全为0。 （1）1

2,1

2

n j k xx yy u +从n->n+，求u 对向后差分,u 向前差分，构造出差分格式为：

11112222,,1,,1,,1,,1

2

2

1222

,,2-221=16

2

1()16n n n n n n n n j k

j k

j k

j k j k

j k j k j k n n x j k y j k h

h

h

τ

δδ+++++-+-+-+-+=+u u u u u

u u u （+

）

u u

从而得到：

111222

1,,1,,1,,1111111(1)(1)321632321632

n n n n n n j k j k j k j k j k j k r r r r r r ++++-+--++-=+--u u u u u u ，其中

1,2,,1,1,2,,1j J k K =-=-

即按行用追赶法求解一系列下面的三对角方程组：

121,122,123,123,

122,121,(1)(1)111163211113216321111321632111132163211113216321113216n k n k n k n J k n J k n J k J J r r r r r r r

r r r r r r r r r ++++-+-+--?-???+-????

????-+-?????

????-+-?????????-+-??????-+-??????-+????

?u u u u u u 123321(1)1(1)1

J J J J J f f f f f f ----?-???

???????????????????=????????????????????????????

????

? 又根据边值条件得：,0,1,1,,,0,1,,n

n

n

n

j j j K j K u u u u j J -===，解出第0行,0n

j u 和第

K 行,,(0,1

,,)n j K u j J =。

（2）第二步：1

2,1

2

n j k xx yy u +从n+->n+1，求u 对向前差分,u 向后差分，构造出差分格

式为：

1111111222,,1,,1,,1,,1

2

2

12

212,,2-221=16

2

1()16n n n n n n n n j k

j k

j k

j k j k

j k j k j k n n x j k y j k h

h

h

τ

δδ+

+

+

+++++-+-++-+-+=+u

u u u

u

u u u （+

）

u u

从而得到：

111111222

,1,,11,,1,111111(1)(1)321632321632

n n n n n n j k j k j k j k j k j k r r r r r r +++++++-+--++-=+--u u u u u u ，其中1,2,,1,

1,2,j J k K =-=-

又根据边值条件得：,0,1,1,,,0,1,,n

n

n

n

j j j K j K u u u u j J -===，

从而得到：

,0,1,1,0

n n

j j n n

j K j K u u u u -?-=??-+=??其中(0,1,,)j J = 即按列用追赶法求解一系列下面的三对角方程组：

1

,01,11,21,31,31,2111111321632111

1321632111

13216321111321632111

1321632111132163211n j n j n j n j n j K n j K K K

r r r r r r r r r r r r

r r r r r r +++++-+-?-??

????

-+-??

??-+-????

??-+-????

??????-+-??

??-+-????

??

-+-?

??

?-???

?u u u u u u u 12343211,11

1,1

K K n K j K K K n j K K f f f f f f f f --+--?+???????????????????????????

????=??

????????????????????????????

u （3） 求解结果

从而得到误差的范数为：

1- 范数：0.233770443573713； 2-范数：0.196807761631447； ∞-范数：0.327253314506086

（3.4）图像

（3.4.1）数值解图像：

（3.4.2）精确解图像：

（5）主要程序

（5.1）主程序

%*************************************************************

%main_chapter主函数

%信息10-2——张道德

%学号：10071223

clc

clear

a = 0; b=1; %x取值范围

c=0; d=1; %y取值范围

tfinal = 1; %最终时刻

t=1/1600;%时间步长；

h=1/40;%空间步长

r=t/h^2;%网比

x=a:h:b;

y=c:h:d;

%**************************************************************

%精确解

m=40;

u1=zeros(m+1,m+1);

for i=1:m+1,

for j=1:m+1

u1(j,i) = uexact(x(i),y(j),1);

end

end

%数值解

u=ADI(a,b,c,d,t,h,tfinal);

%***************************************************************** %绘制图像

figure(1) ;mesh(x,y,u1)

figure(2); mesh(x,y,u)

%误差分析

error=u-u1;

norm1=norm(error,1);

norm2=norm(error,2);

norm00=norm(error,inf);

%*****************************************************************（5.2）ADI函数

%****************************************************************

% 用ADI法求解二维抛物方程的初边值问题

% u_t = 1/16（u_{xx} + u_{yy}）（0,1）*（0,1）

% 精确解: u(t,x,y) = sin(pi*x) sin(pi*y)exp(-pi*pi*t/8) %****************************************************************

function [u]=ADI(a,b,c,d,t,h,tfinal )

%(a , b) x取值范围

%(c, d) y取值范围

%tfinal最终时刻

%t时间步长；

%h空间步长

r=t/h^2;%网比

m=(b-a)/h;%

n=tfinal/t; %

x=a:h:b;

y=c:h:d;

%****************************************************************** %初始条件

u=zeros(m+1,m+1);

for i=1:m+1,

for j=1:m+1

u(j,i) = uexact(x(i),y(j),0);

end

end

%****************************************************************** u2=zeros(m+1,m+1);

a=-1/32*r*ones(1,m-2);

b=(1+r/16)*ones(1,m-1);

aa=-1/32*r*ones(1,m);

cc=aa;aa(m)=-1;cc(1)=-1;

bb=(1+r/16)*ones(1,m+1);

bb(1)=1;bb(m+1)=1;

for i=1:n

%************************************************************** %从n->n+1/2,u_{xx}向后差分，u_{yy}向前差分

for j=2:m

for k=2:m

d(k-1)=1/32*r*(u(j,k+1)-2*u(j,k)+u(j,k-1))+u(j,k);

end

% 修正第一项与最后一项,但由于第一项与最后一项均为零，可以省略

%d(1)=d(1)+u1(j,1);d(m-1)=d(m-1)+u1(j,m+1);

u2(j,2:m)=zhuiganfa(a,b,a,d);

end

u2(1,:)=u2(2,:);

u2(m+1,:)=u2(m,:);

%************************************************************** %从n->n+1,u_{xx}向前差分，u_{yy}向后差分

for k=2:m

dd(1)=0;dd(m+1)=0;

for j=2:m

dd(j)=1/32*r*(u2(j+1,k)-2*u2(j,k)+u2(j-1,k))+u2(j,k);

end

u(:,k)=zhuiganfa(aa,bb,cc,dd);

end

%**************************************************************** u2=u;

end

%********************************************************************（5.3）“追赶法”程序

%******************************************************************** %追赶法

function [x]=zhuiganfa(a,b,c,d)

%对角线下方的元素，个数比A少一个

% %对角线元素

%对角线上方的元素，个数比A少一个

%d为方程常数项

%用追赶法解三对角矩阵方程

r=size(a);

m=r(2);

r=size(b);

n=r(2);

if size(a)~=size(c)|m~=n-1|size(b)~=size(d)

error('变量不匹配,检查变量输入情况!');

end

%%

%LU分解

u(1)=b(1);

for i=2:n

l(i-1)=a(i-1)/u(i-1);

u(i)=b(i)-l(i-1)*c(i-1);

v(i-1)=(b(i)-u(i))/l(i-1);

end

%追赶法实现

%%

%求解Ly=d，"追"的过程

y(1)=d(1);

for i=2:n

y(i)=d(i)-l(i-1)*y(i-1);

end

%%

%求解Ux=y,"赶"的过程

x(n)=y(n)/u(n);

for i=n-1:-1:1

x(i)=y(i)/u(i);

x(i)=(y(i)-c(i)*x(i+1))/u(i);

end

%********************************************************************（5.4）精确解函数

%t时刻，u的取值；

function [ f]=uexact(x,y,t)

f=sin(x*pi)*cos(y*pi)*exp(-pi*pi/8*t);

%********************************************************************