suanli5_2_1.m

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

%******%表示原来程序编写错误
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                      %
%    Example of ADI Method for 2D heat equation                        %
%                                                                      %
%          u_tt = u_{xx} + u_{yy} + f(x,t)                             %
%          a<x<b;c<y<d                                                 %
%    Test problme:                                                     %
%      精确解: u(t,x,y) =     exp( 1/2*(x+y) - t )                     %
%      右端函数:  f(t,x,y) =   1/2*exp( 1/2*(x+y) - t )                %
%      初始条件:  g1 =    exp( 1/2*(x+y) )                             %
%      初始导数条件: g2 = -exp( 1/2*(x+y) )                            %
%      需要函数                                                        %
%      g1.m g2.m    初始条件                                           %
%     f2.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 = 10;  tfinal = 1;

   m = n;
   h = (b-a)/n;        dt=h;   
   s=(dt/h)^2;
   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) = g1(x(i),y(j));
      end
   end
%--  Initial deriavtive condition: 
   for i=1:m+1,
   for j=1:m+1,
         %******% u2(i,j) = g1(x(i),y(j)) + dt * g2( x(i),y(j) ) + dt^2/4 * (  g1(x(i),y(j)) + f2 (t, x(i),y(j) )  );
         u2(i,j) = g1(x(i),y(j)) + dt * g2( x(i),y(j) ) + dt^2/2 * (  g1(x(i),y(j))/2 + f2 (t, x(i),y(j) )  );
   end
   end
%---------- Big loop for time t --------------------------------------

k_t = fix(tfinal/dt);

for k=2:k_t

t1 = t + dt; t2 = t1 + dt;

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

for i=1:m+1,                              % Boundary condition.
  u4(1,i) = ( uexact(t,x(1),y(i)) - 2*uexact(t1,x(1),y(i)) + uexact(t2,x(1),y(i)) )/dt^2  ;
  u4(m+1,i) = ( uexact(t,x(m+1),y(i)) - 2*uexact(t1,x(m+1),y(i)) + uexact(t2,x(m+1),y(i)) )/dt^2  ;
  u4(i,1) = ( uexact(t,x(i),y(1)) - 2*uexact(t1,x(i),y(1)) + uexact(t2,x(i),y(1)) )/dt^2  ;
  u4(i,m+1) = ( uexact(t,x(i),y(m+1)) - 2*uexact(t1,x(i),y(m+1)) + uexact(t2,x(i),y(m+1)) )/dt^2  ;
end

for i=2:m
    u3(1,i)=-s/2*u4(1,i-1) + (1+s)*u4(1,i) - s/2*u4(1,i+1); 
    u3(m+1,i)=-s/2*u4(m+1,i-1) + (1+s)*u4(m+1,i) - s/2*u4(m+1,i+1);
end

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

   b=zeros(m-1,1);
   for i=2:m,
      b(i-1) = ( u2(i-1,j) + u2(i+1,j) + u2(i,j-1) + u2(i,j+1) -4*u2(i,j) )/h^2 + f2 (t1, x(i),y(j) );
   end  

   b(1) = b(1) + s/2 * u3(1,j);
   b(m-1) = b(m-1) + s/2 * u3(m+1,j);
   A = diag( (1+s)*ones(m-1,1) ) + diag( -s/2*ones(m-2,1),1 ) ...
   + diag(-s/2*ones(m-2,1),-1);

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

 end                                    % Finish x-sweep.

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


for i = 2:m,

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

   b(1) = b(1) + s/2*u4(i,1);
   b(m-1) = b(m-1) + s/2*u4(i,m+1);

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

 end                             % Finish y-sweep.
 
 for i=1:m+1,                                % Boundary condition
  u5(i,1) = uexact(t2,x(i),y(1));
  u5(i,m+1) = uexact(t2,x(i),y(m+1));
  u5(1,i) = uexact(t2,x(1),y(i));
  u5(m+1,i) = uexact(t2,x(m+1),y(i));
 end
 
    for i=2:m,
      for j=2:m,
         u5(i,j) = dt^2 * u4(i,j) + 2*u2(i,j) - u1(i,j);
      end
   end

 u1=u2;
 u2=u5; 
 %******%t = t2 + dt;
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,u2)
title('数值解图')
pause
surf(xx,yy,ue)
title('精确解图')
pause
%******%surf(xx,yy,u2-ue)
surf(xx,yy,u2-ue)
title('误差图') 

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