# Errors of the definite integral term on the code

 δ = 10; ρ = 28; β = (8/3); x0 = -7; y0 = 0; z0 = 4;

s = NDSolve[{x'[t] == δ*(y[t] - x[t]), y'[t] == ρ*x[t] - y[t] - x[t] z[t],
z'[t] == x[t] y[t] - β*z[t], x == x0, y == y0, z == z0},
{x, y, z}, {t, 0, 50}, WorkingPrecision -> 16, MaxSteps -> Infinity]

Plot[{x[t] /. s, y[t] /. s, z[t] /. s}, {t, 0, 10},
AxesLabel -> Automatic, PlotLegends -> "Expressions"]

p[t_?NumericQ] := NIntegrate[δ ((y[tt] /. s) - (x[tt] /. s)), {tt, 0, t},
WorkingPrecision -> 16, AccuracyGoal -> 100]

q[t_?NumericQ] := NIntegrate[ρ (x[tt] /. s) - (y[tt] /. s) - (x[tt] /. s)*(z[tt] /. s),
{tt, 0, t}, WorkingPrecision -> 16, AccuracyGoal -> 100]

w[t_?NumericQ] := NIntegrate[(x[tt] /. s)*(z[tt] /. s) - βz[tt] /. s, {tt, 0, t},
WorkingPrecision -> 16, AccuracyGoal -> 100]

res[t_] := (x[t] - x0) + (y[t] - y0) + (z[t] - z0) - (p[t] + q[t] + w[t]);

Plot[Evaluate[RealExponent[res[t] /. s]], {t, 0, 4},
PlotStyle -> {GrayLevel, RGBColor[2, 0, 0]}, PlotPoints -> 10]


I am trying to generate the residual function graphs for Lorentz system but code has some errors on multiplication of function from definite integral term.

• \[Beta]z[tt] is an undefined quantity. Perhaps, you mean \[Beta] z[tt]. With this correction, the code is painfully slow. – bbgodfrey Feb 28 '19 at 11:42
• The definition of w is incorrect. – bbgodfrey Feb 28 '19 at 12:30
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) 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, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Feb 28 '19 at 12:58

## 2 Answers

To obtain literally the quantity requested in the question (with errors in w corrected) in a reasonable amount of time, use

Clear[x, y, z, p, q, w]
{x, y, z} = NDSolveValue[{x'[t] == \[Delta]*(y[t] - x[t]),
y'[t] == \[Rho]*x[t] - y[t] - x[t] z[t], z'[t] == x[t] y[t] - \[Beta]*z[t],
x == x0, y == y0, z == z0}, {x, y, z}, {t, 0, 10}];


to determine {x, y, z} essentially as in the question, and then use NDSolve again but with higher precision to determine {p, q, w} quickly.

{p, q, w} = NDSolveValue[{p'[t] == \[Delta] (y[t] - x[t]),
q'[t] == \[Rho] x[t] - y[t] - x[t]*z[t], w'[t] == x[t]*y[t] - \[Beta] z[t],
p == 0, q == 0, w == 0}, {p, q, w}, {t, 0, 10},
WorkingPrecision -> 30, MaxSteps -> 100000];

res[t_] := (x[t] - x0) + (y[t] - y0) + (z[t] - z0) - (p[t] + q[t] + w[t]);
Plot[Evaluate[RealExponent[res[t]]], {t, 0, 10}, PlotStyle -> Red,
ImageSize -> Large, AxesLabel -> {t, "res"}, LabelStyle -> {Bold, Black, 15}] • Thanks @bbgodfrey much appreciated – Malesela kekana Mar 1 '19 at 14:35
• This is close to what I am looking for much appreciated – Malesela kekana Mar 1 '19 at 14:38

If I understand your question you try to solve an ode-system and plot the residual?

If so, there is no need to define p[t],...!

Try

Clear[x, y, z]
{x, y, z} =NDSolveValue[{x'[t] == \[Delta]*(y[t] - x[t]),
y'[t] == \[Rho]*x[t] - y[t] - x[t] z[t],
z'[t] == x[t] y[t] - \[Beta]*z[t], x == x0, y == y0,
z == z0}, {x, y, z}, {t, 0, 50} ];

Plot[{x[t], y[t], z[t]}, {t, 0, 50}, AxesLabel -> Automatic] Now plot the three residuals

Plot[Evaluate[ {x'[t] == \[Delta]*(y[t] - x[t]),y'[t] == \[Rho]*x[t] - y[t] -x[t] z[t], z'[t] == x[t] y[t] - \[Beta]*z[t]} /. Equal -> Subtract], {t, 0,50}, PlotRange-> {-.01, .01}] That's it!

• Thanks Ulrich the idea is not to plot individual residuals. The idea is to plot the sum of those three residual function – Malesela kekana Mar 1 '19 at 14:36
• Ok, it is no problem to easy modify the plot and sum the Abs of the individual residuals. One point I want to note is your definition of res[t]=... If NIntegrate works without error, you get p[t]=x[t],q[t]=y[t],w[t]=z[t] and res==x0+y0+z0 is the sum of initial conditions. Curios! In my answer the residual is the error obtained by substituting the solution of NDSolve into the ode. – Ulrich Neumann Mar 1 '19 at 15:02