ADI交替方向隐格式求解二维抛物方程含matlab程序Word格式文档下载.docx
《ADI交替方向隐格式求解二维抛物方程含matlab程序Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《ADI交替方向隐格式求解二维抛物方程含matlab程序Word格式文档下载.docx（11页珍藏版）》请在冰豆网上搜索。
（1）问题
用ADI法求解二维抛物方程的初边值问题：
已知（精确解为：
）
设差分解为，则边值条件为：
初值条件为：
取空间步长，时间步长网比。
用ADI法分别计算到时间层。
（2）计算过程
根据边值条件：
，已经知道第0列和第K列数值全为0。
（1），构造出差分格式为：
从而得到：
，其中
即按行用追赶法求解一系列下面的三对角方程组：
又根据边值条件得：
，解出第0行和第行。
，
其中
即按列用追赶法求解一系列下面的三对角方程组：
（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
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）;
升级会员 