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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
let f mu t y f =
f.(0) <- y.(1) ;
f.(1) <- -. y.(0) -. mu *. y.(1) *. (y.(0) *. y.(0) -. 1.)
let jac mu t y dfdy dfdt =
dfdy.{0, 0} <- 0. ;
dfdy.{0, 1} <- 1. ;
dfdy.{1, 0} <- -. 2. *. mu *. y.(0) *. y.(1) -. 1. ;
dfdy.{1, 1} <- -. mu *. (y.(0) *. y.(0) -. 1.) ;
dfdt.(0) <- 0. ;
dfdt.(1) <- 0.
let integ mu =
let step = Gsl_odeiv.make_step Gsl_odeiv.RK8PD ~dim:2 in
let control = Gsl_odeiv.make_control_y_new ~eps_abs:1e-6 ~eps_rel:0. in
let evolve = Gsl_odeiv.make_evolve 2 in
let system = Gsl_odeiv.make_system (f mu) ~jac:(jac mu) 2 in
let (t, t1, h, y) = (0., 100., 1e-6, [| 1.; 0. |]) in
let rec loop data t h =
if t < t1 then begin
let (t, h) =
Gsl_odeiv.evolve_apply evolve control step system
~t ~t1 ~h ~y in
loop ((t, y.(0), y.(1)) :: data) t h
end
else List.rev data
in
loop [] t h
let integ2 mu =
let step = Gsl_odeiv.make_step Gsl_odeiv.RK8PD ~dim:2 in
let control = Gsl_odeiv.make_control_y_new ~eps_abs:1e-6 ~eps_rel:0. in
let evolve = Gsl_odeiv.make_evolve 2 in
let system = Gsl_odeiv.make_system (f mu) ~jac:(jac mu) 2 in
let t1 = 100. in
let y = [| 1.; 0. |] in
let state = ref (0., 1e-2) in
let rec loop ti y = function
| (t, h) when t < ti ->
let new_state =
Gsl_odeiv.evolve_apply evolve control step system
~t ~t1:ti ~h ~y in
loop ti y new_state
| state -> state
in
let data = ref [] in
for i=1 to 100 do
let ti = float i *. t1 /. 100. in
state := loop ti y !state ;
let (t, _) = !state in
data := (t, y.(0), y.(1)) :: !data
done ;
!data
let integ3 mu =
let step = Gsl_odeiv.make_step Gsl_odeiv.RK4 ~dim:2 in
let system = Gsl_odeiv.make_system (f mu) ~jac:(jac mu) 2 in
let t1 = 100. in
let t = ref 0. in
let h = 1e-2 in
let y = [| 1.; 0. |] in
let yerr = Array.make 2 0. in
let dydt_in = Array.make 2 0. in
let dydt_out = Array.make 2 0. in
let dfdy = Gsl_matrix.create 2 2 in
let data = ref [] in
jac mu t y dfdy dydt_in ;
while !t < t1 do
Gsl_odeiv.step_apply step ~t:!t ~h ~y ~yerr ~dydt_in ~dydt_out system ;
Array.blit dydt_out 0 dydt_in 0 2 ;
t := !t +. h ;
data := (!t, y.(0), y.(1)) :: !data ;
done ;
List.rev !data
let with_temp_file f =
let tmp = Filename.temp_file "gnuplot_" ".tmp" in
let res = try f tmp with exn -> Sys.remove tmp ; raise exn in
Sys.remove tmp ; res
let gnuplot script =
with_temp_file
(fun tmp ->
let script_c = open_out tmp in
List.iter
(fun s -> Printf.fprintf script_c "%s\n" s)
script ;
close_out script_c ;
match Unix.system
(Printf.sprintf "gnuplot %s" tmp) with
| Unix.WEXITED 0 -> ()
| _ -> prerr_endline "hmm, problems with gnuplot ?" )
let main () =
Gsl_error.init () ;
if Array.length Sys.argv < 2 then exit 1 ;
let data =
match int_of_string Sys.argv.(1) with
| 3 -> integ3 10.
| 2 -> integ2 10.
| _ -> integ 10.
in
let points = List.map
(fun (_, x, y) -> Printf.sprintf "%.5f %.5f" x y) data in
gnuplot (List.concat
[ [ "plot '-' with lines" ] ;
points ;
[ "e" ; "pause -1" ] ])
let _ =
main ()
|
