home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / PROGRAMM.PAK / RUNGEKUT.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  1.3 KB  |  45 lines
  1.  
  2. (*********************************************************************
  3.  
  4.         Adapted from
  5.         Roman E. Maeder: Programming in Mathematica,
  6.         Second Edition, Addison-Wesley, 1991.
  7.  
  8.  *********************************************************************)
  9.  
  10.  
  11. BeginPackage["RungeKutta`"]
  12.  
  13. RungeKutta::usage = "RungeKutta[{e1,e2,..}, {y1,y2,..}, {a1,a2,..}, {t1, dt}]
  14.     numerically integrates the ei as functions of the yi with inital values ai.
  15.     The integration proceeds in steps of dt from 0 to t1.
  16.     RungeKutta[{e1,e2,..}, {y1,y2,..}, {a1,a2,..}, {t, t0, t1, dt}] integrates
  17.     a time-dependent system from t0 to t1."
  18.  
  19. Begin["`Private`"]
  20.  
  21. RKStep[f_, y_, y0_, dt_] :=
  22.     Module[{ k1, k2, k3, k4 },
  23.         k1 = dt N[ f /. Thread[y -> y0] ];
  24.         k2 = dt N[ f /. Thread[y -> y0 + k1/2] ];
  25.         k3 = dt N[ f /. Thread[y -> y0 + k2/2] ];
  26.         k4 = dt N[ f /. Thread[y -> y0 + k3] ];
  27.         y0 + (k1 + 2 k2 + 2 k3 + k4)/6
  28.     ]
  29.  
  30. RungeKutta[f_List, y_List, y0_List, {t1_, dt_}] :=
  31.     NestList[ RKStep[f, y, #, N[dt]]&, N[y0], Round[N[t1/dt]] ] /;
  32.         Length[f] == Length[y] == Length[y0]
  33.  
  34. RungeKutta[f_List, y_List, y0_List, {t_, t0_, t1_, dt_}] :=
  35.     Module[{res},
  36.         res = RungeKutta[ Append[f, 1], Append[y, t], Append[y0, t0], {t1 - t0, dt} ];
  37.         Drop[#, -1]& /@ res
  38.     ]  /;  Length[f] == Length[y] == Length[y0]
  39.  
  40. End[]
  41.  
  42. Protect[RungeKutta]
  43.  
  44. EndPackage[]
  45.