[R-sig-dyn-mod] lorenz equation?
Thomas Petzoldt
thomas.petzoldt at tu-dresden.de
Tue Mar 27 22:08:57 CEST 2018
Hi
the order of your states was mixed up. They must be identical in init
and in the return list of main. This is a feature of deSolve (in the
interest of efficiency) and a FAQ ;-)
I admit that it is not trivial to find this within the the amount of
available docs.
Here a version that worked
(word wrapped, because of long variable names):
library("simecol")
library("scatterplot3d")
lorenz.model <- new("odeModel",
main = function(times, y, parms) {
with(as.list(c(parms, y)), {
horizTempDt <- convectiveFlowDelta *
(rayleigh - verticalTempDelta) - horizTempDelta
verticalTempDt <- convectiveFlowDelta * horizTempDelta -
height * verticalTempDelta
convectiveFlowDt <- prandtl * (horizTempDelta -
convectiveFlowDelta)
list(c(convectiveFlowDt, horizTempDt, verticalTempDt))
})
},
times = seq(0, 50, .0078125),
parms = c(prandtl=10, rayleigh=28, height=8/3),
init = c(convectiveFlowDelta=0, horizTempDelta=1,
verticalTempDelta=0),
solver = "rk4"
)
lorenz.sim <- sim(lorenz.model)
plot(lorenz.sim)
scatterplot3d(out(lorenz.sim)[,2:4], type="l")
Two notes:
1) I recommend to use "lsoda" instead of "rk4". It is faster and more
precise.
2) The Lorenz equation is THE main example of package deSolve, the
differential equation solver package that simecol depends on.
see:
?deSolve
example(deSolve)
Hope it helps,
Thomas
Am 27.03.2018 um 18:31 schrieb RC Phelps:
> Hi,
>
> I can't duplicate the vensim expression of the lorenz equation with
> simecol. I wonder could there be a difference in the RK4 integration
> implementation?
>
> If anyone is interested in lorenz, maybe you can see my error on:
>
> https://github.com/cordphelps/sim/tree/master/porting-lorenz
> <https://github.com/cordphelps/sim/tree/master/porting-lorenz>
>
> Thank you for any suggestions.
>
> Cord
>
More information about the R-sig-dynamic-models
mailing list