home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / ODE23.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.6 KB  |  97 lines
  1. function [tout, yout] = ode23(ypfun, t0, tfinal, y0, tol, trace)
  2. %ODE23    Solve differential equations, low order method.
  3. %    ODE23 integrates a system of ordinary differential equations using
  4. %    2nd and 3rd order Runge-Kutta formulas.
  5. %    [T,Y] = ODE23('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] = ODE23(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-3).
  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 ODE45, 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. % Initialization
  35. pow = 1/3;
  36. if nargin < 5, tol = 1.e-3; end
  37. if nargin < 6, trace = 0; end
  38.  
  39. t = t0;
  40. hmax = (tfinal - t)/16;
  41. h = hmax/8;
  42. y = y0(:);
  43. chunk = 128;
  44. tout = zeros(chunk,1);
  45. yout = zeros(chunk,length(y));
  46. k = 1;
  47. tout(k) = t;
  48. yout(k,:) = y.';
  49.  
  50. if trace
  51.    clc, t, h, y
  52. end
  53.  
  54. % The main loop
  55.  
  56. while (t < tfinal) & (t + h > t)
  57.    if t + h > tfinal, h = tfinal - t; end
  58.  
  59.    % Compute the slopes
  60.    s1 = feval(ypfun, t, y); s1 = s1(:);
  61.    s2 = feval(ypfun, t+h, y+h*s1); s2 = s2(:);
  62.    s3 = feval(ypfun, t+h/2, y+h*(s1+s2)/4); s3 = s3(:);
  63.  
  64.    % Estimate the error and the acceptable error
  65.    delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');
  66.    tau = tol*max(norm(y,'inf'),1.0);
  67.  
  68.    % Update the solution only if the error is acceptable
  69.    if delta <= tau
  70.       t = t + h;
  71.       y = y + h*(s1 + 4*s3 + s2)/6;
  72.       k = k+1;
  73.       if k > length(tout)
  74.          tout = [tout; zeros(chunk,1)];
  75.          yout = [yout; zeros(chunk,length(y))];
  76.       end
  77.       tout(k) = t;
  78.       yout(k,:) = y.';
  79.    end
  80.    if trace
  81.       home, t, h, y
  82.    end
  83.  
  84.    % Update the step size
  85.    if delta ~= 0.0
  86.       h = min(hmax, 0.9*h*(tau/delta)^pow);
  87.    end
  88. end
  89.  
  90. if (t < tfinal)
  91.    disp('Singularity likely.')
  92.    t
  93. end
  94.  
  95. tout = tout(1:k);
  96. yout = yout(1:k,:);
  97.