ADI交替方向隐格式求解二维抛物方程含matlab程序.docx
- 文档编号:11093844
- 上传时间:2023-05-29
- 格式:DOCX
- 页数:11
- 大小:119.21KB
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.142057658092578
0.202199866713484
0.142057658092578
2/4
2.16292994886484e-15
3.03768181457584e-15
2.12330312762773e-15
3/4
-0.142057658092571
-0.202199866713473
-0.142057658092570
〔3.2〕准确解
yx
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〕数值解-准确解〔即误差〕
yx
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);
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_{**}+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_{**}向后差分,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_{**}向前差分,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);
%********************************************************************
教育之通病是教用脑的人不用手,不教用手的人用脑,所以一无所能。
教育革命的对策是手脑联盟,结果是手与脑的力量都可以大到不可思议。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- ADI 交替 方向 格式 求解 二维 方程 matlab 程序