ADI交替方向隐格式求解二维抛物方程含matlab程序

1、ADI法求解二维抛物方程学校：中国石油大学（华东） 学院：理学院 姓名：张道德 时间：1、ADI法介绍作为模型，考虑二维热传导方程的边值问题：（3.6.1）取空间步长，时间步长，作两族平行于坐标轴的网线：将区域分割成个小矩形。第一个ADI算法（交替方向隐格式）是Peaceman和Rachford（1955）提出的。方法：由第n层到第n+1层计算分为两步：（1） 第一步： ，构造出差分格式为：（2） 第二步：，构造出差分格式为：其中。假定第n层的已求得，则由求出，这只需按行解一些具有三对角系数矩阵的方程组；再由求出，这只需按列解一些具有三对角系数矩阵的方程组，所以计算时容易实现的。2、数值例子（

2、1）问题用ADI法求解二维抛物方程的初边值问题：已知（精确解为：）设差分解为，则边值条件为：初值条件为：取空间步长，时间步长网比。用ADI法分别计算到时间层。（2）计算过程根据边值条件：，已经知道第0列和第K列数值全为0。（1），构造出差分格式为：从而得到：，其中即按行用追赶法求解一系列下面的三对角方程组：又根据边值条件得：，解出第0行和第行。（2）第二步：，构造出差分格式为：从而得到：，其中又根据边值条件得：，从而得到：其中即按列用追赶法求解一系列下面的三对角方程组：（3） 求解结果(3.1)数值解yx1/42/43/41/40.2008998667134842/43.03768181457

3、584e-153/4-0.200899866713473（3.2）精确解yx1/42/43/41/42/41.26088801585392e-171.26088801585392e-173/4（3.3）数值解-精确解（即误差）yx1/42/43/41/4-0.00354880851443196-0.003548808514432732/43/40.003548808514439730.00354880851444026从而得到误差的范数为：1- 范数：0.233770443573713； 2-范数：；-范数：0.327253314506086（3.4）图像（3.4.1）数值解图像：（3.4.2

4、） 精确解图像：（5）主要程序（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);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);%*%绘

5、制图像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)

6、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+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)

7、*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,:); %*

8、%从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,cc,dd); end %* u2=u;end%*（5.3）“追赶法”程序%*%追赶法function x=zhuiganfa(a,b,c,d)%对角线下方的元素，个数比A少一个% %对角线元素%对角线上方的元素，个数比A少一个%d为方程常数项%用追赶法解三对角矩阵方程r=size(a);m=r(2)

9、;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)-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