suanli5_3_1.m

《偏微分方程数值解法》孙志忠版中算例的一些MATLAB程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                      %
%    Example of ADI Method for 2D heat equation                        %
%                                                                      %
%          u_t = u_{xx} + u_{yy} + f(x,t)                              %
%          a<x<b;c<y<d                                                 %
%    Test problme:                                                     %
%      Exact solution: u(t,x,y) = exp(-t) sin(pi*x) sin(pi*y)          %
%      Source term:  f(t,x,y) = exp(-t) sin(pi*x) sin(pi*y) (2pi^2-1)  %
%                                                                      %
%    Files needed for the test:                                        %
%                                                                      %
%     f.m:        The file defines the f(t,x,y)                        %
%     uexact.m:    The exact solution.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   clear; close all;

   a = 0; b=1;  c=0; d=1; n = 20;  tfinal = 1;

   m = n;
   h = (b-a)/n;           
   h1 = h*h;  dt=h1;
   r=dt/h1;
   x=a:h:b; y=c:h:d;

%-- Initial condition:

   t = 0;
   for i=1:m+1,
      for j=1:m+1,
         u1(i,j) = uexact(t,x(i),y(j));
      end
   end

%---------- Big loop for time t --------------------------------------

k_t = ceil(tfinal/dt);

for k=1:k_t

t1 = t + dt; t2 = t + dt/2;

%--- sweep in x-direction --------------------------------------

for j=2:m,                              % Boundary condition.
  u2(1,j) = (1/12-r/2)*uexact(t1,x(1),y(j-1)) + (5/6+r)*uexact(t1,x(1),y(j)) + (1/12-r/2) * uexact(t1,x(1),y(j+1));
  u2(m+1,j) = (1/12-r/2)*uexact(t1,x(m+1),y(j-1)) + (5/6+r)*uexact(t1,x(m+1),y(j)) + (1/12-r/2) * uexact(t1,x(m+1),y(j+1));
end

for j = 2:n,                             % Look for fixed y(j) 

b=zeros(m-1,1);
   for i=2:m,
      b(i-1) = (1/144+r/12+r^2/4) * ( u1(i-1,j-1) + u1(i-1,j+1) + u1(i+1,j-1) + u1(i+1,j+1) ) ...
	  + ( 10/144+r/3-r^2/2 )*( u1(i-1,j) + u1(i+1,j) + u1(i,j-1) + u1(i,j+1) )...
      + ( 100/144 - 40 * r/24 + r^2 )*u1(i,j)...
      + dt * (   f(t2,x(i-1),y(j-1)) + f(t2,x(i-1),y(j+1)) + f(t2,x(i+1),y(j-1)) + f(t2,x(i+1),y(j+1))...
      + 10 * (  f(t2,x(i-1),y(j)) + f(t2,x(i+1),y(j)) + f(t2,x(i),y(j-1)) + f(t2,x(i),y(j+1)) )...
      + 100 * f(t2,x(i),y(j))   ) / 144;
   end  

   b(1) = b(1) + ( r/2 -1/12 ) * u2(1,j);
   b(m-1) = b(m-1) + ( r/2 -1/12 ) * u2(m+1,j);
   A = diag( (5/6+r)*ones(m-1,1) ) + diag( (1/12-r/2) * ones(m-2,1),1 ) ...
   + diag( (1/12-r/2) * ones(m-2,1),-1);

     ut = A\b;                          % Solve the diagonal matrix.
     for i=1:m-1,
	     u2(i+1,j) = ut(i);
     end

 end                                    % Finish x-sweep.

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

for i=1:m+1,                                % Boundary condition
  u1(i,1) = uexact(t1,x(i),y(1));
  u1(i,n+1) = uexact(t1,x(i),y(m+1));
  u1(1,i) = uexact(t1,x(1),y(i));
  u1(m+1,i) = uexact(t1,x(m+1),y(i));
end

for i = 2:m,

   for j=2:m,
     b(j-1) = u2(i,j);
   end 

   b(1) = b(1) + ( r/2 -1/12 ) * u1(i,1);
   b(m-1) = b(m-1) + ( r/2 -1/12 ) * u1(i,m+1);

     ut = A\b;
     for j=1:n-1,
        u1(i,j+1) = ut(j);
     end

 end                             % Finish y-sweep.

 t = t + dt;

%--- finish ADI method at this time level, go to the next time level.
      
end       %-- Finished with the loop in time

%----------- Data analysis ----------------------------------

  for i=1:m+1,
    for j=1:n+1,
       ue(i,j) = uexact(tfinal,x(i),y(j));
    end
  end

[xx,yy]=meshgrid(x,y);
xx=xx';yy=yy';
mesh(xx,yy,u1)
title('数值解图')
pause
surf(xx,yy,ue)
title('精确解图')
pause
surf(xx,yy,u1-ue)
title('误差图') 

e = max(max(abs(u1-ue)))        % The infinity error.