ADI交替方向隐格式求解二维抛物方程含matlab程序Word格式文档下载.docx
- 文档编号:15019901
- 上传时间:2022-10-26
- 格式:DOCX
- 页数:11
- 大小:223.39KB
ADI交替方向隐格式求解二维抛物方程含matlab程序Word格式文档下载.docx
《ADI交替方向隐格式求解二维抛物方程含matlab程序Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《ADI交替方向隐格式求解二维抛物方程含matlab程序Word格式文档下载.docx(11页珍藏版)》请在冰豆网上搜索。
〔1〕问题
用ADI法求解二维抛物方程的初边值问题:
〔准确解为:
〕
设差分解为,那么边值条件为:
初值条件为:
取空间步长,时间步长网比。
用ADI法分别计算到时间层。
〔2〕计算过程
根据边值条件:
,已经知道第0列和第K列数值全为0。
〔1〕,构造出差分格式为:
从而得到:
,其中
即按行用追赶法求解一系列下面的三对角方程组:
又根据边值条件得:
,解出第0行和第行。
〔2〕第二步:
,
其中
即按列用追赶法求解一系列下面的三对角方程组:
(3)求解结果
(3.1)数值解
yx
1/4
2/4
3/4
0.2578
0.3484
2.484e-15
3.584e-15
2.773e-15
-0.2571
-0.3473
-0.2570
〔3.2〕准确解
0.7010
0.4859
1.392e-17
1.431e-17
-0.7010
-0.4859
〔3.3〕数值解-准确解〔即误差〕
-0.
2.631e-15
3.929e-15
2.919e-15
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:
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
%数值解
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空间步长
m=(b-a)/h;
%
n=tfinal/t;
%
%******************************************************************
%初始条件
u=zeros(m+1,m+1);
u(j,i)=uexact(x(i),y(j),0);
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;
n
%从n->
n+1/2,u_{xx}向后差分,u_{yy}向前差分
forj=2:
m
fork=2:
d(k-1)=1/32*r*(u(j,k+1)-2*u(j,k)+u(j,k-1))+u(j,k);
%修正第一项与最后一项,但由于第一项与最后一项均为零,可以省略
%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);
u2(1,:
)=u2(2,:
);
u2(m+1,:
)=u2(m,:
n+1,u_{xx}向前差分,u_{yy}向后差分
dd
(1)=0;
dd(m+1)=0;
dd(j)=1/32*r*(u2(j+1,k)-2*u2(j,k)+u2(j-1,k))+u2(j,k);
u(:
k)=zhuiganfa(aa,bb,cc,dd);
u2=u;
%********************************************************************
〔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('
变量不匹配,检查变量输入情况!
'
%%
%LU分解
u
(1)=b
(1);
fori=2:
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);
%追赶法实现
%求解Ly=d,"
追"
的过程
y
(1)=d
(1);
y(i)=d(i)-l(i-1)*y(i-1);
%求解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);
〔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 程序