home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / ODE23P.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.9 KB  |  97 lines
  1. function y = ode23p(FunFcn, t0, tfinal, y0, tol)
  2. %ODE23P Solve differential equations, low order method, displaying plot.
  3. %    Integrate a system of ordinary differential equations using
  4. %    2nd and 3rd order Runge-Kutta formulas and plot the time evolution
  5. %    of the first three components of the solution.
  6. %    [T,Y] = ODE23('yprime', T0, Tfinal, Y0) integrates the system
  7. %    of ordinary differential equations described by the M-file
  8. %    YPRIME.M over the interval T0 to Tf and using initial
  9. %    conditions Y0.
  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. %
  22. %    In addition, setup the graphics with:
  23. %       axis([y1min y1max y2min y2max y3min y3max]);
  24. %       hold
  25. %
  26. %    INTERMEDIATE RESULTS:
  27. %
  28. %    A 3-D plot of the orbit of y(1:3) as functions of t.
  29. %
  30. %    FINAL OUTPUT:
  31. %    y     - The solution at tfinal.
  32.  
  33. %    C.B. Moler, 3-25-87, 10-5-91.
  34. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  35.  
  36. % Initialization
  37. pow = 1/3;
  38. if nargin < 5, tol = 0.001; end
  39.  
  40. t = t0;
  41. hmax = (tfinal - t)/5;
  42. hmin = (tfinal - t)/20000;
  43. h = (tfinal - t)/100;
  44. y = y0(:);
  45. tau = tol * max(norm(y, 'inf'), 1);
  46.  
  47. % Save L steps and plot like a comet tail.
  48. L = 50;
  49. Y = y*ones(1,L);
  50.  
  51. cla;
  52. head = line('color','c','linestyle','o','erase','xor','xdata',y(1),'ydata',y(2), ...
  53.        'zdata',y(3));
  54. body = line('color','y','linestyle','-','erase','none','xdata',[],'ydata',[], ...
  55.        'zdata',[]);
  56. tail = line('color','m','linestyle','-','erase','none','xdata',[],'ydata',[], ...
  57.        'zdata',[]);
  58.  
  59. % The main loop
  60.    while (t < tfinal) & (h >= hmin)
  61.       if t + h > tfinal, h = tfinal - t; end
  62.  
  63.       % Compute the slopes
  64.       s1 = feval(FunFcn, t, y);
  65.       s2 = feval(FunFcn, t+h, y+h*s1);
  66.       s3 = feval(FunFcn, t+h/2, y+h*(s1+s2)/4);
  67.  
  68.       % Estimate the error and the acceptable error
  69.       delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');
  70.       tau = tol*max(norm(y,'inf'),1.0);
  71.  
  72.       % Update the solution only if the error is acceptable
  73.       ts = t;
  74.       ys = y;
  75.       if delta <= tau
  76.          t = t + h;
  77.          y = y + h*(s1 + 4*s3 + s2)/6;
  78.  
  79.          % Update the plot
  80.          Y = [y Y(:,1:L-1)];
  81.          set(head,'xdata',Y(1,1),'ydata',Y(2,1),'zdata',Y(3,1))
  82.          set(body,'xdata',Y(1,1:2),'ydata',Y(2,1:2),'zdata',Y(3,1:2))
  83.          set(tail,'xdata',Y(1,L-1:L),'ydata',Y(2,L-1:L),'zdata',Y(3,L-1:L))
  84.          drawnow;
  85.       end
  86.  
  87.       % Update the step size
  88.       if delta ~= 0.0
  89.          h = min(hmax, 0.9*h*(tau/delta)^pow);
  90.       end
  91.    end;
  92.  
  93.    if (t < tfinal)
  94.       disp('Singularity likely.')
  95.       t
  96.    end
  97.