heat2.m

使用ADI方法解偏微分方程（二微加热器方程）

function [x,y,U] = heat2 (T,a,b,m,n,p,beta,f,g)
%----------------------------------------------------------------------
% Usage:       [x,y,U] = heat2 (T,a,b,m,n,p,beta,f,g)
%          
% Description: Use the alternating direction implicit (ADI) method to
%              solve the following partial differential equation called
%              the two-dimensional heat equation (also the diffusion
%              equation):
%
%                 (d/dt)u(t,x,y) = beta*[(d/dx)(d/dx)u(t,x,y)
%                                  + (d/dy)(d/dy)u(t,x,y)]
%
%              The heat equation is to be solved over the region:
%
%                 0 < t <= T,
%                 0 < x < a,
%                 0 < y < b
%         
%              subject to the initial and boundary conditions:
%
%                 u(0,x,y) = f(x,y),
%                 u(t,x,y) = g(t,x,y)
%
% Inputs:      T    = final solution time (T >= 0)
%              a    = upper boundary of x (a > 0)
%              b    = upper boundary of y (b > 0)
%              m    = number of steps in t (m >= 1)
%              n    = number of steps in x (n >= 1)
%              p    = number of steps in y (p >= 1)
%              beta = heat equation coefficient (beta > 0)
%              f    = string containing name of user-supplied 
%                     initial condition function
%              g    = string containing name of user-supplied 
%                     boundary condition function
%
%                     The functions f and g are of the form:
%
%                        function z = f(x,y)
%                        function z = g(t,x,y)
%
%                     When f is called it must return the value u(0,x,y).
%                     When g is called with (x,y) on the boundary, it
%                     must return the boundary value u(t,x,y). The functions
%                     f and g are optional in the sense that either can be
%                     replaced by the empty string, ''. When this is done,
%                     the corresponding function is assumed to be zero 
%                     which means that a user-supplied function is NOT 
%                     required in this case.
%
% Outputs:     x = n by 1 vector of x grid points
%              y = p by 1 vector of y grid points     
%              U = n by p matrix containing final solution estimates
%                  U(i,j) = u(T,x(i),y(j)).
%
% Notes:       1. The function heat2 uses the functions trifac and trisub 
%                 to solve n+p tridiagonal linear algebraic systems at each
%                 time step.
%
%              2. The ADI method is a stable method. 
%----------------------------------------------------------------------
   
% Check arguments

   T = args (T,eps,T,1,'heat2');
   a = args (a,eps,a,2,'heat2');
   b = args (b,eps,b,3,'heat2');
   m = args (m,1,m,4,'heat2');
   n = args (n,1,n,5,'heat2');
   p = args (p,1,p,6,'heat2');
   beta = args (beta,eps,beta,7,'heat2');
   funf = getfun (f);
   fung = getfun (g);

% Initialize

   A1 = zeros (3,p);
   A2 = zeros (3,n);
   Q1 = zeros (3,p);
   Q2 = zeros (3,n);
   W = zeros (n,p);
   c1 = zeros (p,1);
   c2 = zeros (n,1);
   v2 = zeros (p,1);   
   x = zeros (n,1);
   y = zeros (p,1);
   U = zeros (n,p);
   dt = T/m;
   dx = a/(n+1);
   dy = b/(p+1);
   gx = beta*(dt/2)/(dx*dx);   
   gy = beta*(dt/2)/(dy*dy);
   for i = 1 : n 
      x(i) = i*dx;
   end
   for j = 1 : p
      y(j) = j*dy;   
   end
   
%* Compute coefficient matrices 

      g0 = 1 + 2*gy;
      for i = 1 : p
         A1(1,i) = -gy;
         A1(2,i) =  g0;
         A1(3,i) = -gy; 
      end
      [Q1,delta] = trifac (A1);
      if abs(delta) < eps/100
         disp ('Tridiagonal coefficient matrix in heats is singular.') 
         return
      end
      g0 = 1 + 2*gx;
      for i = 1 : n
         A2(1,i) = -gx;
         A2(2,i) =  g0;
         A2(3,i) = -gx;
      end
      [Q2,delta] = trifac (A2);
      if abs(delta) < eps/100
         disp ('Tridiagonal coefficient matrix in heats is singular.') 
         return
      end

