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

ADI法求解二维抛物方程
学校：中国石油大学（华东） 学院：理学院 姓名：张道德 时间：2013.4.27
1、ADI法介绍
作为模型，考虑二维热传导方程的边值问题：
（3.6.1）
取空间步长，时间步长，作两族平行于坐标轴的网线：将区域分割成个小矩形。第一个ADI算法（交替方向隐格式）是Peaceman和Rachford（1955）提出的。
方法：
由第n层到第n+1层计算分为两步：
（1） 第一步： ，构造出差分格式为：

（2） 第二步：，构造出差分格式为：

其中。假定第n层的已求得，则由求出，这只需按行解一些具有三对角系数矩阵的方程组；再由求出，这只需按列解一些具有三对角系数矩阵的方程组，所以计算时容易实现的。
2、数值例子
（1）问题
用ADI法求解二维抛物方程的初边值问题：

已知（精确解为：）
设差分解为，则边值条件为：
初值条件为：
取空间步长，时间步长网比。用ADI法分别计算到时间层。
（2）计算过程

根据边值条件：，已经知道第0列和第K列数值全为0。
（1），构造出差分格式为：

从而得到：
，其中
即按行用追赶法求解一系列下面的三对角方程组：

又根据边值条件得：，解出第0行和第行。
（2）第二步：，构造出差分格式为：

从而得到：
，其中
又根据边值条件得：，
从而得到：
其中
即按列用追赶法求解一系列下面的三对角方程组：

（3） 求解结果
(3.1)数值解
y x
1/4
2/4
3/4
1/4
0.142057658092578
0.200899866713484
0.142057658092578
2/4
2.16292994886484e-15
3.03768181457584e-15
2.12330312762773e-15
3/4
-0.142057658092571
-0.200899866713473
-0.142057658092570

（3.2）精确解
y x
1/4
2/4
3/4
1/4
0.145606466607010
0.205918639844859
0.145606466607010
2/4
1.26088801585392e-17
1.78316493265431e-17
1.26088801585392e-17
3/4
-0.145606466607010
-0.205918639844859
-0.145606466607010

（3.3）数值解-精确解（即误差）
y x
1/4
2/4
3/4
1/4
-0.00354880851443196
-0.00501877313137564
-0.00354880851443273
2/4
2.15032106870631e-15
3.01985016524929e-15
2.11069424746919e-15
3/4
0.00354880851443973
0.00501877313138652
0.00354880851444026

从而得到误差的范数为：

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);
%********************************************************************