File: rk4.ml

package info (click to toggle)
planets 0.1.13-19
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 512 kB
  • sloc: ml: 4,541; makefile: 207; ansic: 38
file content (64 lines) | stat: -rw-r--r-- 1,680 bytes parent folder | download | duplicates (9)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
let dydx = ref (Array.make 0 0.)
let dym = ref (Array.make 0 0.)
let dyt = ref (Array.make 0 0.)
let yt = ref (Array.make 0 0.)

let make_static_arrays dims = 
  begin
    dydx := Array.make dims 0. ;
    dym := Array.make dims 0. ;
    dyt := Array.make dims 0. ;
    yt := Array.make dims 0. ;
  end

let dims = ref 0  

(** Given values for the variables y[1..n] and their derivatives dydx[1..n]
  known at x, use the fourth-order Runge-Kutta method to advance the solution
  over an interval h and return the incremented variables as yout[1..n],
  which need not be a distinct array from y. The user supplies the routine
  derivs(x,y,dydx), which returns derivatives dydx at x. 
*)
let step ~y ~x ~h ~yout ~derivs = 

  (* We store the static arrays used by the routine globally, so they don't
     get reallocated each call.  The dimensions of the problem are not likely
     to change between calls...  We save 3 Array.make-s at the expense of a
     dereference and a compare. *)

  let n = Array.length y in
  let hh = h *. 0.5 
  and h6 = h /. 6.0 in
  let xh = x +. hh in
  let _ = (if n <> !dims then (dims:=n;make_static_arrays n)) in
  let dydx = !dydx
  and dym = !dym
  and dyt = !dyt
  and yt = !yt 
  in 
  derivs x y dydx ;

  for i = 0 to n-1 do 
    yt.(i) <- y.(i) +. hh *. dydx.(i) 
  done ;
  
  derivs xh yt dyt ;
  
  for i = 0 to n-1 do 
    yt.(i) <- y.(i) +. hh *. dyt.(i) 
  done ;
  
  derivs xh yt dym ; 

  for i = 0 to n-1 do 
    yt.(i) <- y.(i) +. h *. dym.(i) ; 
    dym.(i) <- dym.(i) +. dyt.(i)
  done ;
  
  derivs (x +. h) yt dyt ; 
  
  for i = 0 to n-1 do 
    yout.(i) <- y.(i) +. h6 *. (dydx.(i) +. dyt.(i) +. 2.0 *. dym.(i))
  done