Generating Lorenz

17 views
Skip to first unread message

idontgetoutmuch

unread,
Jun 15, 2018, 5:44:11 AM6/15/18
to LibBi Users
I am trying to generate values for a Lorenz system. Here's my system in Haskell

alpha = 16.0
rho = 45.92
beta = 4

dzdt :: Double -> [Double] -> [Double]
dzdt _t [x, y, z] = [ alpha * (y - x)
                      , x * (rho - z) - y
                      , x * y - beta * z
                      ]

m = odeSolve dzdt [1.0, 1.0, 1.0] (LA.vector $ take 1000 [1.0e-4, 2.0e-4..])

And I get a solution which looks reasonable

*Main> m
(1000><3)
 [                1.0,                1.0,                1.0
 , 1.0000035116198551, 1.0043918006561992, 0.9997002796734803
 , 1.0000140385819727, 1.0087832236309844, 0.9994011189915194
 ,  1.000031569092141, 1.0131743003919085, 0.9991025184087792
 , 1.0000560914252676, 1.0175650623511352,   0.99880447839092
 , 1.0000875939251799, 1.0219555408657492, 0.9985069994145762
 ,  1.000126065004428, 1.0263457672380667, 0.9982100819673314
 ,  1.000171493144087, 1.0307357727159434, 0.9979137265476951
 ,  1.000223866893561, 1.0351255884930834, 0.9976179336650787
...

If I try this in libbi

model Lorenz {
  const rho = 45.92
  const beta = 4.0
  const alpha = 16.0

  state X, Y, Z

  sub initial {
    X ~ log_normal(log(1.0), 0.00002)
    Y ~ log_normal(log(1.0), 0.00002)
    Z ~ log_normal(log(1.0), 0.00002)
  }

  sub transition(delta = 0.0001) {
    ode {
      dX/dt = alpha * (Y - X)
      dY/dt = X * (rho - Z) - Y
      dZ/dt = X * Y - beta * Z
    }
  }
}


then I get nothing like the results I expect. Even if X, Y and Z are not exactly 1.0 initially, I think at least at time t = 0.0001 I should get values close to 1.

> T <- 1000
> init_parameters <- list(X = 1, Y = 1, Z = 1)

> synthetic_dataset <- bi_generate_dataset(end_time=T, model=Lorenz,
                                           init=init_parameters)

> synthetic_data <- bi_read(synthetic_dataset)

> head(synthetic_data$X)
  time     value
1    0   1.00000
2    1 -12.65339
3    2 -13.72596
4    3 -13.57361
5    4 -11.26844
6    5 -15.36319


Can anyone tell me what I am doing wrong? Perhaps I need to state somewhere that I want results at t=0.0001, 0.0002, ... but I couldn't see how to do this.

Lawrence Murray

unread,
Jun 15, 2018, 5:52:40 AM6/15/18
to libbi...@googlegroups.com

The end_time in LibBi is a real time, not a number of time steps. So it looks like your Haskell code simulates from time 0 to 0.001 with 1000 outputs, while your LibBi code simulates from time 0 to 1000 with 1000 outputs, which is bound to look different. To match your Haskell results in LibBi, you will want to set end_time to 0.001, and set the number of outputs to 1000 (I'm not sure what the exact argument for this is in RBi, with the command line interface it is --noutputs).

--
You received this message because you are subscribed to the Google Groups "LibBi Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to libbi-users...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.









När du har kontakt med oss på Uppsala universitet med e-post så innebär det att vi behandlar dina personuppgifter. För att läsa mer om hur vi gör det kan du läsa här: http://www.uu.se/om-uu/dataskydd-personuppgifter/

E-mailing Uppsala University means that we will process your personal data. For more information on how this is performed, please read here: http://www.uu.se/om-uu/dataskydd-personuppgifter/

idontgetoutmuch

unread,
Jun 15, 2018, 5:56:49 AM6/15/18
to LibBi Users
In fact I just ran my Haskell code to 10,000 steps i.e. t=1.0 and got -12.652290724996961 which pretty much agrees with libbi so now I just need to convince libbi to return me solutions at 0.0001, 0.0002, ...

idontgetoutmuch

unread,
Jun 15, 2018, 6:04:45 AM6/15/18
to LibBi Users
Thanks for this - I have stared hard but can find nothing in https://cran.r-project.org/web/packages/rbi/rbi.pdf that allows me to set that

Sebastian Funk

unread,
Jun 15, 2018, 6:09:18 AM6/15/18
to idontgetoutmuch, LibBi Users
Any libbi option can just be passed as an argument to `bi_generate_dataset` (or sample/filter/predict). So you'd pass `noutputs=1000`.
>> email to libbi-users...@googlegroups.com <javascript:>.

idontgetoutmuch

unread,
Jun 15, 2018, 6:37:25 AM6/15/18
to LibBi Users
Great - thanks
Lorenz.png
Reply all
Reply to author
Forward
0 new messages