4
$\begingroup$

I am thrying to solve the Lorrentz system using both the Heun's numerical method given by the following code

Clear["Global`*"]
S[a_, b_, c_, h_, N_] := (x[0] = a; y[0] = b; z[0] = c; 
  Do[{x[n + 1] = 
     x[n] + h*
       f[x[n] + h/2*f[x[n], y[n], z[n]], 
        y[n] + h/2*g[x[n], y[n], z[n]], 
        z[n] + h/2*p[x[n], y[n], z[n]]], 
    y[n + 1] = 
     y[n] + h*
       g[x[n] + h/2*f[x[n], y[n], z[n]], 
        y[n] + h/2*g[x[n], y[n], z[n]], 
        z[n] + h/2*p[x[n], y[n], z[n]]], 
    z[n + 1] = 
     z[n] + h*
       p[x[n] + h/2*f[x[n], y[n], z[n]], 
        y[n] + h/2*g[x[n], y[n], z[n]], 
        z[n] + h/2*p[x[n], y[n], z[n]]]}, {n, 0, N}])

f[x_, y_, z_] = s*(y - x);
g[x_, y_, z_] = x*(r - z) - y;
p[x_, y_, z_] = x*y - b*z;
s = 5;
r = 28;
b = 8/3;

S[0.001, 0.001, 0.001, 0.05, 100/0.05]
X = Interpolation[Table[{n, x[n]}, {n, 0, 100/0.05}]]
Y = Interpolation[Table[{n, y[n]}, {n, 0, 100/0.05}]]
Z = Interpolation[Table[{n, z[n]}, {n, 0, 100/0.05}]]


fig1 = ParametricPlot3D[{X[n], Y[n], Z[n]}, {n, 0, 500}]

and Mathematica's NDSolve

Clear["Global`*"]
f[a_, b_, c_] := 
 NDSolve[{x[0] == a, y[0] == b, z[0] == c, x'[t] == 5*(y[t] - x[t]), 
   y'[t] == x[t]*(28 - z[t]) - y[t], 
   z'[t] == x[t]*y[t] - 8/3*z[t]}, {x[t], y[t], z[t]}, {t, 0, 100}]
s1 = f[0.001, 0.001, 0.001]
ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. s1], {t, 0, 50}, 
 PlotRange -> All]
s2 = f[0.00101, 0.00101, 0.00101]

Plot[{x[t] /. s1, x[t] /. s2}, {t, 0, 50}]
Plot[{y[t] /. s1, y[t] /. s2}, {t, 0, 50}]
Plot[{z[t] /. s1, z[t] /. s2}, {t, 0, 50}]

However the solutions are a lot different for those two methods. Why is this happening? Is there a problem why my code or with the parameters (like time step) of Heun's method?

$\endgroup$
2
  • 2
    $\begingroup$ You could check your Heun's method code against Heun's method in NDSolve: HeunCoefficients[2, p_] := N[{{{1}}, {1/2, 1/2}, {1}}, p]; Block[{h = 0.05}, sol = NDSolve[ ..., Method -> {"FixedStep", Method -> {"ExplicitRungeKutta", "Coefficients" -> HeunCoefficients, "DifferenceOrder" -> 2}}, StartingStepSize -> h] ] $\endgroup$
    – Michael E2
    Feb 6, 2022 at 1:04
  • $\begingroup$ Try e.g. h=0.04 for Heun's method. $\endgroup$
    – xzczd
    Feb 6, 2022 at 4:22

1 Answer 1

6
$\begingroup$

There is typos in the code with N using as a limit in Do loop and with X,Y,Z definitions. Also we need to localize variables in NDSolve with using Module. The code for comparison two methods can be written as follows

Clear["Global`*"]
f[x_, y_, z_] := s*(y - x);
g[x_, y_, z_] := x*(r - z) - y;
p[x_, y_, z_] := x*y - b*z;
S[x0_, y0_, z0_, h_, nn_] := (x[0] = x0; y[0] = y0; z[0] = z0;
  Do[{x[n + 1] = 
     x[n] + h*
       f[x[n] + h/2*f[x[n], y[n], z[n]], 
        y[n] + h/2*g[x[n], y[n], z[n]], 
        z[n] + h/2*p[x[n], y[n], z[n]]], 
    y[n + 1] = 
     y[n] + h*
       g[x[n] + h/2*f[x[n], y[n], z[n]], 
        y[n] + h/2*g[x[n], y[n], z[n]], 
        z[n] + h/2*p[x[n], y[n], z[n]]], 
    z[n + 1] = 
     z[n] + h*
       p[x[n] + h/2*f[x[n], y[n], z[n]], 
        y[n] + h/2*g[x[n], y[n], z[n]], 
        z[n] + h/2*p[x[n], y[n], z[n]]]}, {n, 0, nn}])





    HeunCoefficients[2, par_] := 
 N[{{{1}}, {1/2, 1/2}, {1}}, par]; op = {"ExplicitRungeKutta", 
  "Coefficients" -> HeunCoefficients, "DifferenceOrder" -> 2}; 
S1[x0_, y0_, z0_, tmax_] := 
 Module[{tm = tmax, x, y, z}, 
  sol = NDSolveValue[{x[0] == x0, y[0] == y0, z[0] == z0, 
     x'[t] == s*(y[t] - x[t]), y'[t] == x[t]*(r - z[t]) - y[t], 
     z'[t] == x[t]*y[t] - b*z[t]}, {x[t], y[t], z[t]}, {t, 0, tm}, 
    Method -> op, StartingStepSize -> h]; sol]

s = 5;
r = 28;
b = 8/3; h = 1/200; tmax = 100; nmax = Round[tmax/h];

S[0.001, 0.001, 0.001, h, nmax];
X = Interpolation[Table[{h n, x[n]}, {n  , 0, nmax}]];
Y = Interpolation[Table[{h n, y[n]}, {n, 0, nmax}]];
Z = Interpolation[Table[{h n, z[n]}, {n, 0, nmax}]];  
s1 = S1[0.001, 0.001, .001, 100];  

Visualization

{ParametricPlot3D[Evaluate[s1], {t, 0, tmax}, PlotRange -> All, 
  ColorFunction -> Hue, PlotLabel -> "NDSolve", 
  PlotTheme -> "Scientific"], 
 ParametricPlot3D[Evaluate[{X[t], Y[t], Z[t]}], {t, 0, tmax}, 
  PlotRange -> All, ColorFunction -> Hue, 
  PlotLabel -> "Heun's method", PlotTheme -> "Scientific"]}

Figure 1 Even two solutions look identical, nevertheless there is a difference of h order at t about 1 as it shown below

Plot[Norm[s1 - {X[t], Y[t], Z[t]}], {t, 0, tmax}, PlotPoints -> 200, 
 PlotRange -> All]

Figure 2

$\endgroup$

Your Answer

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