% Initialize solution 

   if funf
      for i = 1 : n
         for j = 1 : p
            U(i,j) = feval(f,x(i),y(j));
         end
      end 
   end

% Solve heat equation  
         
   hwbar = waitbar (0,'Solving Heat Equation: heat2');
   for k = 0 : m-1

      waitbar (k/(m-1))
      t0 = k*dt;
      t1 = (k+1)*dt;
      tm = (t0 + t1)/2;
   
% First pass: i = 1 

      g0 = 1 - 2*gx;
      for j = 1 : p
         c1(j) = g0*U(1,j) + gx*U(2,j);
         if fung 
            c1(j) = c1(j) + gx*feval(g,t0,0,y(j));
         end
      end
      if fung      
         c1(1) = c1(1) + gy*feval(g,tm,x(1),0);
         c1(p) = c1(p) + gy*feval(g,tm,x(1),b);
      end
      [v2,delta] = trisub(Q1,c1);
      W(1,1:p) = v2';   

% First pass: 1 < i < n 

      for i = 2 : n-1
         for j = 1 : p
            c1(j) = gx*U(i-1,j) + g0*U(i,j) + gx*U(i+1,j);
         end 
         if fung   
            c1(1) = c1(1) + gy*feval(g,tm,x(i),0);
            c1(p) = c1(p) + gy*feval(g,tm,x(i),b);
         end
         [v2,delta] = trisub(Q1,c1);
         W(i,1:p) = v2'; 
      end
      
% First pass: i = n 

      for j = 1 : p
         c1(j) = gx*U(n-1,j) + g0*U(n,j);
         if fung
            c1(j) = c1(j) + gx*feval(g,t0,a,y(j));
         end
      end
      if fung      
         c1(1) = c1(1) + gy*feval(g,tm,x(n),0);
         c1(p) = c1(p) + gy*feval(g,tm,x(n),b);
      end
      [v2,delta] = trisub(Q1,c1);
      W(n,1:p) = v2';
   
% Second pass: j = 1 

      g0 = 1 - 2*gy;
      for i = 1 : n
         c2(i) = g0*W(i,1) + gx*W(i,2);
         if fung
            c2(i) = c2(i) + gy*feval(g,t0,x(i),0); 
         end
      end   
      if fung  
         c2(1) = c2(1) + gx*feval(g,t1,0,y(1));
         c2(n) = c2(n) + gx*feval(g,t1,a,y(1)); 
      end
      [U(1:n,1),delta] = trisub (Q2,c2);

% Second pass: 1 < j < p 

      for j = 2 : p-1
         for i = 1 : n
            c2(i) = gy*W(i,j-1) + g0*W(i,j) + gy*W(i,j+1);
         end 
         if fung 
            c2(1) = c2(1) + gx*feval(g,t1,0,y(j));
            c2(n) = c2(n) + gx*feval(g,t1,a,y(j));
         end
         [U(1:n,j),delta] = trisub (Q2,c2);
      end

% Second pass: j = p 

      for i = 1 : n
         c2(i) = gy*W(i,p-1) + g0*W(i,p);
         if fung 
            c2(i) = c2(i) + gx*feval(g,tm,x(i),b);
         end
      end
      if fung  	
         c2(1) = c2(1) + gx*feval(g,t1,0,y(p));
         c2(n) = c2(n) + gx*feval(g,t1,a,y(p)); 
      end 
      [U(1:n,p)delta] = trisub (Q2,c2);
   end
   fprintf ('\n')
   close (hwbar);
%----------------------------------------------------------------------