3
$\begingroup$

I'm attempting to add noise to a set of ODE's with two state variables. $$\frac{dx}{dt} = 10 -(x-1)\left(1+\frac{exp\left(\frac{x-1}{5y}\right)}{50y}\right)$$ $$\frac{dy}{dt} = 2(1-y) -y\cdot exp\left(\frac{x-1}{5y}\right) $$

The numerical solution for the ODE is below

enter image description here

As u can see around $t=0.4$ there is a sharp drop in $y$ (orange curve) which causes the term $exp(\frac{x-1}{5y})$ to explode, nevertheless NDSolve handles the integration well.

Before adding the noise, I'm attempting to solve the equations using RandomFunction without any addition of noise, the result is below

enter code here

As u can see, this solver does not integrate the solution properly, and the result diverges.

If I reduce $dt$ to be very small it sometimes works out the problem, however, my original problem is far more complex than here, and there is no way that solving with $dt=0.0000001$ or lower can work out (tried that).

I've attempting both ItoProcess, StratanovichProcess, all the integration methods and working precision possible of RandomFunction.

Can anyone please advise on how to handle such a situation with mathematica?

I would like to keep $dt$ reasonably small, or maybe varying such that it is close to $n=0$ it will become small.

I don't see any reason why such a calculation will produce this result, without even having an additive noise.

here is the sample code below

f[x_, y_] := 10 - (x - 1) (1 + Exp[(x - 1)/(5 y)]/(50 y  )); 
g[x_, y_] := 2 (1 - y) - y Exp[(x - 1)/(5 y)];     
sol = NDSolveValue[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]], 
    x[0] == 2.1, y[0] == 0.4}, {x[t], y[t]}, {t, 0, 0.5}];
sol2 = RandomFunction[
   ItoProcess[{\[DifferentialD]x[t] == 
      f[x[t], y[t]] \[DifferentialD]t, \[DifferentialD]y[t] == 
      g[x[t], y[t]] \[DifferentialD]t}, {x[t], 
     y[t]}, {{x, y}, {2.1, 0.4}}, t, 
    w \[Distributed] WienerProcess[0, 1]], {0, 0.5, 0.00001}];
Plot[{sol[[1]], sol[[2]]}, {t, 0, 0.5}]
ListLinePlot[sol2, PlotRange -> All]
Improve this question
$\endgroup$
2
  • $\begingroup$ What solution do you want to get? Numerical solution of ODE and random process are different objects. You can't just mix them up. $\endgroup$
    – Alex Trounev
    Commented Nov 19, 2018 at 15:42
  • $\begingroup$ @AlexTrounev, at this point I'm trying to solve using mathematica's tools the stochastic differential equations above with a noise term that is identically zero. This is in order to see that the integration works fine before fluctuations are added. If you have any suggestions of solving these equations in another way besides RandomFunction, please post your solution and demonstrate that the integration holds before and after the addition of noise. $\endgroup$
    – jarhead
    Commented Nov 19, 2018 at 15:47

1 Answer 1

Reset to default
3
$\begingroup$

Here's a crazy idea: maybe it's easier to add noise to NDSolve's adaptive step size algorithms than to deal with RandomFunction[ItoProcess[]] 's fixed step size. You could use WhenEvent to perturb the solution by a normally distributed amount every δ time steps.

δ = 10^-3;
{σx, σy} = {1, 10^-8};
sol = NDSolveValue[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]], 
  WhenEvent[Mod[t, δ] == 0, {
   x[t] -> RandomVariate[NormalDistribution[x[t], Sqrt[δ] σx]],
   y[t] -> RandomVariate[NormalDistribution[y[t], Sqrt[δ] σy]]}], 
  x[0] == 2.1, y[0] == 0.4}, {x[t], y[t]}, {t, 0, 0.5}, MaxSteps -> Infinity];

Plot[{sol[[1]], sol[[2]]}, {t, 0, 0.5}]

Mathematica graphics

Note I scaled the standard deviation of the noise by Sqrt[δ], which I think is correct but someone who knows more about numerically solving stochastic differential equations should check it out. In fact, I am not sure this approach is even legit at all.

Adding larger noise to y[t] might be tricky when y[t] gets close to zero.

Here's an example with a larger noise time step δ = 10^-2 to give a better idea what's happening:

Mathematica graphics

Improve this answer
$\endgroup$
5
  • $\begingroup$ This looks great, but I'm a bit worried with the time step. 1) you fix δ, but how do u know that this is the time step used in the solution, also u use MaxSteps->infinity, can u explain please if it contradicts. 2) in the sde using simple euler-mayurana we have dx=f(x,y)*dt+σxsqrt(dt), it seems here that what that is calculated when WhenEvent applies is dx=d(x,y)*dt+σxsqrt(dt)*dt, can you please clarify this issue. $\endgroup$
    – jarhead
    Commented Nov 19, 2018 at 18:42
  • $\begingroup$ There are two distinct time steps here: the fixed δ where noise is added and the automatic time steps chosen by NDSolve between noise events. Check out the plot I added with δ = 10^-2 for a better idea what it's doing. I don't see the extra dt in the noise term you mention in the second question. I think the Euler-Maruyama is $dx=f(x,y)*dt+N(0,σx \sqrt{dt})$ where $N()$ is a normal variate. $\endgroup$
    – Chris K
    Commented Nov 19, 2018 at 18:53
  • $\begingroup$ To simulate the ito process with a noise for x (that's only what I'm interested in) we need to solve: dx=f(x,y)dt+dW, dy=g(x,y)dt, where dW=z(t)*sqrt(dt), where z(t) is an uncorrelated white noise. The sqrt(dt) is normalization otherwise the solutions depend on the time step. to have your solution work (in which I wish it will work), the above must be satisfied. so first the time steps must be equal, such that the noise is added in every steps. maybe you can also provide a simple known problem to compare with your method and RandomFunction to show that both ways are equivalent? $\endgroup$
    – jarhead
    Commented Nov 19, 2018 at 19:55
  • $\begingroup$ It seems like a plausible approximation to me, but I'm not an expert in SDEs. If you want to test it on a known problem, let us know if it works out. $\endgroup$
    – Chris K
    Commented Nov 19, 2018 at 20:10
  • $\begingroup$ Since it really helped with my problem I'll accept the answer, but before please clarify that in the body of your answer that this is an approximation and explain exactly what you do with the time steps, so other people will not get confused. good job. $\endgroup$
    – jarhead
    Commented Nov 19, 2018 at 20:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.