suanli5_1_1dy.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:                                        %
%                                                                      %
%     adi.m:      This file, the main calling code.                    %
%     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 = 40;  tfinal = 1;

   m = n;
   h = (b-a)/n;        dt=h;   
   h1 = h*h;
   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 = fix(tfinal/dt);

for k=1:k_t

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

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

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

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

b=zeros(m-1,1);
   for i=2:m,
      b(i-1) = r^2/4*( u1(i-1,j-1) + u1(i-1,j+1) + u1(i+1,j-1) + u1(i+1,j+1) ) ...
	  + r/2*( 1-r )*( u1(i-1,j) + u1(i+1,j) + u1(i,j-1) + u1(i,j+1) )...
      + ( 1-2*r+r^2 )*u1(i,j) - 3/2*dt*exp( 1/2*( (i-1)*h + (j-1)*h ) -t2 );
   end  

   b(1) = b(1) + r/2 * u2(1,j);
   b(m-1) = b(m-1) + r/2 * u2(m+1,j);
   A = diag( (1+r)*ones(m-1,1) ) + diag( -r/2*ones(m-2,1),1 ) ...
   + diag(-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*u1(i,1);
   b(m-1) = b(m-1) + r/2*u1(i,m+1);
   A = diag( (1+r)*ones(m-1,1) ) + diag( -r/2*ones(m-2,1),1 ) ...
   + diag(-r/2*ones(m-2,1),-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.