So I am a beginner with stochastic differential equations and came across Mathematica's capabilities for solving them. I am solving the stochastic Mathieu equation with a harmonic forcing term that has a stochastic fluctuation imposed.

The linear set of ODEs to solve this are (ready for Ito form):

$ dx(t)=y(t) dt,\\ dy(t)=\lambda (t) dt \cos (\nu t) \sin (x(t))-\epsilon \,\textbf{dw(t)} \sin (x(t))$

Where $\lambda$ is time-varying as in the code below. $\nu$ possesses some physical characteristic I cannot place at the moment but takes real values. $\epsilon$ is a small number $<<1$. $w$ is the Wiener process distribution with zero drift and volatility (Ariaratnam et al.) (google drive link to paper. Only introduction is relevant for now.).

tmax = 100;
\[Lambda][t_] := 
 Piecewise[{{.95, t <= tmax/4}, {1.1, tmax/4 < t <= tmax/2}, {1.25, 
    t > tmax/
      2}}]; \[Nu] = 105.; \[Epsilon] = 0.001; x0 = 0.1; y0 = 0.1; \
drift = 0; volatility = 1.;

sme = ItoProcess[
    \[DifferentialD]x[t] == y[t] \[DifferentialD]t,
       t] == \[Lambda][t] Cos[\[Nu] t] Sin[
        x[t]] \[DifferentialD]t - \[Epsilon] Sin[
        x[t]] \[DifferentialD]w[t]
    }, x[t],
   {{x, y}, {x0, y0}},
   {t, 0}, w \[Distributed] WienerProcess[drift, volatility]
path = RandomFunction[sme, {0., tmax, 1}, 1]

My main issue is with the path = RandomFunction[sme, {0., tmax, 1}, 1] line. Here, it is my understanding that the {0., tmax, 1} is a time-stepping from 0 to tmax in steps of 1.

With a "time-step" of 1, I get the following dynamics which are reasonable, I would think:

ListLinePlot[path, PlotRange -> All, GridLines -> None, 
 ImageSize -> 250 {1, 1}, Frame -> True, FrameStyle -> Black, 
 FrameLabel -> {"time", "position"}, 
 PlotStyle -> {Thickness[0.005], Black}, 
 BaseStyle -> {FontSize -> 14, Bold}, Axes -> False]

enter image description here

However, if I change the time-step to, say, 0.1, I see an instability (or is it?) with some sort of monotonically increasing position.

enter image description here

Am I going wrong somewhere? Any input or assistance, I am grateful for.

Edit: The position plot is for x[t].


I think your dt=1 solution is the wrong one. The Cos[ν t] term is very rapidly varying with your parameter value of ν=105.:

Plot[Cos[ν t], {t, 0, tmax}, PlotPoints -> 200]

Mathematica graphics

The apparent periodic-like solution you get with dt=1. might be due to an aliasing effect of sampling the high frequency cosine with period 1. To get the right solution, I'd pick a time step much less than the period of the forcing.

Going out on a limb here, I'd guess that we could replace that rapidly varying Cos[ν t] term by its average of zero.

sme = ItoProcess[{\[DifferentialD]x[t] == y[t] \[DifferentialD]t,
  \[DifferentialD]y[t] == -\[Epsilon] Sin[x[t]] \[DifferentialD]w[t]},
  x[t], {{x, y}, {x0, y0}}, {t, 0}, 
  w \[Distributed] WienerProcess[drift, volatility]];
path = RandomFunction[sme, {0., tmax, 0.01}, 1]

Mathematica graphics

Further speculation: the time average of -ϵ Sin[x[t]] \[DifferentialD]w[t] is zero, so the net effect of noise is zero. We might get some insights from the resulting ODE:

sol = NDSolve[{x'[t] == y[t], y'[t] == 0, x[0] == x0, y[0] == y0}, {x, y}, {t, 0, tmax}];
Plot[x[t] /. sol, {t, 0, tmax}]

Mathematica graphics

I think for more interesting dynamics, you'll want to lower ν.

  • $\begingroup$ Good points. I'll try and implement these in the next day or so and see what cranks out. $\endgroup$ – drN Jun 25 '18 at 15:15
  • $\begingroup$ Is there a way to choose these parameters $\lambda$, $\nu$? Is there some theory I should study first, that will help me make these decisions? $\endgroup$ – drN Jun 25 '18 at 18:08
  • $\begingroup$ @drN Sorry, this is totally outside my domain of expertise. Good luck! $\endgroup$ – Chris K Jun 25 '18 at 19:59

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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