more ode example problems

[pnkraemer]

Would be great to have more ODE example problems. Currently we only have logistic, lotkavolterra, fitzhughnagumo.

Further, I'd like to also have

• three-body (I will post code below)
• seir (https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model)
• rigid-body (http://www.ece.northwestern.edu/local-apps/matlabhelp/techdoc/math_anal/diffeq8.html)
• van-der-pol (also http://www.ece.northwestern.edu/local-apps/matlabhelp/techdoc/math_anal/diffeq8.html); the default here should probably a non-stiff version...

and perhaps (if necessary, in some simplified version)

• hodgkin-huxley (https://gist.github.com/giuseppebonaccorso/60ce3eb3a829b94abf64ab2b7a56aaef) because I catch myself doing every example with lotka-volterra. That is super boring, but probably due to its availability.

[pnkraemer]

Three body code (variables may require naming that is consistent with the rest of the probnum source; the snippet is shamelessly stolen)

def threebody(t,Y):
    # defining the ODE:
    # assume Y = [y1,y2,y1',y2']
    mu = 0.012277471 #  (standardised moon mass)
    mp = 1 - mu    D1 = ((Y[0] + mu)**2 + Y[1]**2)**(3/2)
    D2 = ((Y[0] - mp)**2 + Y[1]**2)**(3/2)    y1p = Y[0] + 2 * Y[3] - mp * (Y[0] + mu) / D1 - mu * (Y[0] - mp) / D2
    y2p = Y[1] - 2 * Y[2] - mp * Y[1] / D1 - mu * Y[1] / D2    return np.array([Y[2],Y[3],y1p,y2p])# defining the initial value
iv = np.array([ 0.994, 0, 0, -2.00158510637908252240537862224 ])
initrv = rvs.Dirac(iv)# defining the time span
Tmax = 17.0652165601579625588917206249
timespan = (0, Tmax)ivp = IVP(timespan, initrv, threebody)

[pnkraemer]

@schmidtjonathan, while you're working on this, what do you think about moving the examples into a separate file, called examples.py or something along these lines? It might slim down ivp.py a little bit...

[schmidtjonathan]

@pnkraemer I wanted to suggest the same thing, I think it's good to separate functionality from examples. So, yes, I think that's a good idea, and I can include this in the pull request, if you want

[pnkraemer]

Please do it then, thank you lots!

[schmidtjonathan]

@pnkraemer When writing tests for the IVPs I noticed that you can create IVPs with invalid parameters. What do you think about checks (we would have to talk about, how to do that), whether or not the provided params are valid? I am not necessarily talking about value ranges or anything like that, but at least maybe we could check the number of provided parameters in the params tuple?

EDIT: By the way, same thing applies for the initial conditions RV.

[pnkraemer]

Hmmm good point. Though I am in favour of just allowing this at the moment and postponing the solution. I am not sure whether we are at the stage yet, where it is not a waste of our time to do adversarial testing. What do you think?

Initial condition checks, on the other hand, I find more urgent, because if this fails it will be a cryptic error down the line (I'd guess in ivp2filter) which may be hard to interpret. Would you be so kind and open a separate issue for this? Let's solve this in a different PR.

[pnkraemer]

Another thing: if you have the time, can you quickly brush over the docs for the ODE examples and check that the docstring says which dimensionality the initial random variable should have? E.g. write something along the lines of (for a 2d problem such as Lotka-Volterra):

initrv: RandomVariable
    2-dimensional random variable describing the initial state of the ODE. 
    Usually, `pnrv.Constant` supported on a 2-dimensional array. 
    Alternatively, `pnrv.Normal` with a 2-dimensional mean and 2-dimensional array

[schmidtjonathan]

Sure, I'll do that today.