@@ -12,9 +12,35 @@ function lotka(du, u, p, t)
1212 du[1 ] = p[1 ] * x - p[2 ] * x * y
1313 du[2 ] = - p[3 ] * y + p[4 ] * x * y
1414end
15- prob = ODEProblem {true, SciMLBase.FullSpecialize} (
15+ prob_lotka = ODEProblem {true, SciMLBase.FullSpecialize} (
1616 lotka, [1.0 , 1.0 ], (0.0 , 10.0 ), [1.5 , 1.0 , 3.0 , 1.0 ])
17- sol = solve (prob, Vern7 (), abstol = 1 / 10 ^ 14 , reltol = 1 / 10 ^ 14 )
17+ sol_lotka = solve (prob_lotka, Vern7 (), abstol = 1e-14 , reltol = 1e-14 )
18+
19+ function fitzhugh (du, u, p, t)
20+ v = u[1 ]
21+ w = u[2 ]
22+ a = p[1 ]
23+ b = p[2 ]
24+ τinv = p[3 ]
25+ l = p[4 ]
26+ du[1 ] = v - v^ 3 / 3 - w + l
27+ du[2 ] = τinv * (v + a - b * w)
28+ end
29+ prob_fitzhugh = ODEProblem {true, SciMLBase.FullSpecialize} (
30+ fitzhugh, [1.0 ; 1.0 ], (0.0 , 10.0 ), [0.7 , 0.8 , 1 / 12.5 , 0.5 ])
31+ sol_fitzhugh = solve (prob_fitzhugh, Vern7 (), abstol = 1e-14 , reltol = 1e-14 )
32+
33+ function rigid (du, u, p, t)
34+ I₁ = p[1 ]
35+ I₂ = p[2 ]
36+ I₃ = p[3 ]
37+ du[1 ] = I₁ * u[2 ] * u[3 ]
38+ du[2 ] = I₂ * u[1 ] * u[3 ]
39+ du[3 ] = I₃ * u[1 ] * u[2 ] + 0.25 * sin (t)^ 2
40+ end
41+ prob_rigid = ODEProblem {true, SciMLBase.FullSpecialize} (
42+ rigid, [1.0 ; 0.0 ; 0.9 ], (0.0 , 10.0 ), [- 2.0 , 1.25 , - 0.5 ])
43+ sol_rigid = solve (prob_rigid, Vern7 (), abstol = 1e-14 , reltol = 1e-14 );
1844
1945# Make work-precision plot
2046setups = [Dict (:alg => DP5 ())
@@ -26,6 +52,13 @@ for order in 6:2:12
2652 push! (names, " Taylor $(order) " )
2753 push! (setups, Dict (:alg => ExplicitTaylor (order = Val (order + 1 ))))
2854end
29- wp = WorkPrecisionSet ([prob], abstols, reltols, setups; names = names, appxsol = [sol],
30- save_everystep = false , numruns = 100 , maxiters = 10000 )
31- plot (wp)
55+
56+ function make_plot (prob, sol)
57+ wp = WorkPrecisionSet ([prob], abstols, reltols, setups; names = names, appxsol = [sol],
58+ save_everystep = false , numruns = 100 , maxiters = 10000 )
59+ p = plot (wp)
60+ return p
61+ end
62+ p_lotka = make_plot (prob_lotka, sol_lotka)
63+ p_fitzhugh = make_plot (prob_fitzhugh, sol_fitzhugh)
64+ p_rigid = make_plot (prob_rigid, sol_rigid)
0 commit comments