Runge-Kutta methods for ODE integration in OCaml¶
I want to implement and illustrate the Runge-Kutta method (actually, different variants), in the OCaml programming language.
The Runge-Kutta methods are a family of numerical iterative algorithms to approximate solutions of Ordinary Differential Equations. I will simply implement them, for the mathematical descriptions, I let the interested reader refer to the Wikipedia page, or any good book or course on numerical integration of ODE.
I will start with the order 1 method, then the order 2 and the most famous order 4.
They will be compared on different ODE.
Preliminary¶
In [1]:
Sys.command "ocaml -version";;
The OCaml toplevel, version 4.04.2
Out[1]:
- : int = 0
In [147]:
#thread ;;
#require "jupyter.notebook" ;;
#require "jupyter.archimedes" ;;
I don't want to try to find a reference implementation of an Ordinary Differential Equations solver in OCaml. I'm sure there is, just don't care.
I will use as a first example the one included in the scipy (Python) documentation for this odeint function.
If $\omega(t) := \theta'(t)$, this gives $$ \begin{cases} \theta'(t) = \omega(t) \\ \omega'(t) = -b \omega(t) - c \sin(\theta(t)) \end{cases} $$
Vectorially, if $y(t) = [\theta(t), \omega(t)]$, then the equation is $y' = f(t, y)$ where $f(t, y) = [y_2(t), -b y_2(t) - c \sin(y_1(t))]$.
We assume the values of $b$ and $c$ to be known, and the starting point to be also fixed:
In [11]:
let b = 0.25 ;;
let c = 5.0 ;;
Out[11]:
val b : float = 0.25
Out[11]:
val c : float = 5.
In [14]:
let y0 = [| 3.14156 -. 0.1; 0.0 |];;
Out[14]:
val y0 : float array = [|3.04156; 0.|]
In [148]:
let f_pend (y : float array) (_ : float) : float array =
[|
y.(1);
(-.b) *. y.(1) +. (-.c) *. sin(y.(0))
|]
;;
Out[148]:
val f_pend : float array -> float -> float array = <fun>
Runge-Kutta method of order 1, or the Euler method¶
The approximation is computed using this update: $$y_{n+1} = y_n + (t_{n+1} - t_n) f(y_n, t_n).$$
The math behind this formula are the following: if $g$ is a solution to the ODE, and so far the approximation is correct, $y_n \simeq g(t_n)$, then a small step $h = t_{n+1} - t_n$ satisfy $g(t_n + h) \simeq g(t_n) + h g'(t_n) \simeq y_n + h f(g(t_n), t_n) + \simeq y_n + h f(y_n, t_n)$.
In [32]:
let rungekutta1 (f : float array -> float -> float array) (y0 : float array) (t : float array) =
let d = Array.length y0 in
let n = Array.length t in
let y = Array.make_matrix n d 0. in
for j = 0 to d-1 do
y.(0).(j) <- y0.(j);
done;
for i = 0 to n-2 do
let h = t.(i+1) -. t.(i) in
let fyt = f y.(i) t.(i) in
for j = 0 to d-1 do
y.(i+1).(j) <- y.(i).(j) +. h *. fyt.(j);
done;
done;
y
;;
Out[32]:
val rungekutta1 : (float array -> float -> float array) -> float array -> float array -> float array array = <fun>
We have to redefine ourselfves most of the useful functions:
In [28]:
let linspace (a : float) (b : float) (n : int) =
let t = Array.make n a in
let h = (b -. a) /. (float_of_int n) in
for i = 0 to n-1 do
t.(i) <- a +. (float_of_int i) *. h;
done;
t
;;
Out[28]:
val linspace : float -> float -> int -> float array = <fun>
In [30]:
let t = linspace 0. 10. 101 ;;
Out[30]:
val t : float array =
[|0.; 0.0990099009900990146; 0.198019801980198029; 0.297029702970297071;
0.396039603960396058; 0.495049504950495045; 0.594059405940594143;
0.69306930693069313; 0.792079207920792117; 0.891089108910891103;
0.99009900990099009; 1.08910891089108919; 1.18811881188118829;
1.28712871287128716; 1.38613861386138626; 1.48514851485148514;
1.58415841584158423; 1.68316831683168333; 1.78217821782178221;
1.8811881188118813; 1.98019801980198018; 2.07920792079207928;
2.17821782178217838; 2.27722772277227747; 2.37623762376237657;
2.47524752475247523; 2.57425742574257432; 2.67326732673267342;
2.77227722772277252; 2.87128712871287162; 2.97029702970297027;
3.06930693069306937; 3.16831683168316847; 3.26732673267326756;
3.36633663366336666; 3.46534653465346532; 3.56435643564356441;
3.66336633663366351; 3.76237623762376261; 3.86138613861386171;
3.96039603960396036; 4.05940594059405946; 4.15841584158415856;
4.25742574257425765; 4.35643564356435675; 4.45544554455445585;
4.55445544554455495; 4.65346534653465405; 4.75247524752475314;
4.85148514851485135; 4.95049504950495045; 5.04950495049504955;
5.14851485148514865; 5.24752475247524774; 5.34653465346534684;
5.44554455445544594; 5.54455445544554504; 5.64356435643564414;
5.74257425742574323; 5.84158415841584144; 5.94059405940594054;
6.03960396039604; 6.13861386138613874; 6.23762376237623783;
6.33663366336633693; 6.43564356435643603; 6.53465346534653513;
6.63366336633663423; 6.73267326732673332; 6.83168316831683242;
6.93069306930693063; 7.02970297029703; 7.12871287128712883;
7.22772277227722793; 7.32673267326732702; 7.42574257425742612;
7.52475247524752522; 7.62376237623762432; 7.72277227722772341;
7.82178217821782251; 7.92079207920792072; 8.01980198019802;
8.11881188118811892; 8.21782178217821802; 8.31683168316831711;
8.41584158415841621; 8.51485148514851531; 8.61386138613861441;
8.7128712871287135; 8.8118811881188126; 8.9108910891089117;
9.00990099009901; 9.10891089108911; 9.20792079207920899;
9.30693069306930809; 9.40594059405940719; 9.50495049504950629;
9.60396039603960361; 9.70297029702970271; 9.8019801980198018;
9.9009900990099009|]
In [31]:
let sol = rungekutta1 f_pend y0 t ;;
Out[31]:
val sol : float array array =
[|[|3.04156; 0.|]; [|3.04156; -0.0494385678472543874|];
[|3.03666509229235126; -0.097653408767596539|];
[|3.02699643795892603; -0.147085318849099422|];
[|3.01243353510257972; -0.200051306879579782|];
[|2.99262647501549273; -0.258862071967601193|];
[|2.96699656689988878; -0.325927783092890666|];
[|2.93472648936593927; -0.403855500959909908|];
[|2.89474079620159186; -0.495539222697703119|];
[|2.84567750682558174; -0.604239966604331213|];
[|2.78585176755782626; -0.73364756279708776|];
[|2.71321339500365921; -0.887906266173617|];
[|2.62530188350132088; -1.07157049578384078|];
[|2.51920579480985163; -1.28943151630173158|];
[|2.39153930804730397; -1.54611693249596316|];
[|2.23845842364176306; -1.84531096455946542|];
[|2.05575436774478604; -2.18838313254248|];
[|1.83908277046335211; -2.57218281242493196|];
[|1.5844112048767256; -2.98585479258093889|];
[|1.28878201749247401; -3.40695111899184289|];
[|0.951460124522985; -3.79811412259147829|];
[|0.575409221296105611; -4.1072023706189853|];
[|0.168755521234819572; -4.27493407456400387|];
[|-0.254505278226963338; -4.25226525047144666|];
[|-0.675521639659780293; -4.02237420284389824|];
[|-1.07377651122848161; -3.61325381198249129|];
[|-1.43152441340496628; -3.08866463882397024|];
[|-1.7373327934865479; -2.52195643825148741|];
[|-1.98703145073917065; -1.97133136657366337|];
[|-2.18221277416230564; -1.46975502925005852|];
[|-2.3277330740880533; -1.0280107170053443|];
[|-2.42951631339551311; -0.642692161055145705|];
[|-2.49314920062869572; -0.303315136253584139|];
[|-2.52318040223796158; 0.00317609544709579472|];
[|-2.52286593734220954; 0.29009856550137425|];
[|-2.49414330709454868; 0.570045824573125581|];
[|-2.43770312644374432; 0.85452682168123717|];
[|-2.35309651043570112; 1.15376510886464168|];
[|-2.23886234124118211; 1.47634220367394819|];
[|-2.09268984582791973; 1.82842342496740917|];
[|-1.91165782355391944; 2.21231233647640257|];
[|-1.69261699816021594; 2.62411996065747877|];
[|-1.4328031406693762; 3.05054720612826546|];
[|-1.13076876382499325; 3.46538219564712158|];
[|-0.787661615741119592; 3.82749645947745609|];
[|-0.408701570248301838; 4.08360074067417056|];
[|-0.00438466523105685; 4.17926318538280217|];
[|0.409403768965260539; 4.07798669623058441|];
[|0.813164827997992; 3.77998581698302|];
[|1.18742084948145621; 3.32678533270321619|];
[|1.51680553588771549; 2.78532662951843868|];
[|1.79258044970142261; 2.22205475687363219|];
[|2.01258587117405963; 1.68412936790090439|];
[|2.17933135314444648; 1.19492419674890016|];
[|2.29764067955522888; 0.759165492009172826|];
[|2.37280557975415718; 0.370436101655445|];
[|2.40948242150222125; 0.0170768149151984128|];
[|2.41117319525620122; -0.314257086424339|];
[|2.38005863224389058; -0.636766636731305624|];
[|2.31701243058732631; -0.962605357601812539|];
[|2.22170496943863194; -1.30227496505272144|];
[|2.09276685408687735; -1.66386893711239425|];
[|1.92802735536287773; -2.05181197311360775|];
[|1.72487765505459967; -2.46482092024770427|];
[|1.48083597978254944; -2.89299512443404883|];
[|1.19440081894749484; -3.31443400560251211|];
[|0.866239036214572544; -3.69278741132267418|];
[|0.500616520242030294; -3.97855919213080167|];
[|0.106699768545910956; -4.1176871798846113|];
[|-0.300992031442664787; -4.06848572681177778|];
[|-0.703812400433926677; -3.82101440440216811|];
[|-1.08213065829552768; -3.40607400666509186|];
[|-1.41936570846038856; -2.88465631149967239|];
[|-1.70497524425243552; -2.32386963189152063|];
[|-1.9350613464199129; -1.77574835044584578|];
[|-2.1108780147808881; -1.26922692614652322|];
[|-2.236544047072623; -0.813222952841671676|];
[|-2.3170611711163529; -0.403759425509745029|];
[|-2.35703735185989194; -0.0302852883400764328|];
[|-2.36003589525989943; 0.320222039037892758|];
[|-2.32833074287991026; 0.66100133863006616|];
[|-2.26288506578782433; 1.00430942404913726|];
[|-2.16344848914929599; 1.36059620737647924|];
[|-2.02873599336944643; 1.73754293425696282|];
[|-1.85670203948261836; 2.13857666169778726|];
[|-1.64496177594818382; 2.56059536373928109|];
[|-1.39143748250865085; 2.99090290597198161|];
[|-1.0953084819173653; 3.40397871344527214|];
[|-0.758280886526744; 3.75985523320688042|];
[|-0.386017992149824796; 4.00722235732067755|];
[|0.0107366966938070019; 4.0944210195236419|];
[|0.416124916448623372; 3.98775887032095389|];
[|0.810952527371490373; 3.68894350003373317|];
[|1.1761944580678998; 3.23875004120867871|];
[|1.49686277897965048; 2.70157816128681105|];
[|1.76434576524567155; 2.14101030581653617|];
[|1.97632698364334858; 1.60220920537008604|];
[|2.13496155843246349; 1.1076529482996047|];
[|2.24463016717499864; 0.661901511814660393|];
[|2.31016497032496515; ...|]; ...|]
In [34]:
let column sol i =
Array.map (fun x -> x.(i)) sol
;;
Out[34]:
val column : 'a array array -> int -> 'a array = <fun>
Let see a first plot!
In [151]:
let vp = A.init ~w:800. ~h:360. ["jupyter"] in
A.Axes.box vp ;
A.Viewport.xlabel vp "Time t";
A.Viewport.title vp "Solution to the pendulum ODE with Runge-Kutta 1";
A.set_color vp A.Color.red ;
A.Viewport.text vp 7. 3. "o theta(t)" ;
A.Array.xy ~style:(`Linesmarkers "o") vp t (column sol 0);
A.set_color vp A.Color.blue ;
A.Viewport.text vp 7. 2.5 "+ omega(t)" ;
A.Array.xy ~style:(`Linesmarkers "+") vp t (column sol 1);
A.close vp;;
Out[151]:
- : unit = ()
With the same number of points, the Euler method (i.e. the Runge-Kutta method of order 1) is less precise than the reference solve method. With more points, it can give a satisfactory approximation of the solution:
In [55]:
let t = linspace 0. 10. 1001 ;;
let sol = rungekutta1 f_pend y0 t ;;
Out[55]:
val t : float array =
[|0.; 0.00999000999000999; 0.01998001998001998; 0.0299700299700299683;
0.03996003996003996; 0.0499500499500499517; 0.0599400599400599365;
0.0699300699300699352; 0.07992007992007992; 0.0899100899100899;
0.0999000999000999; 0.109890109890109888; 0.119880119880119873;
0.129870129870129858; 0.13986013986013987; 0.149850149850149855;
0.15984015984015984; 0.169830169830169825; 0.17982017982017981;
0.189810189810189822; 0.199800199800199807; 0.209790209790209792;
0.219780219780219777; 0.229770229770229761; 0.239760239760239746;
0.249750249750249759; 0.259740259740259716; 0.269730269730269756;
0.279720279720279741; 0.289710289710289726; 0.29970029970029971;
0.309690309690309695; 0.31968031968031968; 0.329670329670329665;
0.33966033966033965; 0.349650349650349634; 0.359640359640359619;
0.369630369630369604; 0.379620379620379644; 0.389610389610389629;
0.399600399600399614; 0.409590409590409599; 0.419580419580419584;
0.429570429570429568; 0.439560439560439553; 0.449550449550449538;
0.459540459540459523; 0.469530469530469508; 0.479520479520479492;
0.489510489510489533; 0.499500499500499517; 0.509490509490509447;
0.519480519480519431; 0.529470529470529416; 0.539460539460539512;
0.549450549450549497; 0.559440559440559482; 0.569430569430569467;
0.579420579420579451; 0.589410589410589436; 0.599400599400599421;
0.609390609390609406; 0.61938061938061939; 0.629370629370629375;
0.63936063936063936; 0.649350649350649345; 0.65934065934065933;
0.669330669330669314; 0.679320679320679299; 0.689310689310689284;
0.699300699300699269; 0.709290709290709254; 0.719280719280719238;
0.729270729270729223; 0.739260739260739208; 0.749250749250749304;
0.759240759240759289; 0.769230769230769273; 0.779220779220779258;
0.789210789210789243; 0.799200799200799228; 0.809190809190809213;
0.819180819180819197; 0.829170829170829182; 0.839160839160839167;
0.849150849150849152; 0.859140859140859137; 0.869130869130869121;
0.879120879120879106; 0.889110889110889091; 0.899100899100899076;
0.909090909090909061; 0.919080919080919; 0.929070929070929;
0.939060939060939; 0.949050949050949; 0.959040959040959;
0.96903096903096908; 0.979020979020979065; 0.989010989010989;
0.999000999000999; 1.00899100899100902; 1.01898101898101889;
1.02897102897102899; 1.03896103896103886; 1.04895104895104896;
1.05894105894105883; 1.06893106893106893; 1.07892107892107902;
1.0889110889110889; 1.09890109890109899; 1.10889110889110887;
1.11888111888111896; 1.12887112887112884; 1.13886113886113893;
1.14885114885114881; 1.1588411588411589; 1.16883116883116878;
1.17882117882117887; 1.18881118881118875; 1.19880119880119884;
1.20879120879120872; 1.21878121878121881; 1.22877122877122869;
1.23876123876123878; 1.24875124875124865; 1.25874125874125875;
1.26873126873126862; 1.27872127872127872; 1.28871128871128882;
1.29870129870129869; 1.30869130869130879; 1.31868131868131866;
1.32867132867132876; 1.33866133866133863; 1.34865134865134872;
1.3586413586413586; 1.36863136863136869; 1.37862137862137857;
1.38861138861138866; 1.39860139860139854; 1.40859140859140863;
1.41858141858141851; 1.4285714285714286; 1.43856143856143848;
1.44855144855144857; 1.45854145854145845; 1.46853146853146854;
1.47852147852147842; 1.48851148851148851; 1.49850149850149861;
1.50849150849150848; 1.51848151848151858; 1.52847152847152845;
1.53846153846153855; 1.54845154845154842; 1.55844155844155852;
1.56843156843156839; 1.57842157842157849; 1.58841158841158836;
1.59840159840159846; 1.60839160839160833; 1.61838161838161843;
1.6283716283716283; 1.63836163836163839; 1.64835164835164827;
1.65834165834165836; 1.66833166833166824; 1.67832167832167833;
1.68831168831168821; 1.6983016983016983; 1.7082917082917084;
1.71828171828171827; 1.72827172827172837; 1.73826173826173824;
1.74825174825174834; 1.75824175824175821; 1.76823176823176831;
1.77822177822177818; 1.78821178821178828; 1.79820179820179815;
1.80819180819180825; 1.81818181818181812; 1.82817182817182822;
1.83816183816183809; 1.84815184815184819; 1.85814185814185806;
1.86813186813186816; 1.87812187812187803; 1.88811188811188813;
1.898101898101898; 1.9080919080919081; 1.91808191808191797;
1.92807192807192807; 1.93806193806193816; 1.94805194805194803;
1.95804195804195813; 1.968031968031968; 1.9780219780219781;
1.98801198801198797; 1.99800199800199807; 2.00799200799200817;
2.01798201798201804; 2.02797202797202791; 2.03796203796203779;
2.0479520479520481; 2.05794205794205798; 2.06793206793206785;
2.07792207792207773; 2.08791208791208804; 2.09790209790209792;
2.10789210789210779; 2.11788211788211767; 2.12787212787212798;
2.13786213786213786; 2.14785214785214773; 2.15784215784215805;
2.16783216783216792; 2.1778221778221778; 2.18781218781218767;
2.19780219780219799; 2.20779220779220786; 2.21778221778221774;
2.22777222777222761; 2.23776223776223793; 2.2477522477522478;
2.25774225774225767; 2.26773226773226755; 2.27772227772227787;
2.28771228771228774; 2.29770229770229761; 2.30769230769230749;
2.31768231768231781; 2.32767232767232768; 2.33766233766233755;
2.34765234765234787; 2.35764235764235774; 2.36763236763236762;
2.37762237762237749; 2.38761238761238781; 2.39760239760239768;
2.40759240759240756; 2.41758241758241743; 2.42757242757242775;
2.43756243756243762; 2.4475524475524475; 2.45754245754245737;
2.46753246753246769; 2.47752247752247756; 2.48751248751248744;
2.49750249750249731; 2.50749250749250763; 2.5174825174825175;
2.52747252747252737; 2.53746253746253725; 2.54745254745254757;
2.55744255744255744; 2.56743256743256731; 2.57742257742257763;
2.58741258741258751; 2.59740259740259738; 2.60739260739260725;
2.61738261738261757; 2.62737262737262744; 2.63736263736263732;
2.64735264735264719; 2.65734265734265751; 2.66733266733266738;
2.67732267732267726; 2.68731268731268713; 2.69730269730269745;
2.70729270729270732; 2.7172827172827172; 2.72727272727272707;
2.73726273726273739; 2.74725274725274726; 2.75724275724275714;
2.76723276723276701; 2.77722277722277733; 2.7872127872127872;
2.79720279720279708; 2.80719280719280739; 2.81718281718281727;
2.82717282717282714; 2.83716283716283701; 2.84715284715284733;
2.85714285714285721; 2.86713286713286708; 2.87712287712287695;
2.88711288711288727; 2.89710289710289715; 2.90709290709290702;
2.91708291708291689; 2.92707292707292721; 2.93706293706293708;
2.94705294705294696; 2.95704295704295683; 2.96703296703296715;
2.97702297702297702; ...|]
Out[55]:
val sol : float array array =
[|[|3.04156; 0.|]; [|3.04156; -0.00498830704552716622|];
[|3.04151016676278196; -0.00996415578174981824|];
[|3.04141062474698032; -0.0149300540462087225|];
[|3.04126147335790709; -0.0198884971903821491|];
[|3.04106278707228883; -0.0248419693190828489|];
[|3.04081461555061949; -0.0297929445205624109|];
[|3.04051698373722745; -0.0347438880877792081|];
[|3.04016989194813858; -0.0396972577312808828|];
[|3.0397733159468272; -0.0446555047841474054|];
[|3.03932720700792469; -0.0496210753994339171|];
[|3.03883149196896918; -0.0545964117405467325|];
[|3.03828607327026257; -0.059583953164978172|];
[|3.0376908289829; -0.0645861374018181639|];
[|3.0370456128250396; -0.0696054017234524758|];
[|3.03635025416646354; -0.0746441841118481442|];
[|3.03560455802149; -0.0797049244198173|];
[|3.0348083050302832; -0.0847900655276403459|];
[|3.03396125142860829; -0.0899020544954184653|];
[|3.03306312900607677; -0.0950433437115140861|];
[|3.03211364505291492; -0.10021639203742537|];
[|3.03111248229529817; -0.105423665949428341|];
[|3.03005929881928; -0.110667640677305584|];
[|3.02895372798334295; -0.115950801340466647|];
[|3.0277953783196021; -0.121275644081749304|];
[|3.02658383342368031; -0.126644677199174549|];
[|3.02531865183327886; -0.132060422275910555|];
[|3.02399936689545745; -0.137525415308682619|];
[|3.02262548662264363; -0.143042207834846169|];
[|3.02119649353738051; -0.148613368058319567|];
[|3.01971184450582886; -0.154241481974550704|];
[|3.01817097056002925; -0.159929154494668657|];
[|3.01657327670893372; -0.165679010568946294|];
[|3.01491814173821515; -0.171493696309673593|];
[|3.01320491799885781; -0.177375880113513379|];
[|3.01143293118453714; -0.183328253783381329|];
[|3.00960148009779; -0.189353533649860462|];
[|3.00770983640498413; -0.195454461692126963|];
[|3.00575724438008773; -0.201633806658328057|];
[|3.00374292063724724; -0.207894365185315|];
[|3.00166605385217933; -0.214238962917594228|];
[|2.99952580447238315; -0.220670455625315509|];
[|2.99732130441618638; -0.227191730321072205|];
[|2.99505165676063134; -0.233805706375238592|];
[|2.99271593541822156; -0.240515336629518517|];
[|2.99031318480254216; -0.247323608508324466|];
[|2.98784241948277884; -0.254233545127548088|];
[|2.98530262382715916; -0.261248206400220451|];
[|2.98269275163534875; -0.268370690138495627|];
[|2.98001172575983908; -0.275604133151319375|];
[|2.97725843771636933; -0.282951712337071504|];
[|2.97443174728343163; -0.290416645770390724|];
[|2.97153048209092; -0.298002193782306302|];
[|2.96855343719799; -0.305711660032710719|];
[|2.96549937466020053; -0.313548392574114043|];
[|2.9623670230860335; -0.321515784905517|];
[|2.9591550771828814; -0.329617277015135157|];
[|2.95586219729262023; -0.337856356410588754|];
[|2.95248700891689; -0.346236559135054911|];
[|2.94902810223222422; -0.354761470767747611|];
[|2.94548403159518379; -0.363434727406955427|];
[|2.94185331503767156; -0.372260016633721047|];
[|2.9381344337526194; -0.381241078454092552|];
[|2.93432583157026095; -0.390381706217713731|];
[|2.93042591442522893; -0.39968574751034619|];
[|2.92643304981473618; -0.409157105017733391|];
[|2.92234556624812525; -0.418799737358020729|];
[|2.91816175268810518; -0.428617659879740431|];
[|2.91387985798401195; -0.438614945422150737|];
[|2.90949809029747719; -0.448795725034487669|];
[|2.90501461652090898; -0.459164188650442|];
[|2.9004275616892361; -0.469724585713915843|];
[|2.89573500838540099; -0.480481225751838|];
[|2.89093499614012783; -0.491438478889528|];
[|2.88602552082654595; -0.502600776303792296|];
[|2.88100453405028434; -0.513972610608612324|];
[|2.87586994253571282; -0.525558536167942347|];
[|2.87061960750906; -0.537363169329777|];
[|2.86525134407919202; -0.549391188575263811|];
[|2.85976292061690174; -0.56164733457624072|];
[|2.85415205813362283; -0.574136410154152621|];
[|2.84841642966055453; -0.58686328013286|];
[|2.8425536596292571; -0.599832871077383589|];
[|2.83656132325485766; -0.613050170910138359|];
[|2.83043694592308803; -0.62652022839569621|];
[|2.82417800258247187; -0.640248152484572786|];
[|2.81778191714306558; -0.654239111505968|];
[|2.81124606188326576; -0.668498332198795|];
[|2.804567756866295; -0.68303109856970845|];
[|2.79774426936809606; -0.697842750566194492|];
[|2.7907728133184837; -0.712938682552101|];
[|2.7836505487575236; -0.728324341572281098|];
[|2.77637458130924886; -0.744005225392280645|];
[|2.76894196167496043; -0.759986880298231471|];
[|2.76134968514850465; -0.776274898641312|];
[|2.75359469115608402; -0.792874916110307359|];
[|2.74567386282331372; -0.8097926087149383|];
[|2.73758402657241495; -0.827033689461742272|];
[|2.72932195175261727; -0.844603904703366326|];
[|2.72088435030702902; ...|]; ...|]
In [153]:
let vp = A.init ~w:800. ~h:360. ["jupyter"] in
A.Axes.box vp ;
A.Viewport.xlabel vp "Time t";
A.Viewport.title vp "Solution to the pendulum ODE with Runge-Kutta 1";
A.set_color vp A.Color.red ;
A.Viewport.text vp 8.5 3. "o theta(t)" ;
A.Array.xy ~style:(`Linesmarkers "o") vp t (column sol 0);
A.set_color vp A.Color.blue ;
A.Viewport.text vp 8.5 2.5 "+ omega(t)" ;
A.Array.xy vp ~style:(`Linesmarkers "+") t (column sol 1);
A.close vp;;
Out[153]:
- : unit = ()
In [72]:
let t = linspace 0. 10. 10001 ;;
let sol = rungekutta1 f_pend y0 t ;;
Out[72]:
val t : float array =
[|0.; 0.000999900009999; 0.001999800019998; 0.00299970002999700022;
0.003999600039996; 0.004999500049995; 0.00599940005999400044;
0.00699930006999300094; 0.007999200079992; 0.008999100089991;
0.00999900009999; 0.0109989001099890012; 0.0119988001199880009;
0.012998700129987; 0.0139986001399860019; 0.0149985001499850015;
0.015998400159984; 0.016998300169983; 0.017998200179982;
0.018998100189981; 0.01999800019998; 0.0209979002099790028;
0.0219978002199780025; 0.0229977002299770021; 0.0239976002399760018;
0.024997500249975; 0.025997400259974; 0.026997300269973;
0.0279972002799720038; 0.0289971002899710034; 0.0299970002999700031;
0.0309969003099690027; 0.031996800319968; 0.032996700329967;
0.033996600339966; 0.034996500349965; 0.035996400359964;
0.036996300369963; 0.037996200379962; 0.038996100389961;
0.03999600039996; 0.040995900409959006; 0.0419958004199580057;
0.0429957004299570053; 0.0439956004399560049; 0.0449955004499550046;
0.0459954004599540042; 0.0469953004699530039; 0.0479952004799520035;
0.048995100489951; 0.04999500049995; 0.050994900509949;
0.051994800519948; 0.052994700529947; 0.053994600539946;
0.054994500549945; 0.0559944005599440076; 0.0569943005699430072;
0.0579942005799420068; 0.0589941005899410065; 0.0599940005999400061;
0.0609939006099390058; 0.0619938006199380054; 0.062993700629937;
0.063993600639936; 0.064993500649935; 0.065993400659934;
0.066993300669933; 0.067993200679932; 0.068993100689931;
0.06999300069993; 0.070992900709929; 0.071992800719928;
0.072992700729927; 0.073992600739926; 0.074992500749925;
0.075992400759924; 0.076992300769923; 0.077992200779922;
0.078992100789921; 0.07999200079992; 0.0809919008099190124;
0.0819918008199180121; 0.0829917008299170117; 0.0839916008399160113;
0.084991500849915011; 0.0859914008599140106; 0.0869913008699130103;
0.0879912008799120099; 0.0889911008899110095; 0.0899910008999100092;
0.0909909009099090088; 0.0919908009199080084; 0.0929907009299070081;
0.0939906009399060077; 0.0949905009499050074; 0.095990400959904007;
0.096990300969903; 0.097990200979902; 0.098990100989901; 0.0999900009999;
0.100989901009899; 0.101989801019898; 0.102989701029897;
0.103989601039896; 0.104989501049895; 0.105989401059894;
0.106989301069893; 0.107989201079892; 0.108989101089891;
0.10998900109989; 0.110988901109889; 0.111988801119888015;
0.112988701129887015; 0.113988601139886014; 0.114988501149885014;
0.115988401159884014; 0.116988301169883013; 0.117988201179882013;
0.118988101189881013; 0.119988001199880012; 0.120987901209879012;
0.121987801219878012; 0.122987701229877011; 0.123987601239876011;
0.12498750124987501; 0.125987401259874; 0.126987301269873;
0.127987201279872; 0.128987101289871; 0.12998700129987;
0.130986901309869; 0.131986801319868; 0.132986701329867;
0.133986601339866; 0.134986501349865; 0.135986401359864;
0.136986301369863; 0.137986201379862; 0.138986101389861;
0.13998600139986; 0.140985901409859; 0.141985801419858;
0.142985701429857; 0.143985601439856; 0.144985501449855;
0.145985401459854; 0.146985301469853; 0.147985201479852;
0.148985101489851; 0.14998500149985; 0.150984901509849;
0.151984801519848; 0.152984701529847; 0.153984601539846;
0.154984501549845; 0.155984401559844; 0.156984301569843;
0.157984201579842; 0.158984101589841; 0.15998400159984;
0.160983901609839025; 0.161983801619838025; 0.162983701629837024;
0.163983601639836024; 0.164983501649835024; 0.165983401659834023;
0.166983301669833023; 0.167983201679832023; 0.168983101689831022;
0.169983001699830022; 0.170982901709829022; 0.171982801719828021;
0.172982701729827021; 0.173982601739826021; 0.17498250174982502;
0.17598240175982402; 0.176982301769823019; 0.177982201779822019;
0.178982101789821019; 0.179982001799820018; 0.180981901809819018;
0.181981801819818018; 0.182981701829817017; 0.183981601839816017;
0.184981501849815017; 0.185981401859814016; 0.186981301869813016;
0.187981201879812015; 0.188981101889811015; 0.189981001899810015;
0.190980901909809014; 0.191980801919808014; 0.192980701929807;
0.193980601939806; 0.194980501949805; 0.195980401959804;
0.196980301969803; 0.197980201979802; 0.198980101989801; 0.1999800019998;
0.200979902009799; 0.201979802019798; 0.202979702029797;
0.203979602039796; 0.204979502049795; 0.205979402059794;
0.206979302069793; 0.207979202079792; 0.208979102089791;
0.20997900209979; 0.210978902109789; 0.211978802119788;
0.212978702129787; 0.213978602139786; 0.214978502149785;
0.215978402159784; 0.216978302169783; 0.217978202179782;
0.218978102189781; 0.21997800219978; 0.220977902209779;
0.221977802219778; 0.222977702229777; 0.22397760223977603;
0.22497750224977503; 0.22597740225977403; 0.226977302269773029;
0.227977202279772029; 0.228977102289771028; 0.229977002299770028;
0.230976902309769028; 0.231976802319768027; 0.232976702329767027;
0.233976602339766027; 0.234976502349765026; 0.235976402359764026;
0.236976302369763026; 0.237976202379762025; 0.238976102389761025;
0.239976002399760024; 0.240975902409759024; 0.241975802419758024;
0.242975702429757023; 0.243975602439756023; 0.244975502449755023;
0.245975402459754022; 0.246975302469753022; 0.247975202479752022;
0.248975102489751021; 0.249975002499750021; 0.250974902509749;
0.251974802519748; 0.252974702529747; 0.253974602539746;
0.254974502549745; 0.255974402559744; 0.256974302569743;
0.257974202579742; 0.258974102589741; 0.25997400259974;
0.260973902609739; 0.261973802619738; 0.262973702629737;
0.263973602639736; 0.264973502649735; 0.265973402659734;
0.266973302669733; 0.267973202679732; 0.268973102689731;
0.26997300269973; 0.270972902709729; 0.271972802719728;
0.272972702729727; 0.273972602739726; 0.274972502749725;
0.275972402759724; 0.276972302769723; 0.277972202779722;
0.278972102789721; 0.27997200279972; 0.280971902809719;
0.281971802819718; 0.282971702829717; 0.283971602839716;
0.284971502849715; 0.285971402859714; 0.286971302869713;
0.287971202879712; 0.288971102889711; 0.28997100289971;
0.290970902909709; 0.291970802919708; 0.292970702929707;
0.293970602939706; 0.294970502949705; 0.295970402959704;
0.296970302969703; 0.297970202979702; ...|]
Out[72]:
val sol : float array array =
[|[|3.04156; 0.|]; [|3.04156; -0.000499279607296539596|];
[|3.04155950077031578; -0.000998434407171997374|];
[|3.04155850243574211; -0.00149746691424664175|];
[|3.04155700511855942; -0.00199637964189098346|];
[|3.04155500893853548; -0.00249517510223822188|];
[|3.04155251401292581; -0.00299385580619668289|];
[|3.04154952045647509; -0.00349242426346224686|];
[|3.0415460283814193; -0.00399088298253077144|];
[|3.04154203789748534; -0.0044892344707105|];
[|3.04153754911189322; -0.00498748123413446723|];
[|3.04153256212935741; -0.00548562577777289128|];
[|3.04152707705208725; -0.00598367060544555563|];
[|3.0415210939797892; -0.00648161821983419227|];
[|3.04151461300966641; -0.00697947112249484164|];
[|3.04150763423642134; -0.00747723181387021513|];
[|3.0415001577522558; -0.00797490279330204457|];
[|3.04149218364687313; -0.00847248655904342569|];
[|3.0414837120074778; -0.00896998560827114676|];
[|3.04147474291877851; -0.009467402437098019|];
[|3.04146527646298681; -0.00996473954058518796|];
[|3.04145531271982072; -0.0104619994127544448|];
[|3.04144485176650337; -0.0109591845466005186|];
[|3.04143389367776562; -0.0114562974341033807|];
[|3.04142243852584659; -0.0119533405662405161|];
[|3.04141048638049494; -0.0124503164329992054|];
[|3.04139803730896929; -0.0129472275233887858|];
[|3.04138509137603918; -0.0134440763254529139|];
[|3.04137164864398679; -0.0139408653262818177|];
[|3.04135770917260784; -0.0144375970120245306|];
[|3.04134327301921115; -0.0149342738679011343|];
[|3.04132834023862131; -0.0154308983782149838|];
[|3.04131291088317868; -0.0159274730263649252|];
[|3.04129698500274026; -0.0164240002948575076|];
[|3.04128056264468105; -0.0169204826653191823|];
[|3.04126364385389492; -0.0174169226185085027|];
[|3.04124622867279459; -0.017913322634328311|];
[|3.04122831714131348; -0.0184096851918379191|];
[|3.04120990929690604; -0.0189060127692652763|];
[|3.04119100517454921; -0.0194023078440191399|];
[|3.04117160480674187; -0.0198985728927012287|];
[|3.04115170822350755; -0.0203948103911183801|];
[|3.04113131545239357; -0.0208910228142946755|];
[|3.04111042651847274; -0.0213872126364835977|];
[|3.04108904144434389; -0.0218833823311801422|];
[|3.04106716025013224; -0.0223795343711329434|];
[|3.04104478295349079; -0.0228756712283563896|];
[|3.04102190956960072; -0.0233717953741427245|];
[|3.04099854011117232; -0.0238679092790741523|];
[|3.04097467458844539; -0.0243640154130349181|];
[|3.04095031300919; -0.0248601162452234056|];
[|3.04092545537870818; -0.0253562142441642047|];
[|3.04090010169983183; -0.0258523118777201875|];
[|3.0408742519729266; -0.0263484116131045693|];
[|3.04084790619589107; -0.0268445159168929672|];
[|3.04082106436415733; -0.0273406272550354504|];
[|3.04079372647069146; -0.0278367480928685788|];
[|3.04076589250599527; -0.0283328808951274387|];
[|3.040737562458105; -0.0288290281259576779|];
[|3.04070873631259353; -0.0293251922489275313|];
[|3.04067941405257081; -0.0298213757270398248|];
[|3.04064959565868298; -0.0303175810227439939|];
[|3.04061928110911506; -0.0308138105979480874|];
[|3.04058847037959; -0.0313100669140307583|];
[|3.04055716344336968; -0.0318063524318532506|];
[|3.0405253602712552; -0.0323026696117714|];
[|3.04049306083158744; -0.0327990209136475863|];
[|3.04046026509024792; -0.0332954087968627271|];
[|3.04042697301065923; -0.0337918357203282269|];
[|3.0403931845537846; -0.0342883041424979496|];
[|3.04035889967812967; -0.0347848165213801436|];
[|3.04032411833974203; -0.0352813753145494249|];
[|3.04028884049221215; -0.0357779829791586845|];
[|3.04025306608667334; -0.0362746419719510368|];
[|3.04021679507180309; -0.0367713547492717549|];
[|3.04018002739382176; -0.0372681237670801635|];
[|3.04014276299649433; -0.037764951480961588|];
[|3.04010500182113086; -0.0382618403461392406|];
[|3.04006674380658604; -0.0387587928174861415|];
[|3.04002798888926; -0.039255811349537|];
[|3.03998873700309913; -0.0397528983965001123|];
[|3.0399489880795949; -0.0402500564122692667|];
[|3.03990874204778594; -0.0407472878504355693|];
[|3.03986799883425673; -0.0412445951642993955|];
[|3.03982675836313954; -0.0417419808068822|];
[|3.03978502055611344; -0.0422394472309384067|];
[|3.03974278533240483; -0.0427369968889672602|];
[|3.03970005260878828; -0.0432346322332246707|];
[|3.03965682229958611; -0.0437323557157350828|];
[|3.03961309431666882; -0.0442301697883033|];
[|3.03956886856945507; -0.0447280769025263092|];
[|3.03952414496491308; -0.0452260795098051518|];
[|3.03947892340755921; -0.0457241800613567|];
[|3.03943320379945847; -0.0462223810082255174|];
[|3.03938698604022628; -0.0467206848012956533|];
[|3.03934027002702623; -0.0472190938913024591|];
[|3.03929305565457231; -0.0477176107288443882|];
[|3.03924534281512759; -0.0482162377643948|];
[|3.03919713139850467; -0.0487149774483137538|];
[|3.03914842129206697; ...|]; ...|]
In [154]:
let vp = A.init ~w:800. ~h:360. ["jupyter"] in
A.Axes.box vp ;
A.Viewport.xlabel vp "Time t";
A.Viewport.title vp "Solution to the pendulum ODE with Runge-Kutta 1";
A.set_color vp A.Color.red ;
A.Viewport.text vp 8.5 3. "o theta(t)" ;
A.Array.xy ~style:(`Linesmarkers "o") vp t (column sol 0);
A.set_color vp A.Color.blue ;
A.Viewport.text vp 8.5 2.5 "+ omega(t)" ;
A.Array.xy vp ~style:(`Linesmarkers "+") t (column sol 1);
A.close vp;;
Out[154]:
- : unit = ()
Runge-Kutta method of order 2¶
The order 2 Runge-Method uses this update: $$ y_{n+1} = y_n + h f(t + \frac{h}{2}, y_n + \frac{h}{2} f(t, y_n)),$$ if $h = t_{n+1} - t_n$.
Again, we need some basic functions as OCaml Array are quite poor.
With Julia arrays or NumPy Python arrays, we can write h * f(y, t) to multiply each entry of f(y, t) by a number h.
In [92]:
let aplus a k =
Array.map ( ( +. ) k) a
;;
let ( ++. ) = aplus ;;
Out[92]:
val aplus : float array -> float -> float array = <fun>
Out[92]:
val ( ++. ) : float array -> float -> float array = <fun>
In [93]:
let aaplus a b =
Array.map2 ( +. ) a b
;;
let ( +++. ) = aaplus ;;
Out[93]:
val aaplus : float array -> float array -> float array = <fun>
Out[93]:
val ( +++. ) : float array -> float array -> float array = <fun>
In [94]:
let amul a k =
Array.map ( ( *. ) k ) a
;;
let ( **. ) = amul ;;
Out[94]:
val amul : float array -> float -> float array = <fun>
Out[94]:
val ( **. ) : float array -> float -> float array = <fun>
In [95]:
let rungekutta2 (f : float array -> float -> float array) (y0 : float array) (t : float array) =
let d = Array.length y0 in
let n = Array.length t in
let y = Array.make_matrix n d 0. in
for j = 0 to d-1 do
y.(0).(j) <- y0.(j);
done;
for i = 0 to n-2 do
let h = t.(i+1) -. t.(i) in
let fyt = f y.(i) t.(i) in
let fyt2 = f (y.(i) +++. (fyt **. (h /. 2.))) (t.(i) +. (h /. 2.)) in
for j = 0 to d-1 do
y.(i+1).(j) <- y.(i).(j) +. h *. fyt2.(j);
done;
done;
y
;;
Out[95]:
val rungekutta2 : (float array -> float -> float array) -> float array -> float array -> float array array = <fun>
For our simple ODE example, this method is already quite efficient.
In [87]:
let t = linspace 0. 10. 101 ;;
let sol = rungekutta2 f_pend y0 t ;;
Out[87]:
val t : float array =
[|0.; 0.0990099009900990146; 0.198019801980198029; 0.297029702970297071;
0.396039603960396058; 0.495049504950495045; 0.594059405940594143;
0.69306930693069313; 0.792079207920792117; 0.891089108910891103;
0.99009900990099009; 1.08910891089108919; 1.18811881188118829;
1.28712871287128716; 1.38613861386138626; 1.48514851485148514;
1.58415841584158423; 1.68316831683168333; 1.78217821782178221;
1.8811881188118813; 1.98019801980198018; 2.07920792079207928;
2.17821782178217838; 2.27722772277227747; 2.37623762376237657;
2.47524752475247523; 2.57425742574257432; 2.67326732673267342;
2.77227722772277252; 2.87128712871287162; 2.97029702970297027;
3.06930693069306937; 3.16831683168316847; 3.26732673267326756;
3.36633663366336666; 3.46534653465346532; 3.56435643564356441;
3.66336633663366351; 3.76237623762376261; 3.86138613861386171;
3.96039603960396036; 4.05940594059405946; 4.15841584158415856;
4.25742574257425765; 4.35643564356435675; 4.45544554455445585;
4.55445544554455495; 4.65346534653465405; 4.75247524752475314;
4.85148514851485135; 4.95049504950495045; 5.04950495049504955;
5.14851485148514865; 5.24752475247524774; 5.34653465346534684;
5.44554455445544594; 5.54455445544554504; 5.64356435643564414;
5.74257425742574323; 5.84158415841584144; 5.94059405940594054;
6.03960396039604; 6.13861386138613874; 6.23762376237623783;
6.33663366336633693; 6.43564356435643603; 6.53465346534653513;
6.63366336633663423; 6.73267326732673332; 6.83168316831683242;
6.93069306930693063; 7.02970297029703; 7.12871287128712883;
7.22772277227722793; 7.32673267326732702; 7.42574257425742612;
7.52475247524752522; 7.62376237623762432; 7.72277227722772341;
7.82178217821782251; 7.92079207920792072; 8.01980198019802;
8.11881188118811892; 8.21782178217821802; 8.31683168316831711;
8.41584158415841621; 8.51485148514851531; 8.61386138613861441;
8.7128712871287135; 8.8118811881188126; 8.9108910891089117;
9.00990099009901; 9.10891089108911; 9.20792079207920899;
9.30693069306930809; 9.40594059405940719; 9.50495049504950629;
9.60396039603960361; 9.70297029702970271; 9.8019801980198018;
9.9009900990099009|]
Out[87]:
val sol : float array array =
[|[|3.04156; 0.|]; [|3.03911254614617565; -0.0488267043837982695|];
[|3.03183092241293961; -0.0988404593793910102|];
[|3.01948127772376; -0.152388585215509798|];
[|3.00159483065725174; -0.211923116134953959|];
[|2.97745224188518653; -0.280106126962090252|];
[|2.9460575795158439; -0.359915119550046958|];
[|2.90610187009085497; -0.454746846467090182|];
[|2.85591658999954312; -0.568512488459194332|];
[|2.79341852588212358; -0.705707056702320923|];
[|2.71604980667654772; -0.871417812472915942|];
[|2.62072159504783; -1.07120512144372437|];
[|2.50377854277070533; -1.31073777557490101|];
[|2.3610159594913549; -1.59498872770536404|];
[|2.18780513934095255; -1.92670587301033902|];
[|1.97941455921401221; -2.30382082261406795|];
[|1.73164679588694903; -2.71560548610054786|];
[|1.44191154148313938; -3.13804848237782785|];
[|1.11075480950422323; -3.53044584205407341|];
[|0.743572352514064394; -3.83725457909199763|];
[|0.351758622387601194; -3.99957876824615033|];
[|-0.0477823187833688046; -3.97551260420907848|];
[|-0.435355379078979254; -3.75873893055921959|];
[|-0.792566330837187083; -3.38121852743959517|];
[|-1.10574407550246878; -2.8973969458654607|];
[|-1.36716000728090981; -2.3623768246026|];
[|-1.57416289561318345; -1.81743549025345841|];
[|-1.72737271016564775; -1.2862384265376352|];
[|-1.82893933061530634; -0.777756413931331325|];
[|-1.88129651085227101; -0.291240213889777222|];
[|-1.886439820040684; 0.179192594251308|];
[|-1.84562089865235701; 0.640925073236255471|];
[|-1.75936063343311178; 1.09985381280031236|];
[|-1.62773895223768394; 1.55755230771725195|];
[|-1.45096675062721703; 2.00830812610837084|];
[|-1.23025362332544352; 2.43632820045757681|];
[|-0.968918385066261267; 2.81433103704307896|];
[|-0.67351953612175719; 3.1056585677510693|];
[|-0.354547903212739823; 3.27179885370290968|];
[|-0.0261084401934915111; 3.28443628530758902|];
[|0.295698398658890971; 3.13662613817199531|];
[|0.595270261290242297; 2.84596832523174781|];
[|0.859819841409028673; 2.44743139225582418|];
[|1.0805709012737108; 1.98083632592742709|];
[|1.25264496742762099; 1.48035044177632|];
[|1.37412283593078266; 0.969698968223488889|];
[|1.44490944640484531; 0.462413620202401954|];
[|1.46581287840047891; -0.0352383219593074282|];
[|1.43799464484600303; -0.520516426851603708|];
[|1.3628045813203864; -0.990557587736249801|];
[|1.24196416474831528; -1.43913433247567824|];
[|1.07804479128885755; -1.85411054146652865|];
[|0.875149593837422; -2.216241424541467|];
[|0.639622619355808553; -2.50035497682177432|];
[|0.380498291815628276; -2.67976599575457408|];
[|0.109356994580230626; -2.73341990183247763|];
[|-0.160603905511440936; -2.65307922813173436|];
[|-0.416114927732840034; -2.44671882330057|];
[|-0.6454600443102817; -2.13603282352703161|];
[|-0.839588192768512376; -1.74961989675213392|];
[|-0.992431423756186204; -1.31569413075129171|];
[|-1.10056449840736703; -0.857392350245871682|];
[|-1.16255675909603573; -0.391455573074287411|];
[|-1.17834164809000641; 0.0705805769255422|];
[|-1.14879578161057139; 0.519942157737366184|];
[|-1.07559608109481819; 0.947893575800511368|];
[|-0.961343337819828614; 1.34339665772494898|];
[|-0.809884806196925711; 1.69172414498315193|];
[|-0.626711904469132; 1.9746365504195329|];
[|-0.419249779616737039; 2.17258967945957693|];
[|-0.196827758700094557; 2.26880258893402331|];
[|0.02981868916202704; 2.25392822334985476|];
[|0.249487332939835216; 2.12924443807099406|];
[|0.451643447796762743; 1.90667949135001069|];
[|0.627391080054349626; 1.60566258030707587|];
[|0.770013347948908; 1.24849018107622145|];
[|0.875035674223142; 0.856314657708687643|];
[|0.939958944379134076; 0.446994478711824961|];
[|0.963877532861861286; 0.034882683497337641|];
[|0.947157913711899502; -0.36806262007216034|];
[|0.891272978842545; -0.750618111970555835|];
[|0.798810509166694938; -1.10076437142558659|];
[|0.673612570255083; -1.40484595188823835|];
[|0.520952319112750306; -1.64786573966020122|];
[|0.347619051432104809; -1.81511348099665915|];
[|0.161780333936247855; -1.89485037549946722|];
[|-0.0274542642815124305; -1.88116723211021308|];
[|-0.210730571520729909; -1.7757989665522389|];
[|-0.379249920380139605; -1.58801466145910841|];
[|-0.525459966123299704; -1.33263364367163462|];
[|-0.643477740530508768; -1.02711420011423016|];
[|-0.729209638217741; -0.688915462994294514|];
[|-0.780246110816953; -0.333918707582131769|];
[|-0.79565835672033125; 0.023947572411301854|];
[|-0.775810417760952831; 0.372200606965266079|];
[|-0.722252549738436445; 0.698921727954476935|];
[|-0.637707549933894224; 0.992005726665617704|];
[|-0.52611412873334; 1.23895677403037396|];
[|-0.392656274012012751; 1.42748628240091024|];
[|-0.243692590313146967; ...|]; ...|]
In [156]:
let vp = A.init ~w:800. ~h:360. ["jupyter"] in
A.Axes.box vp ;
A.Viewport.xlabel vp "Time t";
A.Viewport.title vp "Solution to the pendulum ODE with Runge-Kutta 2";
A.set_color vp A.Color.red ;
A.Viewport.text vp 8.5 3. "o theta(t)" ;
A.Array.xy ~style:(`Linesmarkers "o") vp t (column sol 0);
A.set_color vp A.Color.blue ;
A.Viewport.text vp 8.5 2.5 "+ omega(t)" ;
A.Array.xy vp ~style:(`Linesmarkers "+") t (column sol 1);
A.close vp;;
Out[156]:
- : unit = ()
Runge-Kutta method of order 4, "RK4"¶
The order 4 Runge-Method uses this update: $$ y_{n+1} = y_n + \frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4),$$ if $h = t_{n+1} - t_n$, and $$\begin{cases} k_1 &= f(y_n, t_n), \\ k_2 &= f(y_n + \frac{h}{2} k_1, t_n + \frac{h}{2}), \\ k_3 &= f(y_n + \frac{h}{2} k_2, t_n + \frac{h}{2}), \\ k_4 &= f(y_n + h k_3, t_n + h). \end{cases}$$
In [97]:
let rungekutta4 (f : float array -> float -> float array) (y0 : float array) (t : float array) =
let d = Array.length y0 in
let n = Array.length t in
let y = Array.make_matrix n d 0. in
for j = 0 to d-1 do
y.(0).(j) <- y0.(j);
done;
for i = 0 to n-2 do
let h = t.(i+1) -. t.(i) in
let k1 = f y.(i) t.(i) in
let k2 = f (y.(i) +++. (k1 **. (h /. 2.))) (t.(i) +. (h /. 2.)) in
let k3 = f (y.(i) +++. (k2 **. (h /. 2.))) (t.(i) +. (h /. 2.)) in
let k4 = f (y.(i) +++. (k3 **. h)) (t.(i) +. h) in
for j = 0 to d-1 do
y.(i+1).(j) <- y.(i).(j) +. (h/.6.) *. (k1.(j) +. 2.*.k2.(j) +. 2.*.k3.(j) +. k4.(j));
done;
done;
y
;;
Out[97]:
val rungekutta4 : (float array -> float -> float array) -> float array -> float array -> float array array = <fun>
For our simple ODE example, this method is even more efficient.
In [107]:
let t = linspace 0. 10. 31 ;;
let sol = rungekutta4 f_pend y0 t ;;
Out[107]:
val t : float array =
[|0.; 0.322580645161290314; 0.645161290322580627; 0.967741935483871;
1.29032258064516125; 1.61290322580645151; 1.93548387096774199;
2.25806451612903203; 2.58064516129032251; 2.90322580645161299;
3.22580645161290303; 3.54838709677419351; 3.87096774193548399;
4.19354838709677402; 4.51612903225806406; 4.83870967741935498;
5.16129032258064502; 5.48387096774193505; 5.80645161290322598;
6.12903225806451601; 6.45161290322580605; 6.77419354838709697;
7.09677419354838701; 7.41935483870967705; 7.74193548387096797;
8.06451612903225801; 8.38709677419354804; 8.70967741935483808;
9.03225806451612812; 9.35483870967742; 9.67741935483871|]
Out[107]:
val sol : float array array =
[|[|3.04156; 0.|]; [|3.01514459001706; -0.168073982292008067|];
[|2.92524068186212505; -0.40997602608111805|];
[|2.73082595162226482; -0.83986111238677208|];
[|2.34558387188495843; -1.62374367239927175|];
[|1.63147992840831235; -2.86279243145111817|];
[|0.508849624613306517; -3.94638176996873646|];
[|-0.729938057166801491; -3.4423574128899137|];
[|-1.58130089343988933; -1.76934614226155151|];
[|-1.88100723824819749; -0.124998397630503355|];
[|-1.67675988248597196; 1.38189318693121943|];
[|-0.998660995093421544; 2.76245745655193087|];
[|0.00990026972481095; 3.26207897518789913|];
[|0.9349425506047504; 2.28077702467475563|];
[|1.41418850276070662; 0.663566385977470086|];
[|1.36910538191439057; -0.919824478454474903|];
[|0.846941901653241791; -2.23857035893071554|];
[|0.0178465256173153675; -2.70513906586802078|];
[|-0.754166736888580491; -1.91125086669990418|];
[|-1.14564076405576865; -0.47621471666032722|];
[|-1.06175338137134245; 0.966682274531973|];
[|-0.561766705543892408; 2.02888883753261595|];
[|0.146152092751127172; 2.18268237193314318|];
[|0.731805215477561388; 1.32587845081367828|];
[|0.954379467760988254; 0.0337641526393033242|];
[|0.763316350305066615; -1.16858571345420392|];
[|0.257276461206502716; -1.84530056233370932|];
[|-0.327139986143612616; -1.63299256455690345|];
[|-0.715665389025514687; -0.703045696222519889|];
[|-0.759083554196956722; 0.425753036314224231|];
[|-0.467540082468857643; 1.30687972996076196|]|]
In [157]:
let vp = A.init ~w:800. ~h:360. ["jupyter"] in
A.Axes.box vp ;
A.Viewport.xlabel vp "Time t";
A.Viewport.title vp "Solution to the pendulum ODE with Runge-Kutta 4 (31 points)";
A.set_color vp A.Color.red ;
A.Viewport.text vp 8.5 3. "o theta(t)" ;
A.Array.xy ~style:(`Linesmarkers "o") vp t (column sol 0);
A.set_color vp A.Color.blue ;
A.Viewport.text vp 8.5 2.5 "+ omega(t)" ;
A.Array.xy vp ~style:(`Linesmarkers "+") t (column sol 1);
A.close vp;;
Out[157]:
- : unit = ()
In [109]:
let t = linspace 0. 10. 101 ;;
let sol = rungekutta4 f_pend y0 t ;;
Out[109]:
val t : float array =
[|0.; 0.0990099009900990146; 0.198019801980198029; 0.297029702970297071;
0.396039603960396058; 0.495049504950495045; 0.594059405940594143;
0.69306930693069313; 0.792079207920792117; 0.891089108910891103;
0.99009900990099009; 1.08910891089108919; 1.18811881188118829;
1.28712871287128716; 1.38613861386138626; 1.48514851485148514;
1.58415841584158423; 1.68316831683168333; 1.78217821782178221;
1.8811881188118813; 1.98019801980198018; 2.07920792079207928;
2.17821782178217838; 2.27722772277227747; 2.37623762376237657;
2.47524752475247523; 2.57425742574257432; 2.67326732673267342;
2.77227722772277252; 2.87128712871287162; 2.97029702970297027;
3.06930693069306937; 3.16831683168316847; 3.26732673267326756;
3.36633663366336666; 3.46534653465346532; 3.56435643564356441;
3.66336633663366351; 3.76237623762376261; 3.86138613861386171;
3.96039603960396036; 4.05940594059405946; 4.15841584158415856;
4.25742574257425765; 4.35643564356435675; 4.45544554455445585;
4.55445544554455495; 4.65346534653465405; 4.75247524752475314;
4.85148514851485135; 4.95049504950495045; 5.04950495049504955;
5.14851485148514865; 5.24752475247524774; 5.34653465346534684;
5.44554455445544594; 5.54455445544554504; 5.64356435643564414;
5.74257425742574323; 5.84158415841584144; 5.94059405940594054;
6.03960396039604; 6.13861386138613874; 6.23762376237623783;
6.33663366336633693; 6.43564356435643603; 6.53465346534653513;
6.63366336633663423; 6.73267326732673332; 6.83168316831683242;
6.93069306930693063; 7.02970297029703; 7.12871287128712883;
7.22772277227722793; 7.32673267326732702; 7.42574257425742612;
7.52475247524752522; 7.62376237623762432; 7.72277227722772341;
7.82178217821782251; 7.92079207920792072; 8.01980198019802;
8.11881188118811892; 8.21782178217821802; 8.31683168316831711;
8.41584158415841621; 8.51485148514851531; 8.61386138613861441;
8.7128712871287135; 8.8118811881188126; 8.9108910891089117;
9.00990099009901; 9.10891089108911; 9.20792079207920899;
9.30693069306930809; 9.40594059405940719; 9.50495049504950629;
9.60396039603960361; 9.70297029702970271; 9.8019801980198018;
9.9009900990099009|]
Out[109]:
val sol : float array array =
[|[|3.04156; 0.|]; [|3.0391226684992545; -0.0492285629019174747|];
[|3.03177276128229733; -0.0996338062721686679|];
[|3.01927516015746411; -0.153620003759229684|];
[|3.00115241209564187; -0.213696393492316|];
[|2.9766689247853475; -0.282586378652772197|];
[|2.94480437110224091; -0.363336300206605733|];
[|2.90421640674631698; -0.459421802365044596|];
[|2.85319326842347198; -0.574843794764229199|];
[|2.78959810389226659; -0.714194965104298696|];
[|2.71080965865368428; -0.882657987552806|];
[|2.61366933385879374; -1.0858624704672728|];
[|2.49445437940516; -1.32947272354123558|];
[|2.34891337747415152; -1.61829970140246937|];
[|2.17242504413163218; -1.95464379348485062|];
[|1.96037272241410943; -2.33555053709422911|];
[|1.70885073221755923; -2.74888146198779859|];
[|1.41579566533549261; -3.16888687934279778|];
[|1.08249517159997977; -3.55354957515189|];
[|0.715109711947466087; -3.84770377368949612|];
[|0.325441250106985047; -3.99536686046597067|];
[|-0.0699190864919839727; -3.95904774800841208|];
[|-0.452287687444798803; -3.73554535387515418|];
[|-0.804467169392459214; -3.35673269485760128|];
[|-1.11357474679483404; -2.87467196041220285|];
[|-1.37205587780615756; -2.34200178397280734|];
[|-1.57700879364184265; -1.79866877944456216|];
[|-1.72866651635680224; -1.26803574155739063|];
[|-1.82883871793355479; -0.759255125090063|];
[|-1.87972614160690199; -0.27186408672986262|];
[|-1.88319494772193252; 0.199754224892532639|];
[|-1.84045392115182849; 0.662711676307255804|];
[|-1.75206143877875142; 1.12254614316519419|];
[|-1.6182252493199738; 1.58031278585125734|];
[|-1.43939291209877784; 2.02962670364174791|];
[|-1.21712505674593463; 2.4540415503089843|];
[|-0.955162531422122152; 2.82604474993697741|];
[|-0.660430305711281362; 3.10967848132059599|];
[|-0.34353383787510211; 3.26823308268441926|];
[|-0.0182828421951334308; 3.27563108435625772|];
[|0.29988995231147525; 3.12629814733551825|];
[|0.596151233707482575; 2.8375005425207136|];
[|0.858257214380308; 2.44260762434597423|];
[|1.07759690578971656; 1.97974392787833819|];
[|1.24915018030287173; 1.48208982599572381|];
[|1.37072155712265542; 0.973166564272889056|];
[|1.44194967693485165; 0.466767478570425|];
[|1.46345948781373059; -0.0304287144576129154|];
[|1.43633119759964156; -0.515294348898755583|];
[|1.36192787952049943; -0.984670590846349447|];
[|1.24206138750156536; -1.4321665633009637|];
[|1.07944285254233874; -1.84568988317779|];
[|0.878316056391732269; -2.20636299151447624|];
[|0.645090032944561353; -2.48978191303466367|];
[|0.388698736892993613; -2.67030320443958669|];
[|0.120403804144515847; -2.72774383108151941|];
[|-0.14708801964373408; -2.65398358815839597|];
[|-0.401014093611074762; -2.45610407390851693|];
[|-0.629997905437077121; -2.15426998471919617|];
[|-0.825063506953729098; -1.77565126645040539|];
[|-0.979986860814240912; -1.34767078914604421|];
[|-1.0910679227066189; -0.89337466666507992|];
[|-1.1565949674863385; -0.429847881471155435|];
[|-1.17628143989641787; 0.0309220266056009496|];
[|-1.15086377377149418; 0.479837337099708772|];
[|-1.08194373033462554; 0.908093424482677825|];
[|-0.972077812548367781; 1.30488663625732659|];
[|-0.82505430457777329; 1.65601407508140652|];
[|-0.646240433036617379; 1.94394948589238603|];
[|-0.442829622786570176; 2.14981996365949968|];
[|-0.223801781785368276; 2.2571237113826963|];
[|0.000528489775196344658; 2.25607234986715977|];
[|0.219348908747622; 2.14671521943166699|];
[|0.422366789921064; 1.93929511007784194|];
[|0.600710760394784216; 1.65167055463977475|];
[|0.74748362573778071; 1.30514678495061487|];
[|0.857902204805326196; 0.920583984378085307|];
[|0.929122154365651; 0.516052437651637552|];
[|0.959921681840344321; 0.106299134247129612|];
[|0.950405114828223; -0.296423506018091476|];
[|0.901827225513500808; -0.68086468247841|];
[|0.816570279638364926; -1.03528621433651136|];
[|0.698244496501420175; -1.34654412027003123|];
[|0.551832059877745085; -1.60018721507722805|];
[|0.383760349982204596; -1.78179552755520376|];
[|0.201787149120740694; -1.87936224564236398|];
[|0.0146267505409670706; -1.88598419251135474|];
[|-0.168661078918962176; -1.80178276010032268|];
[|-0.339395470111264641; -1.63417945592763658|];
[|-0.489935641633520325; -1.39639794156977781|];
[|-0.614132507559556617; -1.10489127599143466|];
[|-0.7075165529098707; -0.776757359358784427|];
[|-0.767265518869675645; -0.427970848742863064|];
[|-0.79205159527398461; -0.0727208064874359|];
[|-0.781871429406644447; 0.276331513203226298|];
[|-0.737929788874430503; 0.607306508710532156|];
[|-0.662601233432069781; 0.908353861544634|];
[|-0.559448723938056669; 1.16731564334651239|];
[|-0.433242417624848775; 1.37209070227011587|];
[|-0.289903982741078; ...|]; ...|]
In [158]:
let vp = A.init ~w:800. ~h:360. ["jupyter"] in
A.Axes.box vp ;
A.Viewport.xlabel vp "Time t";
A.Viewport.title vp "Solution to the pendulum ODE with Runge-Kutta 4 (101 points)";
A.set_color vp A.Color.red ;
A.Viewport.text vp 8.5 3. "o theta(t)" ;
A.Array.xy ~style:(`Linesmarkers "o") vp t (column sol 0);
A.set_color vp A.Color.blue ;
A.Viewport.text vp 8.5 2.5 "+ omega(t)" ;
A.Array.xy vp ~style:(`Linesmarkers "+") t (column sol 1);
A.close vp;;
Out[158]:
- : unit = ()
Comparisons¶
In [120]:
let methods = [|rungekutta1; rungekutta2; rungekutta4|] ;;
let names = [|"Runge-Kutta 1"; "Runge-Kutta 2"; "Runge-Kutta 4"|] ;;
let markers = [|"o"; "s"; ">"|] ;;
let colors = [|A.Color.red; A.Color.blue; A.Color.green|] ;;
Out[120]:
val methods : ((float array -> float -> float array) -> float array -> float array -> float array array) array = [|<fun>; <fun>; <fun>|]
Out[120]:
val names : string array = [|"Runge-Kutta 1"; "Runge-Kutta 2"; "Runge-Kutta 4"|]
Out[120]:
val markers : string array = [|"o"; "s"; ">"|]
Out[120]:
val colors : A.Color.t array = [|<abstr>; <abstr>; <abstr>|]
In [159]:
let test_1 ?(n=101) () =
let t = linspace 0. 10. n in
let vp = A.init ~w:800. ~h:360. ["jupyter"] in
A.Axes.box vp ;
A.Viewport.xlabel vp "Time t";
A.Viewport.title vp (Format.sprintf "Solution to the pendulum ODE with different methods (%i points)" n);
for i = 0 to (Array.length methods) - 1 do
let sol = methods.(i) f_pend y0 t in
A.set_color vp colors.(i);
A.Viewport.text vp 6.5 (2. -. 2.*.(float_of_int i)) (Format.sprintf "%s %s" markers.(i) names.(i)) ;
A.Array.xy ~style:(`Linesmarkers markers.(i)) vp t (column sol 0);
done;
A.close vp;
;;
Out[159]:
val test_1 : ?n:int -> unit -> unit = <fun>
In [160]:
test_1 ~n:10 ();;
Out[160]:
- : unit = ()
In [161]:
test_1 ~n:30 ();;
Out[161]:
- : unit = ()
In [162]:
test_1 ~n:100 ();;
Out[162]:
- : unit = ()
Conclusion¶
That's it for today, folks! See my other notebooks, available on GitHub.