e831.m

matlab算法集 matlab算法集

%------------------------------------------------------------------
% Example 8.3.1: Runge-Kutta Methods 
%------------------------------------------------------------------
   clc
   clear
   n = 2;
   p = 6;
   r = 3;
   m0 = 5*(2^(p-1))+1;
   alpha = 0;
   beta = 1;
   z = zeros (n,1);
   h = zeros (p,1);
   Q = zeros (n,r);
   Y = zeros (m0,r+1);
   E = zeros (p,r);
   y = inline ('1 - (10/9)*exp(-t) + (1/9)*exp(-10*t)','t');
   f = inline ('[x(2); 10 - 10*x(1) - 11*x(2)]','t','x');
   fprintf ('\nExample 8.3.1: Runge-Kutta Methods\n');
   for j = 1 : p
      m = 5*(2^(j-1)) + 1;
      h(j) = (beta - alpha)/(m - 1);
      
% First-order RK method (Euler) 

      x = zeros(n,1);
      for k = 1 : m
          t = alpha + (k-1)*h(j);
         Q(:,1) = f(t,x);
         for i = 1 : n
            x(i) = x(i) + h(j)*Q(i,1);
         end
         Y(k,1) = x(1);   
         E(j,1) = max(E(j,1),abs(y(t)-x(1)));   
      end

% Second-order RK method 

      x = zeros(n,1);
      for k = 1 : m
         t = alpha + (k-1)*h(j);
         Q(:,1) = f(t,x);
         for i = 1 : n
            z(i) = x(i) + h(j)*Q(i,1);
         end
         Q(:,2) = f(t+h(j),z);
         for i = 1 : n
            x(i) = x(i) + (h(j)/2)*(Q(i,1) + Q(i,2));
         end
         Y(k,2) = x(1);   
         E(j,2) = max(E(j,2),abs(y(t)-x(1)));   
      end

% Third-order RK method 

      x = zeros(n,1);
      for k = 1 : m
         t = alpha + (k-1)*h(j);
         Q(:,1) = f(t,x);
         for i = 1 : n
            z(i) = x(i) + (h(j)/2)*Q(i,1);
         end
         Q(:,2) = f(t+h(j)/2,z);
         for i = 1 : n
            z(i) = x(i) - h(j)*Q(i,1) + 2*h(j)*Q(i,2);
         end
         Q(:,3) = f(t+h(j),z);  
         for i = 1 : n
            x(i) = x(i) + (h(j)/6)*(Q(i,1) + 4*Q(i,2) + Q(i,3));
         end
         Y(k,3) = x(1);
         E(j,3) = max(E(j,3),abs(y(t)-x(1)));   
         Y(k,4) = y(t);    
      end
   end

% Display results 

   show ('step sizes',h) 
   graphmat (Y,'Solutions','k','y')
   graphxy (h,E,'Errors','h','E')
%------------------------------------------------------------------