2
$\begingroup$
 δ = 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.

$\endgroup$
  • $\begingroup$ \[Beta]z[tt] is an undefined quantity. Perhaps, you mean \[Beta] z[tt]. With this correction, the code is painfully slow. $\endgroup$ – bbgodfrey Feb 28 at 11:42
  • $\begingroup$ The definition of w is incorrect. $\endgroup$ – bbgodfrey Feb 28 at 12:30
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey Feb 28 at 12:58
2
$\begingroup$

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

enter image description here

$\endgroup$
  • $\begingroup$ Thanks @bbgodfrey much appreciated $\endgroup$ – Malesela kekana Mar 1 at 14:35
  • $\begingroup$ This is close to what I am looking for much appreciated $\endgroup$ – Malesela kekana Mar 1 at 14:38
1
$\begingroup$

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]

enter image description here

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

enter image description here

That's it!

$\endgroup$
  • $\begingroup$ Thanks Ulrich the idea is not to plot individual residuals. The idea is to plot the sum of those three residual function $\endgroup$ – Malesela kekana Mar 1 at 14:36
  • $\begingroup$ 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. $\endgroup$ – Ulrich Neumann Mar 1 at 15:02

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.