e872.m

matlab算法集 matlab算法集

%------------------------------------------------------------------
% Example 8.7.2: Stiff Differential Equations                       
%------------------------------------------------------------------
   clc
   clear
   q = 5000;
   n = 2;
   m = 11;
   alpha = 0;	
   beta	= 1.0;
   h = (beta-alpha)/(m-1);
   tol = 1.e-6;
   t = zeros(m,1);
   X0 = [1 -1; 1 -1; 1 -1]';
   E = zeros (m,3);
   P = zeros (m,3);
   f  = inline ('[-x(1) - x(2); -x(1) - 10001*x(2)]','t','x');
   fx = inline ('[-1 -1; -1 -10001]','t','x');
   ft = inline ('[0; 0]','t','x');

% Initialize 

   for k = 1 : m
      t(k) = alpha + (k-1)*h;
   end
      
% Solve stiff system 

   fprintf ('Example 8.7.2: Stiff Differential Equations\n');
   hwbar = waitbar (0,'Solving Stiff Differential Equations');
   for k = 1 : m-1
      waitbar(k/(m-1))
      [X0(:,1),E(k,1),P(k,1)] = stiff0 (X0(:,1),t(k),t(k+1),tol,q,f,fx,ft);
      [X0(:,2),E(k,2),P(k,2)] = stiff0 (X0(:,2),t(k),t(k+1),tol,q,f,'','');
      [X0(:,3),E(k,3),h,P(k,3)] = rkf0 (X0(:,3),t(k),t(k+1),h,tol,q,f);
   end
   close (hwbar)

% Compute average values       

   for j = 1 : 3
      for k = 1 : m-1
         E(m,j) = E(m,j) + E(k,j);
         P(m,j) = P(m,j) + P(k,j);
      end
      E(m,j) = E(m,j)/(m-1);
      P(m,j) = P(m,j)/(m-1);
   end    

% Display results 
   
   fprintf ('\n t     BS Extrap.    evaluations       RKF');
   fprintf ('        evaluations\n');
   fprintf ('----------------------------------------------');
   fprintf ('---------------\n');   
   for k = 1 : m-1
      fprintf ('%3.1f  %12.7e   %6i      %12.7e    %6i\n',...
               t(k+1),E(k,1),P(k,1),E(k,3),P(k,3)); 
   end
   fprintf ('----------------------------------------------');
   fprintf ('---------------\n');   
   fprintf ('ave. %12.7e    %6.1f     %12.7e     %6.1f\n',...
            E(m,1),P(m,1),E(m,3),P(m,3));
   fprintf ('----------------------------------------------');
   fprintf ('---------------\n');   
   wait
%------------------------------------------------------------------