# Why is solution of NDSolve different than Heun's method?

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?

• 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] ] Feb 6, 2022 at 1:04
• Try e.g. h=0.04 for Heun's method. Feb 6, 2022 at 4:22

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


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