ADI交替方向隐格式求解二维抛物方程含matlab程序.docx
《ADI交替方向隐格式求解二维抛物方程含matlab程序.docx》由会员分享，可在线阅读，更多相关《ADI交替方向隐格式求解二维抛物方程含matlab程序.docx（11页珍藏版）》请在冰豆网上搜索。
ADI交替方向隐格式求解二维抛物方程含matlab程序
ADI法求解二维抛物方程
学校：
中国石油大学〔华东〕学院：
理学院：
道德时间：
2021.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）数值解
yx
1/4
2/4
3/4
1/4
0.2578
0.3484
0.2578
2/4
2.484e-15
3.584e-15
2.773e-15
3/4
-0.2571
-0.3473
-0.2570
〔3.2〕准确解
yx
1/4
2/4
3/4
1/4
0.7010
0.4859
0.7010
2/4
1.392e-17
1.431e-17
1.392e-17
3/4
-0.7010
-0.4859
-0.7010
〔3.3〕数值解-准确解〔即误差〕
yx
1/4
2/4
3/4
1/4
-0.
-0.
-0.
2/4
2.631e-15
3.929e-15
2.919e-15
3/4
0.
0.
0.
从而得到误差的数为：
1-数：
0.3713；2-数：
0.1447；
∞-数：
0.6086
〔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）;
fori=1:
m+1,
forj=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）;
fori=1:
m+1,
forj=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;
fori=1:
n
%**************************************************************
%从n->n+1/2,u_{xx}向后差分，u_{yy}向前差分
forj=2:
m
fork=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}向后差分
fork=2:
m
dd
（1）=0;dd（m+1）=0;
forj=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）;
ifsize（a）~=size（c）|m~=n-1|size（b）~=size（d）
error（'变量不匹配,检查变量输入情况!
'）;
end
%%
%LU分解
u
（1）=b
（1）;
fori=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）;
fori=2:
n
y（i）=d（i）-l（i-1）*y（i-1）;
end
%%
%求解Ux=y,"赶"的过程
x（n）=y（n）/u（n）;
fori=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）;
%********************************************************************
升级会员 