2

I have the following ODE:

ODE

where p is a probability, and y is a random variable with pdf:

see link

epsilon is a small cut-off value (typically 0.0001).

I am looking to numerically solve this system in Python, for t = 0 to about 500.

Is there a way I can implement this using numpy/scipy?

tel
  • 13,005
  • 2
  • 44
  • 62
Mr. Bump
  • 23
  • 2
  • 1
    A simple google search will lead you [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html). – MaJoR Nov 24 '18 at 17:57
  • 4
    @MaJoR: [`scipy.integrate.odeint`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) does not solve [stochastic differential equations](https://en.wikipedia.org/wiki/Stochastic_differential_equation). – Warren Weckesser Nov 24 '18 at 18:17
  • https://stackoverflow.com/questions/14594290/is-there-a-python-module-to-solve-integrate-a-system-of-stochastic-differential – kevinkayaks Nov 24 '18 at 22:59
  • @WarrenWeckesser, excuse me for that, but when I commented that bit, OP had not mentioned what type of DEs he wanted to solve. He simply blurted out ODEs. And for a generic question, I left a generic comment. – MaJoR Nov 25 '18 at 06:05

1 Answers1

1

There's an annoying lack of quality Python packages for SDE integration. You have kind of an unusual setup, in that your random variable has an explicit dependence on your integrand (at least it appears that y depends on p). Because of this, it'll probably be hard to find a preexisting implementation that meets your needs.

Fortunately, the simplest method for SDE integration, Euler-Maruyama, is very easy to implement, as in the eulmar function below:

from matplotlib import pyplot as plt
import numpy as np

def eulmar(func, randfunc, x0, tinit, tfinal, dt):
    times = np.arange(tinit, tfinal + dt, dt)
    x = np.zeros(times.size, dtype=float)
    x[0] = x0

    for i,t in enumerate(times[1:]):
        x[i+1] = x[i] + func(x[i], t) + randfunc(x[i], t)

    return times, x

You could then use eulmar to integrate your SDE like so:

def func(x, t):
    return 1 - 2*x

def randfunc(x, t):
    return np.random.rand()

times,x = eulmar(func, randfunc, 0, 0, 500, 5)
plt.plot(times, x)

enter image description here

You'll have to supply your own randfunc however. Like above, it should be a function that takes x and t as arguments and returns a single sample from your random variable y. If you're having trouble coming up with a way to generate samples of y, since you know the PDF you can always use rejection sampling (though it does tend to be fairly inefficient).

Notes

  • This is not a particularly efficient implementation of Euler-Maruyama. For example, the random samples are usually generated all at once (eg np.random.rand(500)). However, you can't pre-generate your random samples, since y depends on p.
tel
  • 13,005
  • 2
  • 44
  • 62