1

I have following system of differential equations

x'(t) = (10+σ(t))(y(t)-x(t))
y'(t) = ρx(t)-y(t)-x(t)z(t)
z'(t) = x(t)y(t)-βz(t)

where ρ and β are constants and σ is Ornstein-Uhlenbeck process. In fact, the system is Lorenz system embedded in stochastic environment.

Solution of σ is

σt=σ0e-αt+γe-αt∫0teαsdWs

and I think I can solve it in Python like below

t_max = 1
n = 1000
t, dt = np.linspace(0, t_max, n, endpoint=False, retstep=True)
dW = np.sqrt(dt) * np.random.randn(n)
alpha = 1
gamma = 1
sigma = sigma_0 * np.exp(-alpha * t) + gamma * np.exp(-alpha * t) * np.cumsum(np.exp(alpha * t) * dW)

Now given sigma as above, if it is correct, how can I get the solution for the system?

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
Midess
  • 11
  • 3
  • Now make the function table (t,sigma) into an interpolation function and use that in the ODE function. – Lutz Lehmann Feb 08 '23 at 19:28
  • do you want a symbolic or numeric solution? – anon01 Feb 08 '23 at 19:32
  • @LutzLehmann Should `scipy.interpolate.interp1d` be enough? Or is there more "accurate" method of interpolation? – Midess Feb 08 '23 at 21:36
  • @anon01 I want numeric solution so I can experiment with effects of α and γ to the whole system. – Midess Feb 08 '23 at 21:37
  • Yes, that is the advanced option with several methods to select from in the background. Piecewise linear should be sufficient, over-smoothing with piecewise cubics could have strange effects. Or take a dedicated SDE solver that solves for all 4 dynamic components simultaneously. – Lutz Lehmann Feb 08 '23 at 22:08

0 Answers0