ADI(交替方向隐格式)

上传人:夏** 文档编号:496301270 上传时间:2024-02-25 格式:DOCX 页数:8 大小:335.02KB
返回 下载 相关 举报
ADI(交替方向隐格式)_第1页
第1页 / 共8页
ADI(交替方向隐格式)_第2页
第2页 / 共8页
ADI(交替方向隐格式)_第3页
第3页 / 共8页
ADI(交替方向隐格式)_第4页
第4页 / 共8页
ADI(交替方向隐格式)_第5页
第5页 / 共8页
点击查看更多>>
资源描述

《ADI(交替方向隐格式)》由会员分享,可在线阅读,更多相关《ADI(交替方向隐格式)(8页珍藏版)》请在金锄头文库上搜索。

1、精选优质文档-倾情为你奉上ADI算法的MATLAB编程应用实例胡坤1 ,任兰兰21ADI算法的具体描述ADI算法又称交替方向隐格式,该算法主要考虑二维热传导方程的边值问题,模型如下:在上述模型中,取空间步长,时间步长,作两族平行于坐标轴的网线:将区域分割成个小矩形。具体步骤是将第n层到第n+1层计算分为两步:(1) 第一步: ,构造出差分格式为:(2) 第二步:,构造出差分格式为:其中假定第n层的已求得,则由上述第一步可求出,这只需按行解一些具有三对角系数矩阵的方程组;再由第二步求出,这只需按列解一些具有三对角系数矩阵的方程组。 2以ADI算法分析具体实例(1)考察例子上述方程精确解为:)(2

2、)分析计算过程首先设差分解为,则边值条件为:初值条件为:.取空间步长,时间步长网比,用ADI法分别计算到时间层.根据边值条件:,已经知道第0列和第K列数值全为0,故:第一步:,构造出差分格式为:从而得到:,其中即按行用追赶法求解一系列下面的三对角方程组:又根据边值条件得:,解出第0行和第行.第二步: 从而得到:其中又根据边值条件得:,从而得到:其中再按列用追赶法求解一系列下面的三对角方程组:从而得到新的时间层的数值解.3.MATLAB编程实现上述实例clearclca = 0; b=1; %x取值范围c=0; d=1; %y取值范围tfinal = 1; %最终时刻t=1/1600;%时间步长

3、;h=1/40;%空间步长r=t/h2;%网比x=a:h:b;y=c:h:d;%-%精确解m=40;u1=zeros(m+1,m+1);for i=1:m+1, for j=1:m+1 u1(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);%-编写ADI函数

4、文件% 用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/h2;%网比m=(b-a)/h;%n=tfinal/t; %x=a:h:b;y=c:h:d;%-%初始条件u=zeros(m+1,m+1);for i=1:m+1, for j=1:m

5、+1 u(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;for i=1:n %- %从n-n+1/2,u_xx向后差分,u_yy向前差分 for j=2:m for k=2:m d(k-1)=1/32*r*(u(j,k+1)-2*u(j,k)+u(j,k-1)+u(j,k)

6、; 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向后差分 for k=2:m dd(1)=0;dd(m+1)=0; for j=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

7、,cc,dd); end %- u2=u;end%- “追赶法”解三对角线性方程函数文件%-%追赶法function x=zhuiganfa(a,b,c,d)%对角线下方的元素,个数比A少一个% %对角线元素%对角线上方的元素,个数比A少一个%d为方程常数项%用追赶法解三对角矩阵方程r=size(a);m=r(2);r=size(b);n=r(2);if size(a)=size(c)|m=n-1|size(b)=size(d) error(变量不匹配,检查变量输入情况!);end%LU分解u(1)=b(1);for i=2:n l(i-1)=a(i-1)/u(i-1); u(i)=b(i)-

8、l(i-1)*c(i-1); v(i-1)=(b(i)-u(i)/l(i-1); end%追赶法实现%求解Ly=d,追的过程y(1)=d(1);for i=2:n y(i)=d(i)-l(i-1)*y(i-1);end%求解Ux=y,赶的过程x(n)=y(n)/u(n);for i=n-1:-1:1 x(i)=y(i)/u(i); x(i)=(y(i)-c(i)*x(i+1)/u(i);end%-精确解函数文件%t时刻,u的取值;function f=uexact(x,y,t)f=sin(x*pi)*cos(y*pi)*exp(-pi*pi/8*t);%-图1 以方程精确解所绘出的网格图图2 以ADI算法得出方程数值解绘出的网格图专心-专注-专业

展开阅读全文
相关资源
相关搜索

当前位置:首页 > 办公文档 > 教学/培训

电脑版 |金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号