5
$\begingroup$

I have the following system of differential equations:

v[t_] := RandomVariate[NormalDistribution[]];
sol = NDSolve[{x1'[t] == -0.33 x1[t] - 1.13 x2[t] - 1.84 x3[t] - 
      1.22 x4[t] + v[t], 
    x2'[t] == 
     0.15 x1[t] - 0.57 x2[t] + 0.29 x3[t] + 0.28 x4[t] + v[t], 
    x3'[t] == 
     0.24 x1[t] + 0.34 x2[t] - 0.48 x3[t] + 0.38 x4[t] + v[t], 
    x4'[t] == 
     0.17 x1[t] + 0.18 x2[t] + 0.32 x3[t] - 0.56 x4[t] + v[t], 
    x1[0] == 0, x2[0] == 0, x3[0] == 0, x4[0] == 0}, {x1, x2, x3, 
    x4}, {t, 0, 60}];
graphAll = 
 Plot[Evaluate[{x1[t], x2[t], x3[t], x4[t]} /. sol], {t, 0, 60}, 
  PlotRange -> All]

I want to drive the system with an AR1-Process. So far I was only able to specify some random variate from a normal distribution, but it does not seem to work properly.

Any ideas how to do this?

$\endgroup$
6
$\begingroup$

You need to use RandomFunction[ItoProcess[]] instead of NDSolve. The syntax is a frustratingly a little different than NDSolve. Driven by a WienerProcess as in your example:

sol = RandomFunction[ItoProcess[{
  \[DifferentialD]x1[t] == \[DifferentialD]v[t] + 
    (-0.33 x1[t] - 1.13 x2[t] - 1.84 x3[t] - 1.22 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x2[t] == \[DifferentialD]v[t] + 
    (0.15 x1[t] - 0.57 x2[t] + 0.29 x3[t] +0.28 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x3[t] == \[DifferentialD]v[t] +
    (0.24 x1[t] + 0.34 x2[t] - 0.48 x3[t] + 0.38 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x4[t] == \[DifferentialD]v[t] +
    (0.17 x1[t] + 0.18 x2[t] + 0.32 x3[t] - 0.56 x4[t]) \[DifferentialD]t},
  {x1[t], x2[t], x3[t], x4[t]}, {{x1, x2, x3, x4}, {0, 0, 0, 0}}, t,
  v \[Distributed] WienerProcess[0, 1]], {0, 100, 0.01}];
ListLinePlot[sol, PlotRange -> All]

Mathematica graphics

For the AR-1 part, I think you need to add a first-order decay equation for v[t] driven by a WienerProcess like:

ϕ = 1;
σ = 1;
sol = RandomFunction[ItoProcess[{
  \[DifferentialD]x1[t] == (v[t] - 0.33 x1[t] - 1.13 x2[t] - 1.84 x3[t] - 1.22 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x2[t] == (v[t] + 0.15 x1[t] - 0.57 x2[t] + 0.29 x3[t] + 0.28 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x3[t] == (v[t] + 0.24 x1[t] + 0.34 x2[t] - 0.48 x3[t] + 0.38 x4[t]) \[DifferentialD]t,
  \[DifferentialD]x4[t] == (v[t] + 0.17 x1[t] + 0.18 x2[t] + 0.32 x3[t] - 0.56 x4[t]) \[DifferentialD]t,
  \[DifferentialD]v[t] == -ϕ v[t] \[DifferentialD]t + σ \[DifferentialD]W[t]},
  {x1[t], x2[t], x3[t], x4[t], v[t]}, {{x1, x2, x3, x4, v}, {0, 0, 0, 0, 0}}, t,
  W \[Distributed] WienerProcess[0, 1]], {0, 100, 0.01}];
ListLinePlot[sol]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Works very well, thank you Chris! $\endgroup$ – holistic Oct 24 '18 at 14:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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