ADI(交替方向隐格式)求解二维抛物方程(含m.docx
- 文档编号:27370125
- 上传时间:2023-06-30
- 格式:DOCX
- 页数:9
- 大小:524.63KB
ADI(交替方向隐格式)求解二维抛物方程(含m.docx
《ADI(交替方向隐格式)求解二维抛物方程(含m.docx》由会员分享,可在线阅读,更多相关《ADI(交替方向隐格式)求解二维抛物方程(含m.docx(9页珍藏版)》请在冰豆网上搜索。
ADI法求解二维抛物方程法求解二维抛物方程学校:
中国石油大学(华东)中国石油大学(华东)学院:
理学院理学院姓名:
张道德张道德时间:
2013.4.2711、ADIADI法介绍法介绍作为模型,考虑二维热传导方程的边值问题:
(3.6.1),0,0(,0)(,)(0,)(,)(,0,)(,)0txxyyuuuxyltuxyxyuytulytuxtuxlt=+=取空间步长1hM=,时间步长0t,作两族平行于坐标轴的网线:
0,1,jkxxjhyykhjkM=将区域0,xyl分割成2M个小矩形。
第一个ADIADI算法(交替方向隐格式)算法(交替方向隐格式)是Peaceman和Rachford(1955)提出的。
方法:
方法:
由第n层到第n+1层计算分为两步:
(1)第一步:
12,12njkxxyyu+从n-n+,求u对向后差分,u向前差分,构造出差分格式为:
1(3.6.1)11112222,1,1,1,1221222,2-22=21()nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhh+=+uuuuuuuu(+)uu
(2)第二步:
12,12njkxxyyu+从n+-n+1,求u对向前差分,u向后差分,构造出差分格式为:
2(3.6.1)1111111222,1,1,1,12212212,2-22=21()nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhh+=+uuuuuuuu(+)uu其中1211,1,1,0,1,2,()22njkMnnn+=+=+上表表示在t=t取值。
假定第n层的,njku已求得,则由1(3.6.1)求出12,njku+,这只需按行(1,1)jM=解一些具有三对角系数矩阵的方程组;再由2(3.6.1)求出1,njku+,这只需按列(1,1)kM=解一些具有三对角系数矩阵的方程组,所以计算时容易实现的。
22、数值例子、数值例子(11)问题问题用ADI法求解二维抛物方程的初边值问题:
21(),(,)(0,1)(0,1),0,4(0,)(1,)0,01,0,(,0,)(,1,)0,01,0(,0)sincos.xxyyyyuuuxyGttuytuytytuxtuxtxtuxyxy=+=已知(精确解为:
2(,)sincosexp()8uxytxyt=)设(0,1,),(0,1,),(0,1,)jknxjhjJykhkKtnnN=差分解为,njku,则边值条件为:
0,0,1,1,0,0,1,0,1,nnkJknnnnjjjKjKuukKuuuujJ=初值条件为:
0,sincosjkjkuxy=取空间步长12140hhh=,时间步长11600=网比21rh=。
用ADI法分别计算到时间层1t=。
(22)计算过程计算过程从n到n+1时,根据边值条件:
0,0,0,1,nnkJkuukK=,已经知道第0列和第K列数值全为0。
(1)12,12njkxxyyu+从n-n+,求u对向后差分,u向前差分,构造出差分格式为:
11112222,1,1,1,1221222,2-221=1621()16nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhh+=+uuuuuuuu(+)uu从而得到:
1112221,1,1,1111111
(1)
(1)321632321632nnnnnnjkjkjkjkjkjkrrrrrr+=+uuuuuu,其中1,2,1,1,2,1jJkK=即按行用追赶法求解一系列下面的三对角方程组:
121,122,123,123,122,121,
(1)
(1)111163211113216321111321632111132163211113216321113216nknknknJknJknJkJJrrrrrrrrrrrrrrrr+uuuuuu123321
(1)1
(1)1JJJJJffffff=又根据边值条件得:
0,1,1,0,1,nnnnjjjKjKuuuujJ=,解出第0行,0nju和第K行,(0,1,)njKujJ=。
(2)第二步:
12,12njkxxyyu+从n+-n+1,求u对向前差分,u向后差分,构造出差分格式为:
1111111222,1,1,1,12212212,2-221=1621()16nnnnnnnnjkjkjkjkjkjkjkjknnxjkyjkhhh+=+uuuuuuuu(+)uu从而得到:
111111222,1,11,1,111111
(1)
(1)321632321632nnnnnnjkjkjkjkjkjkrrrrrr+=+uuuuuu,其中1,2,1,1,2,1jJkK=又根据边值条件得:
0,1,1,0,1,nnnnjjjKjKuuuujJ=,从而得到:
0,1,1,00nnjjnnjKjKuuuu=+=其中(0,1,)jJ=即按列用追赶法求解一系列下面的三对角方程组:
1,01,11,21,31,31,21111113216321111321632111132163211113216321111321632111132163211njnjnjnjnjKnjKKKrrrrrrrrrrrrrrrrrr+uuuuuuu12343211,111,1KKnKjKKKnjKKffffffff+=u(3)求解结果求解结果(3.13.1)数值解数值解yx1/42/43/41/40.1420576580925780.2008998667134840.1420576580925782/42.16292994886484e-153.03768181457584e-152.12330312762773e-153/4-0.142057658092571-0.200899866713473-0.142057658092570(3.3.22)精确解)精确解yx1/42/43/41/40.1456064666070100.2059186398448590.1456064666070102/41.26088801585392e-171.78316493265431e-171.26088801585392e-173/4-0.145606466607010-0.205918639844859-0.145606466607010(3.33.3)数值解)数值解-精确解(即误差)精确解(即误差)yx1/42/43/41/4-0.00354880851443196-0.00501877313137564-0.003548808514432732/42.15032106870631e-153.01985016524929e-152.11069424746919e-153/40.003548808514439730.005018773131386520.00354880851444026从而得到误差的范数为:
1-范数:
0.233770443573713;2-范数:
0.196807761631447;-范数:
0.327253314506086(3.43.4)图像)图像(3.4.1)数值解图像:
(3.4.2)精确解图像:
00.20.40.60.8100.51-0.4-0.200.20.4x轴图一、数值解图像y轴t轴00.20.40.60.8100.51-0.4-0.200.20.4x轴图二、精确解图像y轴t轴(55)主要程序)主要程序(5.1)主程序%*%main_chapter主函数%信息10-2张道德%学号:
10071223clccleara=0;b=1;%x取值范围c=0;d=1;%y取值范围tfinal=1;%最终时刻t=1/1600;%时间步长;h=1/40;%空间步长r=t/h2;%网比x=a:
h:
b;y=c:
h:
d;%*%精确解m=40;u1=zeros(m+1,m+1);fori=1:
m+1,forj=1:
m+1u1(j,i)=uexact(x(i),y(j),1);endend%数值解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)%*functionu=ADI(a,b,c,d,t,h,tfinal)%(a,b)x取值范围%(c,d)y取值范围%tfinal最终时刻%t时间步长;%h空间步长r=t/h2;%网比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+1u(j,i)=uexact(x(i),y(j),0);endend%*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:
mfork=2:
md(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);endu2(1,:
)=u2(2,:
);u2(m+1,:
)=u2(m,:
);%*%从n-n+1,u_xx向前差分,u_yy向后差分fork=2:
mdd
(1)=0;dd(m+1)=0;forj=2:
mdd(j)=1/32*r*(u2(j+1,k)-2*u2(j,k)+u2(j-1,k)+u2(j,k);endu(:
k)=zhuiganfa(aa,bb,cc,dd);end%*u2=u;end%*(5.3)“追赶法”程序%*%追赶法functionx=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:
nl(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:
ny(i)=d(i)-l(i-1)*y(i-1);end%求解Ux=y,赶的过程x(n)=y(n)/u(n);fori=n-1:
-1:
1x(i)=y(i)/u(i);x(i)=(y(i)-c(i)*x(i+1)/u(i);end%*(5.4)精确解函数%t时刻,u的取值;functionf=uexact(x,y,t)f=sin(x*pi)*cos(y*pi)*exp(-pi*pi/8*t);%*
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- ADI 交替 方向 格式 求解 二维 方程