Numerical Solution of Partial Differential Equations

ADI method

交替方向隐式格式的每一次步进需要在\(x\)方向和\(y\)方向上分别步进一次. 我们假定\(x\)方向上的隐士格式的微分矩阵分别是\(A_x\)和\(B_x\), \(y\)方向上的隐士格式的微分矩阵分别是\(A_y\)和\(B_y\). 假定xx和yy是网格节点, \(U(i,j)\)是对应于xx(i,j),yy(i,j)的值(不包括边界), 则

\[ (1+.5\mu_y\delta_y^2)U = B_yU \]

\[ (1+.5\mu_x\delta_x^2)U = (B_xU’)’ \]

为了说明为什么\(x\)方向上的微分算子在用矩阵表示的时候需要进行转置, 我们来看下面的一段程序

x = 0:0.2:1;
y = 0:0.2:1;
[xx,yy] = meshgrid(x,y);

它的输出为

xx =

         0    0.2000    0.4000    0.6000    0.8000    1.0000
         0    0.2000    0.4000    0.6000    0.8000    1.0000
         0    0.2000    0.4000    0.6000    0.8000    1.0000
         0    0.2000    0.4000    0.6000    0.8000    1.0000
         0    0.2000    0.4000    0.6000    0.8000    1.0000
         0    0.2000    0.4000    0.6000    0.8000    1.0000


yy =

         0         0         0         0         0         0
    0.2000    0.2000    0.2000    0.2000    0.2000    0.2000
    0.4000    0.4000    0.4000    0.4000    0.4000    0.4000
    0.6000    0.6000    0.6000    0.6000    0.6000    0.6000
    0.8000    0.8000    0.8000    0.8000    0.8000    0.8000
    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000

而矩阵\(U\)每个位置对应的横坐标和纵坐标就是xx以及yy对应位置的元素. 因此, 在\(y\)方向上做差分, 只需要使用\(B_yU\)即可; 而在\(x\)方向上做差分, 需要进行转置. 最后，放上一个M型的source函数的代码

x1 = 0:0.01:1;
y1 = 0:0.01:1;
[xx1,yy1] = meshgrid(x1,y1);
load('M.mat')
u0 = @(x,y) interp2(xx1,yy1,M,x,y);

其中M.mat由蔡大亨提供，可以点此下载