ADI交替方向隐格式.docxVIP

PAGE 1 ADI算法的MATLAB编程应用实例 胡坤1 ，任兰兰2 1．ADI算法的具体描述 ADI算法又称交替方向隐格式，该算法主要考虑二维热传导方程的边值问题，模型如下： 在上述模型中，取空间步长，时间步长，作两族平行于坐标轴的网线：将区域分割成个小矩形。 具体步骤是将第n层到第n+1层计算分为两步： 第一步： ，构造出差分格式为： 第二步：，构造出差分格式为： 其中 假定第n层的已求得，则由上述第一步可求出，这只需按行解一些具有三对角系数矩阵的方程组；再由第二步求出，这只需按列解一些具有三对角系数矩阵的方程组。 2．以ADI算法分析具体实例 （1）考察例子 上述方程精确解为：） （2）分析计算过程 首先设差分解为，则边值条件为： 初值条件为：. 取空间步长，时间步长网比，用ADI法分别计算到时间层. 根据边值条件：，已经知道第0列和第K列数值全为0，故： 第一步： ，构造出差分格式为： 从而得到： ，其中 即按行用追赶法求解一系列下面的三对角方程组： 又根据边值条件得：，解出第0行和第行. 第二步： 从而得到： 其中 又根据边值条件得：，从而得到： 其中 再按列用追赶法求解一系列下面的三对角方程组： 从而得到新的时间层的数值解. 3.MATLAB编程实现上述实例 clear clc 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); for i=1:m+1, for j=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); % 编写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空间步长 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); for i=1:m+1, for j=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; 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); 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