# planets / rk4.ml

 ``` 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 ```