Skip to content

[R-sig-dyn-mod] lorenz equation?

2 messages · RC Phelps, Thomas Petzoldt

Original
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
Thomas Petzoldt
# 
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: