# Generating Lorenz

17 views

### idontgetoutmuch

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)

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

Jun 15, 2018, 5:52:40 AM6/15/18

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.

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/

### idontgetoutmuch

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

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

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`.