《ADI(交替方向隐格式)4521.pdf》由会员分享，可在线阅读，更多相关《ADI(交替方向隐格式)4521.pdf（8页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、ADI 算法的 MATLAB 编程应用实例 胡坤1，任兰兰2 1ADI 算法的具体描述 ADI 算法又称交替方向隐格式，该算法主要考虑二维热传导方程的边值问题，模型如下：,0,0(,0)(,)(0,)(,)(,0,)(,)0txxyyuuux yl tu x yx yuy tu l y tu xtu x l t 在上述模型中，取空间步长1hM，时间步长0，作两族平行于坐标轴的网线：,0,1,jkxxjh yykh j kM将区域0,x yl分割成2M个小矩形。具体步骤是将第 n 层到第 n+1 层计算分为两步：（1）第一步：12,12nj kxxyyu从n-n+，求u 对向后差分,u 向前差分
2、，构造出差分格式为：11112222,1,1,1,1221222,2-22=21=()nnnnnnnnj kj kjkj kjkj kj kj knnxj kyj khhhuuuuuuuu（+）uu（2）第二步：12,12nj kxxyyu从n+-n+1，求u 对向前差分,u 向后差分，构造出差分格式为：1111+11112222,1,1,1,12212212,2-22=21=()nnnnnnnnj kj kjkj kjkj kj kj knnxj kyj khhhuuuuuuuu（+）uu 其中1211,1,1,0,1,2,()22nj kMnnn其中表示在t=t取值 假 定 第 n 层 的
3、,nj ku已 求 得，则 由 上 述第 一 步可求 出12,nj ku，这 只 需 按行(1,1)jM解一些具有三对角系数矩阵的方程组；再由第二步求出1,nj ku，这只需按列(1,1)kM解一些具有三对角系数矩阵的方程组。2以 ADI 算法分析具体实例（1）考察例子 21(),(,)(0,1)(0,1),0,4(0,)(1,)0,01,0,(,0,)(,1,)0,01,0(,0)sincos.xxyyyyuuux yGttuy tuy tytuxtuxtxtu x yxy 上述方程精确解为：2(,)exp()sincos8u x y ttxy）（2）分析计算过程 首先设(0,1,),(0,
4、1,),(0,1,)jknxjh jJykh kKtnnN差分解为,nj ku，则边值条件为：0,0,1,1,0,0,1,0,1,nnkJ knnnnjjj Kj KuukKuuuujJ 初值条件为：0,sincosj kjkuxy.取空间步长121 40hhh，时间步长1 1600网比21rh，用 ADI法分别计算到时间层1t.从n 到n+1时，根据边值条件：0,0,0,1,nnkJ kuukK，已经知道第 0 列和第 K 列数值全为 0，故：第一步：12,12nj kxxyyu从n-n+，求u 对向后差分,u 向前差分，构造出差分格式为：11112222,1,1,1,1221222,2-2
5、21=1621=()16nnnnnnnnj kj kjkj kjkj kj kj knnxj kyj khhhuuuuuuuu（+）uu 从而得到：1112221,1,1,1111111(1)(1)321632321632nnnnnnjkj kjkj kj kj krrrrrruuuuuu，其 中1,2,1,1,2,1jJkK 即按行用追赶法求解一系列下面的三对角方程组：121,122,123,123,122,121,(1)(1)111163211113216321111321632111132163211113216321113216nknknknJknJknJkJJrrrrrrrrrrrr
6、rrrr uuuuuu123321(1)1(1)1JJJJJffffff 又根据边值条件得：,0,1,1,0,1,nnnnjjj Kj KuuuujJ，解出第 0 行,0nju和第K行,(0,1,)nj KujJ.第二步：12,12nj kxxyyu从n+-n+1，求u 对向前差分构造出差,u 向后差分分格式为：1111+11112222,1,1,1,12212212,2-221=1621=()16nnnnnnnnj kj kjkj kjkj kj kj knnxj kyj khhhuuuuuuuu（+）uu 从而得到：111111222,1,11,1,111111(1)(1)32163232
7、1632nnnnnnj kj kj kjkj kjkrrrrrruuuuuu 其中1,2,1,1,2,1jJkK 又根据边值条件得：,0,1,1,0,1,nnnnjjj Kj KuuuujJ，从而得到：,0,1,1,00nnjjnnj Kj Kuuuu其中(0,1,)jJ 再按列用追赶法求解一系列下面的三对角方程组：1,01,11,21,31,31,21111113216321111321632111132163211113216321111321632111132163211njnjnjnjnj Knj KK Krrrrrrrrrrrrrrrrrruuuuuuu12343211,111,1K
8、KnKj KKKnj KKffffffffu 从而得到新的时间层的数值解.编程实现上述实例 clear clc a=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);end end%数值解 u=ADI(a,b,c,d,t,h,tfinal);%-%绘制图像 figure(1);mesh(x
9、,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空间步长
10、 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);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
11、-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(m+1)=0;for j=
12、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%-“追赶法”解三对角线性方程函数文件%-%追赶法 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(
13、变量不匹配,检查变量输入情况!);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=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 以方程精确解所绘出的网格图 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轴 图2 以ADI算法得出方程数值解绘出的网格图