home *** CD-ROM | disk | FTP | other *** search
/ CD Shareware Magazine 1996 December / CD_shareware_12-96.iso / DOS / Programa / NUMANA01.ZIP / DEF / INTEGRAT.DEF < prev    next >
Encoding:
Modula Definition  |  1996-08-16  |  6.2 KB  |  105 lines
  1. DEFINITION MODULE Integration;
  2.  
  3.         (********************************************************)
  4.         (*                                                      *)
  5.         (*              Solving differential equations          *)
  6.         (*                                                      *)
  7.         (*  Programmer:         P. Moylan                       *)
  8.         (*  Last edited:        16 August 1996                  *)
  9.         (*  Status:             Working                         *)
  10.         (*                                                      *)
  11.         (********************************************************)
  12.  
  13. FROM SYSTEM IMPORT (*type*) ADDRESS;
  14.  
  15. TYPE
  16.     (* User-supplied procedure to compute the right-hand side of the    *)
  17.     (* differential equation being solved.  The equation being solved   *)
  18.     (* is dy/dt = f(t,y), where t is a scalar and y is a vector.        *)
  19.     (* The procedure to compute f(t,y) has two input parameters and     *)
  20.     (* one result parameter.                                            *)
  21.  
  22.     RHSproc = PROCEDURE (LONGREAL, ARRAY OF LONGREAL, VAR (*OUT*) ARRAY OF LONGREAL);
  23.  
  24.     (* You have a choice of solution methods.  RK4 is the standard and  *)
  25.     (* popular fourth-order Runge-Kutta method with fixed step size.    *)
  26.  
  27.     (* RK5 is a fifth-order variant of this which adaptively updates    *)
  28.     (* the step size.  Unless you insist on having a fixed step size,   *)
  29.     (* RK5 will in general be more efficient and accurate; it does more *)
  30.     (* computation per step, but (except for particularly non-smooth    *)
  31.     (* problems) it compensates for this by taking bigger steps.        *)
  32.  
  33.     (* MM is the modified midpoint method - generally not as good as    *)
  34.     (* Runge-Kutta, but I threw it in for completeness.                 *)
  35.  
  36.     (* BS is the Bulirsch-Stoer method, which tends to be good when     *)
  37.     (* you want large step sizes.                                       *)
  38.  
  39.     SolutionMethod = (RK4, RK5, MM, BS);
  40.  
  41.     (* User-supplied procedure which is called each time a new solution *)
  42.     (* point has been calculated.  The second and third parameters are  *)
  43.     (* the current values of t and y.  The first parameter can be used  *)
  44.     (* for anything you like - title information, which screen window   *)
  45.     (* to use, etc.                                                     *)
  46.  
  47.     PlotProc = PROCEDURE (ADDRESS, LONGREAL, ARRAY OF LONGREAL);
  48.  
  49. (************************************************************************)
  50.  
  51. PROCEDURE Solve (N: CARDINAL;  RHS: RHSproc;  method: SolutionMethod;
  52.                         Plot: PlotProc;  Extra: ADDRESS;
  53.                         t0, tf, minstep, maxstep: LONGREAL;
  54.                         y0: ARRAY OF LONGREAL;  eps: LONGREAL);
  55.  
  56.     (* Solves an Nth degree differential equation dy/dt = f(t,y) from   *)
  57.     (* t=t0 to t=tf, calling procedure Plot at each step.  The          *)
  58.     (* parameters are:                                                  *)
  59.     (*          N       the number of elements in y                     *)
  60.     (*          RHS     procedure to compute f(t,y)                     *)
  61.     (*          method  see above for the definition of SolutionMethod. *)
  62.     (*          Plot    procedure which is called at each step, to      *)
  63.     (*                  tell the client the current values of t and y.  *)
  64.     (*          Extra   used as the first parameter to Plot.  You can   *)
  65.     (*                  set this to NIL, or you can use it to pass      *)
  66.     (*                  extra information that Plot needs to know,      *)
  67.     (*                  e.g. which of several DEs is being solved, or   *)
  68.     (*                  which screen window to use.                     *)
  69.     (*          t0      initial value of t                              *)
  70.     (*          tf      final value of t - we stop when a step takes us *)
  71.     (*                  up to or beyond tf.                             *)
  72.     (*          minstep minimum step size - needs to be moderately      *)
  73.     (*                  small, but not so small that you get stuck for  *)
  74.     (*                  a long time making negligible progress.  The    *)
  75.     (*                  "best" value depends on the smoothness of the   *)
  76.     (*                  solution, and may have to be established by     *)
  77.     (*                  trial and error.                                *)
  78.     (*          maxstep maximum step size - make this fairly large, so  *)
  79.     (*                  that the solver can cruise quickly along the    *)
  80.     (*                  long straight sections (if any) while still     *)
  81.     (*                  slowing down for the hairpin bends.  You might  *)
  82.     (*                  however have to limit maxstep for the sake of   *)
  83.     (*                  getting a smooth-looking graph plot.            *)
  84.     (*          y0      initial value of y.                             *)
  85.     (*          eps     specification of required relative error.       *)
  86.     (*                  Make this small, subject to the following       *)
  87.     (*                  proviso: if the actual step falls as low as     *)
  88.     (*                  minstep, then it's likely that you're no longer *)
  89.     (*                  getting the specified accuracy.  In this case   *)
  90.     (*                  you should decrease minstep, or increase eps,   *)
  91.     (*                  or simply accept that the error specification   *)
  92.     (*                  is no longer being met.                         *)
  93.     (*                                                                  *)
  94.     (* When using the fixed-step solution method RK4, maxstep is used   *)
  95.     (* as the step size, and minstep and eps are ignored.               *)
  96.     (*                                                                  *)
  97.     (* Method MM is also a fixed-step approach, but each step is        *)
  98.     (* internally subdivided into a sequence of smaller steps.  In      *)
  99.     (* this case maxstep is the step size, and maxstep/minstep is used  *)
  100.     (* as the number of internal divisions.  In this case too parameter *)
  101.     (* eps is ignored.                                                  *)
  102.  
  103. END Integration.
  104.  
  105.