ADI_2D HEAT EQUATION

PDE

网格

ADI方法

步骤

x-direction

%-------------- loop in the x direction-------------------

% ---- Dirichlet boundary condition for j = n+1
for i = 1:m+1      % The upper boundary condition u(t,x,1) = 0
    u2(i,n+1) = 0;
end

for j = 1:n+1    % The left boundary condition u(t,0,y) = 0
    u2(1,j) = 0;
end

% ---- Loop for fixed y(j) 
for j = 1:n     
    A = sparse(m,m); b=zeros(m,1);
    if j == 1 %j = 1 (x-direction)---ghost
        for i=2:m+1
            b(i-1) = (-2*u1(i,j) + 2*u1(i,j+1))/hx1 + ...
		           f(t2,x(i),y(j)) + 2*u1(i,j)/dt - (2*uy(t2,x(i),y(1)))/hy; %b()-after ghost
            if i == 2 %fixed j = 1, i = 2
                b(i-1) = b(i-1) + uexact(t2,x(1),y(j))/hx1; %b()+mu*u(t,0,y)
                A(i-1,i) = -1/hx1; %-mu
            else
	            if i==m+1 %fixed j = 1,i = m+1 ---ghost
                   b(i-1) = b(i-1) + (2*ux(t2,x(m+1),y(j)))/hx; %b()+2*hx*mu*u_x(t,1,y)
                   A(i-1,i-2) =  -2/hx1; %-2*mu
                else %fixed 1< j <n+1, 2< i <m+1
	               A(i-1,i) = -1/hx1; %-mu
	               A(i-1,i-2) = -1/hx1; %-mu
                end
            end
            A(i-1,i-1) = 2/dt + 2/hx1; %1+2*mu
        end
    
        ut = A\b;      % j = 1.
        for i=1:m
            u2(i+1,j) = ut(i);
        end

    else % 1< j <n+1 (x-direction) or j = 2,...,n
        for i=2:m+1
            b(i-1) = (u1(i,j-1) -2*u1(i,j) + u1(i,j+1))/hx1 + ...
		           f(t2,x(i),y(j)) + 2*u1(i,j)/dt ; %b()
            if i == 2 %fixed 1< j <n+1, i = 2
                b(i-1) = b(i-1) + uexact(t2,x(1),y(j))/hx1; %b()+mu*u(t,0,y)
                A(i-1,i) = -1/hx1; %-mu
            else
	            if i==m+1 %fixed 1< j <n+1,i = m+1
                   b(i-1) = b(i-1) + (2*ux(t2,x(m+1),y(j)))/hx; %b()+2*hx*mu*u_x(t,1,y)
                   A(i-1,i-2) =  -2/hx1; %-2*mu
                else %fixed 1< j <n+1, 2< i <m+1
	               A(i-1,i) = -1/hx1; %-mu
	               A(i-1,i-2) = -1/hx1; %-mu
                end
            end
            A(i-1,i-1) = 2/dt + 2/hx1; %1+2*mu
        end
  
        ut = A\b;      % j = 2, ..., n.
        for i=1:m
            u2(i+1,j) = ut(i);
        end
    end
 end     % Finish x-sweep.

当 j = 1时，要先用ghost point method把U^k
{i,j−1} = U^k{i,j+1} − 2h_yu_y替换；然后fixed j = 1的情况下，当i=m+1的时候，还要用ghost point method把U^{k + 1/2}{m+2,j} = 2h_xu_x + U^{k+ 1/2}{m,j}替换了。当j = 1的时候，得到一个矩阵如下：

y-direction

%-------------- loop in y -direction --------------------------------

for i = 1:m+1      % The upper boundary condition u(t,x,1) = 0
    u1(i,n+1) = 0;
end

for j = 1:n+1    % The left boundary condition u(t,0,y) = 0
    u1(1,j) = 0;
end

for i = 2:m+1 % Loop for fixed x(i) 
    A = sparse(n,n); b=zeros(n,1);
    if i == m+1 %i = m+1 (y-direction)
        for j=1:n 
            b(j) = (2*u2(i-1,j) -2*u2(i,j))/hy1 + ...
                    f(t2,x(i),y(j)) + 2*u2(i,j)/dt + (2*ux(t1,x(m+1),y(j)))/hx; %b()
            if j == 1 %fixed i = m+1, j = 1
               b(j) = b(j) - (2*uy(t1,x(i),y(1)))/hy; %b()- 2*hy*mu*u_y(t,x,0)
               A(j,j+1) = -2/hy1; %-2*mu
            else
               if j==n %fixed i = m+1, j = n
                  b(j) = b(j) + uexact(t1,x(i),y(j+1))/hy1; %b()+ mu*u(t,x,1)
                  A(j,j-1) =  -1/hy1; %-mu
               else %fixed i = m+1, 1< j < n
                  A(j,j+1) = -1/hy1; %-mu
                  A(j,j-1) = -1/hy1; %-mu
               end
            end
            
         A(j,j) = 2/dt + 2/hy1;  %1+2*mu
         end
   
         ut = A\b;  % i = m+1
         for j=1:n
            u1(i,j) = ut(j);
         end
    else %1 < i < m+1  (y-direction) i = 2,...,m
         for j=1:n 
            b(j) = (u2(i-1,j) -2*u2(i,j)+ u2(i+1,j))/hy1 + ...
                    f(t2,x(i),y(j)) + 2*u2(i,j)/dt; %b()
            if j == 1 %fixed 1 < i < m+1, j = 1
               b(j) = b(j) - (2*uy(t1,x(i),y(1)))/hy; %b()- 2*hy*mu*u_y(t,x,0)
               A(j,j+1) = -2/hy1; %-2*mu
            else
               if j==n %fixed 1 < i < m+1, j = n
                  b(j) = b(j) + uexact(t1,x(i),y(j+1))/hy1; %b()+ mu*u(t,x,1)
                  A(j,j-1) =  -1/hy1; %-mu
               else %fixed 1 < i < m+1, 1< j < n
                  A(j,j+1) = -1/hy1; %-mu
                  A(j,j-1) = -1/hy1; %-mu
               end
            end
    
         A(j,j) = 2/dt + 2/hy1;  %1+2*mu
         end
    
         ut = A\b;  % i = 2,...,m
         for j=1:n
            u1(i,j) = ut(j);
         end
    end
 end     % Finish y-sweep.

当i = 1时，Dirichlet boundary condition，不用管。

当 1<i<m+1,当j=1时，要用ghost point，具体看note。得到如下矩阵，

Debug笔记

%-- Initial condition:
t = 0; 
for i = 1:m+1
    for j = 1:n+1
        u1(i,j) = uexact(t,x(i),y(j));
    end
end

这次第一次犯错是因为有一些情况只做了一次ghost point method，逻辑没对。后来把逻辑改对了，一直没得到正确的图是因为被initial condition误导了。当时知道左边界和上边界是Dirichlet，但是initial condition出来的矩阵是第一行和第一列全是0，其他都是有数据。这就表示左边界和下边界是Dirichlet，一直没相通为什么明明左边界和上边界才应该全是0才对（即第一行和最后一列全是0），当时最后一列不是0，但是也是约等于0。一直被这个坐标搞的很迷，最后决定不管了，自己按照自己的坐标一个个排好。然后就得到正确的图了。。。只能说别想太多，要相信自己。。。。

t = t + dt;