# Adding noise to nonlinear ODE system

I have a problem with adding a noise to 3d order ODE nonlinear system. I solved the system numericaly and got the periodic solution and now I need to add random pertubations. Due to the lack of experience of working with SDE I couldn't get the noisy plot. Here's the code:

Clear["Global*"]
(******It's all just params******************************)
δ = 0.5100;
κ = 0.7100;
α = 0.5100;
tl = 270;
too = 280;
Ocean = 500;
Land = 25;
γ = 2.256100;
μ = 0.335100;
k = 5100;
M = 0.9999100 10^4;
vo = 0.746326100 10^(17 - 8/κ);
a = SetPrecision[(0.1666 k)/((too - tl) γ), 100];
b = SetPrecision[(δ a (too - tl))/(vo^κ tl), 100];
p = SetPrecision[a (too - tl), 100];
c = SetPrecision[1/too (k α - (γ p)), 100];
α = SetPrecision[(γ p + c too)/k, 10];
β =SetPrecision[-((γ p - c tl)/((1 - vo^κ/M) k)), 100];
\[Sigma] = SetPrecision[(δ a β)/(M b), 100];
(************Computing stability points****************)
bifP = 119;
f[x_, y_, z_] := a δ y - a δ z - b z x^κ /. x -> 10000 x;
g[x_, y_, z_] := (bifP α)/Ocean + (γ a)/Ocean z - (γ a + c)/Ocean y /. x -> 10000 x;
h[x_, y_, z_] := bifP β/Land - bifP β/(M Land) x^κ + ((γ + δ μ) a)/
Land y - ((γ + δ μ) a + c)/Land z - (μ b z x^κ)/Land /. x -> 10000 x;
j = Outer[D, {f[x, y, z], g[x, y, z], h[x, y, z]}, {x, y, z}];
odeSystem = {f[x, y, z], g[x, y, z], h[x, y, z]};
SetPrecision[stabilityPoints = Solve[odeSystem == 0, {x, y, z}], 10];
stablePoint = stabilityPoints[[1]]
m = j /. stablePoint;
time = 600;
(*******Numerical Solution of ODE****************************)
Solution = SetPrecision[
NDSolve[{v'[t] == a δ to[t] - a δ τ[t] - b τ[t] (10000 v[t])^κ ,
to'[t] == (bifP α)/Ocean + (γ a)/Ocean τ[t] - (γ a + c)/Ocean to[t],
τ'[t] == bifP β/Land (1 - (10000 v[t])^κ/M) + ((γ + μ) a)/
Land to[t] - ((γ + μ) a + c)/Land τ[t] - (b μ)/Land (10000 v[t])^κ τ[t],
v[0] == x, to[0] == y, τ[0] == z}, {v[t], to[t], τ[t]},
{t, time}  , Method -> "StiffnessSwitching",
PrecisionGoal -> 8, AccuracyGoal -> 9], 30] /. stablePoint;
(***************Plots****************)
Plot[Evaluate[{v[t], to[t], τ[t]} /. Solution], {t, 0, time}, PlotRange -> All,
PlotStyle -> {{Green, Thick}, {Blue, Thick}, {Red, Thick}},
FrameLabel -> {Style["t", 16, Italic],
Style[Row[{Style["x", Italic], ", ", Style["y", 16, Italic]}], 16]},
Frame -> True, GridLines -> Automatic, AspectRatio -> 1,
PlotPoints -> time, ImageSize -> .9 400 {1, 1},
ImagePadding -> {{40, Automatic}, {40, Automatic}}]

ParametricPlot3D[Evaluate[{v[t], τ[t], to[t]} /. Solution], {t, 0, time},
BoxRatios -> {1, 1, 1}, PlotRange -> Automatic,
PlotStyle -> Blue, ImageSize -> 300 , PlotPoints -> time,
AxesLabel -> {Style["V", Large], Style["Tl", Large],
Style["To", Large], ImageSize -> 500, Axes -> True}]
(***The same system but with random pertubations, white noize precisely**)
sigma = 1
trajectory = ItoProcess[{
\[DifferentialD] v[t] == (a*δ*to[t] - a*δ*τ[t] - b*τ[t]*(10000 v[t])^κ )
\[DifferentialD] t + sigma \[DifferentialD] w[t],
\[DifferentialD] to[t] == ((bifP*α)/Ocean + (γ*a)/Ocean τ[t] - (γ*a + c)/Ocean to[t] )
\[DifferentialD] t + sigma \[DifferentialD] w[t],
\[DifferentialD] τ[t] == (bifP β/Land*(1 - (10000 v[t])^κ/M) + ((γ + μ)*a)/
Land to[t] - ((γ + μ)*a + c)/Land τ[t] - (b*μ)/Land (10000 v[t])^κ*τ[t])
\[DifferentialD] t + sigma \[DifferentialD] w[t]},
{ v[t], to[t], τ[t]}, {{v, to, τ}, {x, y, z}}, t, w \[Distributed] WienerProcess[]]
RandomFunction[trajectory, {0, 100, 1}];
ListLinePlot[%, Filling -> Axis]

• I'd love to help on an SDE question but your code is very messy. You're using Subscripts which isn't a good idea, you're being a bit excessive on 1000 digit precision, you're using a lot of greek characters, line breaks which didn't paste properly and you do x -> 10000 x which is confusing. Can you give us the final SDE to solve after everything has been properly evaluated or a shorter, cleaner version of your code? May 17, 2015 at 18:59
• I reformatted your code and corrected several syntax errors that may have arisen when you transferred the code to StackExchange. I did not, however, correct two coding errors. You have γa where I believe you meant γ*a, and  /. arr should be inside NDSolve, not outside, so that it assigns values to x, y, and z before solving the ODEs. With these corrections, I was able to generate your two plots. However, I am no expert in ItoProcess, so I cannot help you there. Please check your ItoProcess arguments, because I may have copied them incorrectly due to Unicode issues. May 18, 2015 at 4:50
• Thank you , guys, for your comments and editing the code, really appreciate the help. I'm such a typo and now still struggling with getting plot of pertubated system since I don't know how to use ItoProcess properly.
– Elka
May 18, 2015 at 7:51
• It is not clear what is arr? May 18, 2015 at 8:53
• Sorry I changed it to stablePoint but not everywhere in a code, fixed it. It is the stable point, or initial conditions for NDSolve{ v[0]=x, tau[0]=y, tl[0]=z}
– Elka
May 18, 2015 at 9:21

As noted in my comment, {x, y, z} must be defined in order for ItoProcess to work. Even then, the solution is unstable for the parameters chosen in the question. Stable parameters are, for instance, sigma = .01 and a RandomFunction step size of 0.1. In all, the modified code

sigma = .01;
trajectory =
ItoProcess[{\[DifferentialD]v[t] == (a*δ*to[t] - a*δ*τ[t] -
b*τ[t]*(10000 v[t])^κ) \[DifferentialD]t + sigma \[DifferentialD]w[t],
\[DifferentialD]to[t] == ((bifP*α)/Ocean + (γ*a)/Ocean τ[t] - (γ*a + c)/
Ocean to[t]) \[DifferentialD]t + sigma \[DifferentialD]w[t],
\[DifferentialD]τ[t] == (bifP β/Land*(1 - (10000 v[t])^κ/M) + ((γ + μ)*
a)/Land to[t] - ((γ + μ)*a + c)/Land τ[t] - (b*μ)/Land (10000 v[t])^κ*τ[t])
\[DifferentialD]t + sigma \[DifferentialD]w[t]},
{v[t], to[t], τ[t]}, {{v, to, τ}, {x, y, z} /. stablePoint}, t,
w \[Distributed] WienerProcess[]];
ans = RandomFunction[trajectory, {0, time, .1}];
ListLinePlot[ans, Filling -> Axis, PlotRange -> All]


This can be compared with the corresponding noise-free solution obtained from NDSolve.

An alternative way to display the results is

ListPointPlot3D[ans["Values"], BoxRatios -> {1, 1, 1}, PlotRange -> Automatic,
PlotStyle -> Blue, ImageSize -> 300 , AxesLabel -> {Style["V", Large],
Style["Tl", Large], Style["To", Large], ImageSize -> 500, Axes -> True}]


which can be compared with the noise-free ParametricPlot3D mentioned in the question.

• bbgodfrey, thank you soooooo much!!!!!!!!!! You made my day, seriously, thank you!
– Elka
May 18, 2015 at 18:33
• I hope you remain active in Mathematica.SE! I suggest: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. May 18, 2015 at 18:36
• Ok, I hope to stay active and even useful here :) 2) I don't have enough reputation to vote the answer up, but did accepted it.
– Elka
May 18, 2015 at 19:18
• A year ago, my situation on this site was no different from yours. It is amazing how much a person can learn asking and answering questions. Thanks for accepting the answer. Best wishes. May 18, 2015 at 19:21
• Thank you, very inspiring! One more question. I am trying to get 3D plot so I can observe cycles and see what happens according to the noise intensity, by smth like this: Plot3D[trajectory, {t, time}] or , which doesn't work. How do you work
– Elka
May 18, 2015 at 20:33