home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / ODE45.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.1 KB  |  110 lines
  1. function [tout, yout] = ode45(ypfun, t0, tfinal, y0, tol, trace)
  2. %ODE45    Solve differential equations, higher order method.
  3. %    ODE45 integrates a system of ordinary differential equations using
  4. %    4th and 5th order Runge-Kutta formulas.
  5. %    [T,Y] = ODE45('yprime', T0, Tfinal, Y0) integrates the system of
  6. %    ordinary differential equations described by the M-file YPRIME.M,
  7. %    over the interval T0 to Tfinal, with initial conditions Y0.
  8. %    [T, Y] = ODE45(F, T0, Tfinal, Y0, TOL, 1) uses tolerance TOL
  9. %    and displays status while the integration proceeds.
  10. %
  11. %    INPUT:
  12. %    F     - String containing name of user-supplied problem description.
  13. %            Call: yprime = fun(t,y) where F = 'fun'.
  14. %            t      - Time (scalar).
  15. %            y      - Solution column-vector.
  16. %            yprime - Returned derivative column-vector; yprime(i) = dy(i)/dt.
  17. %    t0    - Initial value of t.
  18. %    tfinal- Final value of t.
  19. %    y0    - Initial value column-vector.
  20. %    tol   - The desired accuracy. (Default: tol = 1.e-6).
  21. %    trace - If nonzero, each step is printed. (Default: trace = 0).
  22. %
  23. %    OUTPUT:
  24. %    T  - Returned integration time points (column-vector).
  25. %    Y  - Returned solution, one solution column-vector per tout-value.
  26. %
  27. %    The result can be displayed by: plot(tout, yout).
  28. %
  29. %    See also ODE23, ODEDEMO.
  30.  
  31. %    C.B. Moler, 3-25-87, 8-26-91, 9-08-92.
  32. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  33.  
  34. % The Fehlberg coefficients:
  35. alpha = [1/4  3/8  12/13  1  1/2]';
  36. beta  = [ [    1      0      0     0      0    0]/4
  37.           [    3      9      0     0      0    0]/32
  38.           [ 1932  -7200   7296     0      0    0]/2197
  39.           [ 8341 -32832  29440  -845      0    0]/4104
  40.           [-6080  41040 -28352  9295  -5643    0]/20520 ]';
  41. gamma = [ [902880  0  3953664  3855735  -1371249  277020]/7618050
  42.           [ -2090  0    22528    21970    -15048  -27360]/752400 ]';
  43. pow = 1/5;
  44. if nargin < 5, tol = 1.e-6; end
  45. if nargin < 6, trace = 0; end
  46.  
  47. % Initialization
  48. t = t0;
  49. hmax = (tfinal - t)/16;
  50. h = hmax/8;
  51. y = y0(:);
  52. f = zeros(length(y),6);
  53. chunk = 128;
  54. tout = zeros(chunk,1);
  55. yout = zeros(chunk,length(y));
  56. k = 1;
  57. tout(k) = t;
  58. yout(k,:) = y.';
  59.  
  60. if trace
  61.    clc, t, h, y
  62. end
  63.  
  64. % The main loop
  65.  
  66. while (t < tfinal) & (t + h > t)
  67.    if t + h > tfinal, h = tfinal - t; end
  68.  
  69.    % Compute the slopes
  70.    temp = feval(ypfun,t,y);
  71.    f(:,1) = temp(:);
  72.    for j = 1:5
  73.       temp = feval(ypfun, t+alpha(j)*h, y+h*f*beta(:,j));
  74.       f(:,j+1) = temp(:);
  75.    end
  76.  
  77.    % Estimate the error and the acceptable error
  78.    delta = norm(h*f*gamma(:,2),'inf');
  79.    tau = tol*max(norm(y,'inf'),1.0);
  80.  
  81.    % Update the solution only if the error is acceptable
  82.    if delta <= tau
  83.       t = t + h;
  84.       y = y + h*f*gamma(:,1);
  85.       k = k+1;
  86.       if k > length(tout)
  87.          tout = [tout; zeros(chunk,1)];
  88.          yout = [yout; zeros(chunk,length(y))];
  89.       end
  90.       tout(k) = t;
  91.       yout(k,:) = y.';
  92.    end
  93.    if trace
  94.       home, t, h, y
  95.    end
  96.  
  97.    % Update the step size
  98.    if delta ~= 0.0
  99.       h = min(hmax, 0.8*h*(tau/delta)^pow);
  100.    end
  101. end
  102.  
  103. if (t < tfinal)
  104.    disp('Singularity likely.')
  105.    t
  106. end
  107.  
  108. tout = tout(1:k);
  109. yout = yout(1:k,:);
  110.