home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / diffeq_23.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  5.9 KB  |  214 lines
  1. ; $Id: diffeq_23.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ; Copyright (c) 1991-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. pro diffeq_23, funct, init, start, finish, times, yvalues, tol=tol, $
  7.                report = report, Params= params, Listname = listname, $
  8.                Depvar = depvar
  9. ;+
  10. ; NAME:
  11. ;       DIFFEQ_23
  12. ;
  13. ; PURPOSE:
  14. ;    Solve a system of first-order, ordinary differential equations:
  15. ;                    yi' = fi(t, y1(t), ... yn(y)), i = 1,..., n
  16. ;                    ai  = yi(start),               i = 1,..., n         
  17. ;    using the Runge-Kutta method of order 2 and 3. Step size is selected
  18. ;    automatically and hence is variable. 
  19. ;         
  20. ; CATEGORY:
  21. ;    Mathematical Functions, General
  22. ;
  23. ; CALLING SEQUENCE:
  24. ;       DIFFEQ_23, Funct, Init, Start, Finish, Times, Yvalues, $
  25. ;                  TOL = Tol, PARAMS = Params, REPORT = Report, $
  26. ;                  LISTNAME = Listname, DEPVAR = Depvar
  27. ;
  28. ; INPUTS:
  29. ;    Funct:  A character string containing the name of the user-supplied
  30. ;                  function implementing f = [f1, ...., fn]. This function
  31. ;                  should be written in IDL and have two arguments -- the scalar-
  32. ;        valued time argument t, and the vector argument 
  33. ;        [y1(t), ... ,yn(t)].  Additional constant parameters may be
  34. ;        supplied through the keyword PARAMS.
  35. ;
  36. ;      Init:    The vector [a1, ..., an] of initial values.
  37. ;
  38. ;      Start:    The initial value of t.
  39. ;
  40. ;      Finish:    The final value of t. 
  41. ;  
  42. ; KEYWORD PARAMETERS:
  43. ;      TOL:    The error tolerance. The default is .001.
  44. ;
  45. ;      PARAMS:    A keyword to be passed to the function f.  PARAMS can be used 
  46. ;        to specify constant-parameter values if f is a parametric 
  47. ;        family of functions.  See the example below.
  48. ;
  49. ;        If the IDL function to compute f does accept the keyword 
  50. ;        PARAMS, then PARAMS should not be set in the call to DIFFEQ_23.
  51. ;
  52. ;      REPORT:    If set, this flag signals that, at each step, the time
  53. ;        value, step size, and dependent variable values should
  54. ;        be written to the screen or to a file specified by keyword
  55. ;        LISTNAME.    
  56. ;
  57. ;    LISTNAME:    The name of the file to receive any output. The default is
  58. ;               to write to the screen. 
  59. ;
  60. ;      DEPVAR:     A string array of the names of the dependent variables to
  61. ;               be used in the output. Depvar(i) = name of variable i.
  62. ;
  63. ; OUTPUT PARAMETERS:
  64. ;          Times:    A vector of times at which f is computed.
  65. ;
  66. ;     Yvalues:    An array of y values. If ti = times(i),
  67. ;
  68. ;                 Yvalues(*, i) = f(ti,y1(ti),..., yn(ti))
  69. ;                               = [f1(ti,y1(ti),..., yn(ti)),  ..., 
  70. ;                                              fn(ti,y1(ti), ...,yn(ti))].
  71. ;
  72. ; COMMON BLOCKS:
  73. ;    None.
  74. ;
  75. ; SIDE EFFECTS:
  76. ;    None.
  77. ;
  78. ; RESTRICTIONS:
  79. ;    None.
  80. ;
  81. ; EXAMPLE:
  82. ;       Solve the set of equations:
  83. ;
  84. ;              y1' = -.1 * y1, 
  85. ;              y2' = .1*y1 - .05*y2,
  86. ;              y3' = .05*y2
  87. ; 
  88. ;              y1(0) = 1000, y2(0) = 0, y3(0) = 0
  89. ;
  90. ;       on the interval [0, 5].
  91. ;
  92. ;       First, we define the function RADIO as
  93. ;
  94. ;        FUNCTION RADIO, t, y, PARAMS = params
  95. ;        k = params(0)
  96. ;        kp = params(1)
  97. ;        RETURN, [-k*y(0), k*y(0) - kp*y(1), kp* y(1)]
  98. ;        END
  99. ;
  100. ;          Next, call DIFFEQ_23:
  101. ;
  102. ;                DIFFEQ_23, "radio", [1000, 0, 0], 0, 5., times, yvalues, $
  103. ;                          PARAMS = [.1, .05], /REPORT
  104. ;
  105. ;          The result can be plotted by entering:
  106. ;
  107. ;              PLOT, times, yvalues(0,*)
  108. ;              FOR i = 1,2 DO OPLOT, times, yvalues(i,*)
  109. ;
  110. ; MODIFICATION HISTORY:
  111. ;    CAB, Sept., 1991.
  112. ;-
  113.  
  114. On_Error,2
  115.  
  116. ; Check parameters
  117. if n_params(0) lt 4 THEN $
  118.    message, "Missing parameters"
  119.  
  120. if KEYWORD_SET(listname) THEN $
  121.    openw, unit, /get, listname $
  122. else unit = -1
  123.  
  124. if KEYWORD_SET(tol) eq 0 THEN tol =.001
  125. if KEYWORD_SET(Params) eq 0 THEN  $
  126.   Paramset = 0   $
  127. else Paramset = 1
  128.  
  129.  
  130. if KEYWORD_SET(report) eq 0 THEN report = 0 $
  131. ELSE  BEGIN
  132.    SN = size( depvar)
  133.    n = n_elements(init)
  134.    if (SN(1) eq 0) THEN BEGIN
  135.       I = indgen(n)
  136.       Names = ['Var' + StrTrim(I, 2)]
  137.    ENDIF ELSE  $
  138.       if SN(1) lt n THEN BEGIN
  139.       I = Indgen(N)
  140.       Names = [DepVar, 'Var' + StrTrim(I(SN(1) : N-1),2)]
  141.       ENDIF else Names = depvar
  142.   
  143.    printf, unit, format = '(A13, 2x, A13, 2x, $)', "Times", "Stepsize"
  144.    for i = 0, n-2 do printf,unit,format = '(A13,2x,$)',Names(i)
  145.    printf,unit, format = '(A13)',Names(n-1)
  146.   Printf,unit, " "
  147. ENDELSE
  148.  
  149. ; Initialize
  150. h = (finish - start)   ;h = stepsize
  151. minh  = h/20000        
  152. maxh  = h/5
  153. times = start
  154. h = h/100.
  155. t = start
  156. v = init
  157. yvalues = v
  158.  
  159. if report  ne 0 THEN BEGIN
  160.    printf,unit, format ='(G13.6, 2x, G13.6, 2x, $)', t, h
  161.    for i = 0, n-2 do printf,unit, format ='(G13.6, 2x,$)', v(i)
  162.    printf,unit,format = '(G13.6)', v(n-1)
  163. ENDIF
  164.  
  165.  
  166. errbound = tol * max([sqrt(total(v^2)), 1])
  167.  
  168. ; compute yvalues for variable step sizes
  169.  
  170.  
  171.    while t lt finish and h ge minh DO BEGIN
  172.  
  173.      if t+h gt finish THEN h = finish - t
  174.      
  175.      if Paramset eq 0 THEN BEGIN
  176.        k1 = CALL_FUNCTION(funct,t, v)
  177.        k2 = CALL_FUNCTION(funct, t+h, v + h * k1)
  178.        k3 = CALL_FUNCTION(funct, t + h/2, v + h*(k1 + k2)/4)
  179.     ENDIF ELSE BEGIN
  180.        k1 = CALL_FUNCTION(funct,t, v, Params = Params)
  181.        k2 = CALL_FUNCTION(funct, t+h, v + h * k1, Params = Params)
  182.        k3 = CALL_FUNCTION(funct, t + h/2, v + h*(k1 + k2)/4, Params = Params)
  183.     ENDELSE
  184.  
  185.     err = (h*(k1 - 2*k3 + k2)/3)^2
  186.     err = sqrt(total(err))
  187.     err_bound = tol * max([sqrt(total(v^2)), 1])
  188.  
  189.  if  err le err_bound THEN BEGIN
  190.     t = t + h
  191.     v = v + h*(k1 + 4*k3 + k2)/6
  192.     times   = [times,t]
  193.     yvalues = [[yvalues], [v]]
  194.     if report  ne 0 THEN BEGIN
  195.        printf,unit, format ='(G13.6, 2x, G13.6, 2x, $)', t, h
  196.       for i = 0, n-2 do printf,unit, format ='(G13.6, 2x,$)', v(i)
  197.       printf,unit,format = '(G13.6)', v(n-1)
  198.    ENDIF
  199. ENDIF
  200.  
  201. if err ne 0 THEN  $
  202.    h = min([maxh, .9*h*(err_bound / err)^(1/3.0)])
  203.  
  204.  endwhile
  205.  
  206.  
  207.  if ( t lt finish) THEN $
  208.      print, " Beware of singularity"
  209.  
  210. if unit ne -1 THEN Free_lun, unit
  211.  
  212. return
  213. end
  214.