# 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[0] == x0, y[0] == y0, z[0] == 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[0], 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. Feb 28, 2019 at 11:42
• The definition of w is incorrect. Feb 28, 2019 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. Feb 28, 2019 at 12:58

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[0] == x0, y[0] == y0, z[0] == 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] == 0, q[0] == 0, w[0] == 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 Mar 1, 2019 at 14:35
• This is close to what I am looking for much appreciated Mar 1, 2019 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[0] == x0, y[0] == y0,
z[0] == 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 Mar 1, 2019 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. Mar 1, 2019 at 15:02