3
$\begingroup$

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.5`100;
    κ = 0.7`100;
    α = 0.5`100;
    tl = 270;
    too = 280;
    Ocean = 500;
    Land = 25;
    γ = 2.256`100;
    μ = 0.335`100;
    k = 5100;
    M = 0.9999`100 10^4;
    vo = 0.746326`100 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]    
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – Histograms May 17 '15 at 18:59
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey May 18 '15 at 4:50
  • $\begingroup$ 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. $\endgroup$ – Elka May 18 '15 at 7:51
  • $\begingroup$ It is not clear what is arr? $\endgroup$ – chris May 18 '15 at 8:53
  • $\begingroup$ 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} $\endgroup$ – Elka May 18 '15 at 9:21
5
$\begingroup$

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]

enter image description here

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

enter image description here

Addendum

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}]

enter image description here

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

enter image description here

$\endgroup$
  • $\begingroup$ bbgodfrey, thank you soooooo much!!!!!!!!!! You made my day, seriously, thank you! $\endgroup$ – Elka May 18 '15 at 18:33
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey May 18 '15 at 18:36
  • $\begingroup$ 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. $\endgroup$ – Elka May 18 '15 at 19:18
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey May 18 '15 at 19:21
  • $\begingroup$ 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 $\endgroup$ – Elka May 18 '15 at 20:33

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.