Legend:
Page
Library
Module
Module type
Parameter
Class
Class type
Source
Source file common.ml
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130(*
* OWL - OCaml Scientific and Engineering Computing
* OWL-ODE - Ordinary Differential Equation Solvers
*
* Copyright (c) 2019 Ta-Chu Kao <tck29@cam.ac.uk>
* Copyright (c) 2019 Marcello Seri <m.seri@rug.nl>
*)(* TODO: update implementations of multiple order RK on the line of
* symplectic.ml *)(* TODO: find a better place to place this module *)moduleMake(M:Owl_types_ndarray_algodiff.Sigwithtypeelt=float)=structletstepst0t1dt=(t1-.t0)/.dt|>floor|>int_of_float|>succtypestate_type=|Row|Col|Matrixletget_state_ty0=letdims=M.shapey0inletdim1,dim2=dims.(0),dims.(1)inifdim1=1thenRow,dim2elseifdim2=1thenCol,dim1elseMatrix,dim1*dim2letintegrate~step~tspan:(t0,t1)~dty0=letstate_t,n=get_state_ty0inletn_steps=stepst0t1dtinletys=matchstate_twith|Row->M.empty[|n_steps;n|]|Col->M.empty[|n;n_steps|]|Matrix->M.empty[|n_steps;n|]inletts=ref[]inlett=reft0inlety=refy0infori=0topredn_stepsdoifi>0then(lety',t'=step!y!tiny:=y';t:=t');(matchstate_twith|Row->M.set_slice[[i];[]]ys!y|Col->M.set_slice[[];[i]]ys!y|Matrix->M.set_slice[[i];[]]ysM.(reshape!y[|1;-1|]));ts:=!t::!tsdone;letts=[|!ts|>List.rev|>Array.of_list|]|>M.of_arraysinmatchstate_twith|Row|Matrix->M.(transposets),ys|Col->ts,ysletsymplectic_integrate~step~tspan:(t0,t1)~dt(x0,p0)=ifM.shapex0<>M.shapep0thenraiseOwl_exception.(DIFFERENT_SHAPE(M.shapex0,M.shapep0));letstate_t,n=get_state_tx0inletn_steps=stepst0t1dtinletxs,ps=matchstate_twith|Row->M.empty[|n_steps;n|],M.empty[|n_steps;n|]|Col->M.empty[|n;n_steps|],M.empty[|n;n_steps|]|Matrix->M.empty[|n_steps;n|],M.empty[|n_steps;n|]inletts=ref[]inlett=reft0inletx=refx0inletp=refp0infori=0topredn_stepsdoifi>0then(let(x',p'),t'=step(!x,!p)!tinx:=x';p:=p';t:=t');(matchstate_twith|Row->M.set_slice[[i];[]]xs!x;M.set_slice[[i];[]]ps!p|Col->M.set_slice[[];[i]]xs!x;M.set_slice[[];[i]]ps!p|Matrix->M.set_slice[[i];[]]xsM.(reshape!x[|1;-1|]);M.set_slice[[i];[]]psM.(reshape!p[|1;-1|]));ts:=!t::!tsdone;letts=[|!ts|>List.rev|>Array.of_list|]|>M.of_arraysinmatchstate_twith|Row|Matrix->M.transposets,xs,ps|Col->ts,xs,psletadaptive_integrate~step~tspan:(t0,t1)~dtmaxy0=letstate_t,_=get_state_ty0inletdt=dtmax/.4.0inletrecgo(ts,ys)(t0:float)y0dt=ift0>=t1thents,yselse(letdt=mindt(t1-.t0)inift0+.dt<=t0thenfailwith"Singular ODE";lety,t,dt,err_ok=step~dty0t0iniferr_okthen(* Update solution if error is OK *)go(t::ts,y::ys)tydtelsego(ts,ys)t0y0dt)inletts,ys=go([t0],[y0])t0y0dtinletts=[|ts|>List.rev|>Array.of_list|]|>M.of_arraysinletys=matchstate_twith|Row->ys|>List.rev|>Array.of_list|>M.of_rows(* TODO: add of_cols to types_ndarray_basic *)|Col->ys|>List.rev|>Array.of_list|>M.of_cols|Matrix->ys|>List.rev|>Array.of_list|>Array.map(funy->M.reshapey[|1;-1|])|>M.of_rowsinmatchstate_twith|Row|Matrix->M.transposets,ys|Col->ts,ysend