6
$\begingroup$

I am beginning to read the Differential Equations Laboratory Workbook by Borelli and the preface begins with an image.

enter image description here

I understand how to use NDSolveValue (at a very beginning level), so I am wondering how I would use it to produce 10,000 points and duplicate this picture. I haven't yet learned how to use the Methods option. Perhaps one of the methods and a step size are available that will let me produce 10,000 points and a duplicate image.

$\endgroup$
  • $\begingroup$ Ah Lorenz… just use NDSolve[] like usual, and evaluate the component functions over 10000 equispaced points in your interval, using Table[] or Range[]. $\endgroup$ – J. M.'s technical difficulties May 9 '15 at 17:27
  • 1
    $\begingroup$ On second thought: I haven't used it myself, but this seems like the sort of situation that calls for ParametricNDSolve[]. $\endgroup$ – J. M.'s technical difficulties May 9 '15 at 17:33
  • $\begingroup$ @Guesswhoitis. Are you saying that instead of using {t,0,10}, you can use Table or Range? Can I ask specifically how? What would that look like? $\endgroup$ – David May 9 '15 at 22:32
  • $\begingroup$ Well, you still need the {t, 0, 10} in NDSolve[] or its parametric variant. You can then do something like Table[x1[t], {t, 0, 10, 10/(1*^4-1)}] to get equispaced points. $\endgroup$ – J. M.'s technical difficulties May 9 '15 at 23:09
8
$\begingroup$

Here's a slightly different approach, using ParametricNDSolve:

sols = ParametricNDSolve[
  {x1'[t] == -10 x1[t] + 10 x2[t],
   x2'[t] == a1 x1[t] - x2[t] - x1[t] x3[t],
   x3'[t] == -(8/3) x3[t] + x1[t] x2[t],
   x1[0] == 1,
   x2[0] == 1,
   x3[0] == 1},
  {x1, x2, x3},
  {t, 0, 10},
  {a1}
 ];

This produces ParametricFunction objects, which take a value of the parameter a1 first, but then can be plotted like any other function:

style = {12, FontFamily -> "Courier"};
lines = {
   Directive[Thickness[0.001], Black],
   Directive[Thickness[0.001] , Black, Dashing[0.005]]
  };
grid = Style[
   Grid[
    {{"x1'=-10x1+10x2", "x1(0):1.0", Item["t-zero 0", Alignment -> Right]},
     {"x2'=a1*x1-x2-x1*x3", "x2(0):1.0", Item["t-final 10", Alignment -> Right]},
     {"x3'=-(8/3)x3+x1*x2", "x3(0):1.0", Null}},
    Alignment -> Left, Spacings -> 4],
    Sequence @@ style];
legend = LineLegend[
   lines,
   {"{50}", "{55}"},
   LegendLabel -> "a1",
   LegendLayout -> {"Column", 1},
   LegendMarkerSize -> {40, 1},
   LabelStyle -> style
  ];

Labeled[
 Plot[
  Evaluate[(x1[#][t] & /@ {50, 55}) /. sols], {t, 0, 10},
  BaseStyle -> style,
  Axes -> False,
  FrameLabel -> {"t", "x1"},
  Frame -> {{True, False}, {True, False}},
  PlotStyle -> lines,
  PlotRange -> {{0, 10}, {-30, 30}},
  ImageSize -> 1.5 {5, 3} 72],
Row[{legend, grid}, "\t"], Top]

plot

Pretty close!

Aside: the precision of the solution is important here. If we use a lower precision,

sols = ParametricNDSolve[
  {x1'[t] == -10 x1[t] + 10 x2[t],
   x2'[t] == a1 x1[t] - x2[t] - x1[t] x3[t],
   x3'[t] == -(8/3) x3[t] + x1[t] x2[t],
   x1[0] == 1,
   x2[0] == 1,
   x3[0] == 1},
  {x1, x2, x3},
  {t, 0, 10},
  {a1},
  WorkingPrecision -> 8
 ];

we get a graph that diverges from the first after a few time units!

plot2

$\endgroup$
  • $\begingroup$ In attempting to run your code, I got this error: ReplaceAll::reps: {%18} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >> $\endgroup$ – David May 9 '15 at 19:59
  • $\begingroup$ @David, sorry: I had a typo, which I've fixed now. In the Plot, I had Out[18], which were my solutions, in place of sols, which is the name I gave them for this code. $\endgroup$ – Virgil May 9 '15 at 20:02
  • $\begingroup$ No need to be sorry, as this is a wonderful contribution. There is a lot to learn here. $\endgroup$ – David May 9 '15 at 20:13
  • $\begingroup$ What I am wondering now is the difference in the interval from about 8 to 10. Your solution and Bob Hanlon's solution match, but differ from the image I provided, which makes me think that Borelli may have used a different method? Euler, Runga-Kutta? Something like that? $\endgroup$ – David May 9 '15 at 20:15
  • $\begingroup$ @David, it is probably in the numerics - this diffEQ is really sensitive to initial conditions, so it makes sense that it will also be sensitive to the numerical precision. This can be specified with the option WorkingPrecision in the NDSolve family. Both Bob Hanlon and I used the default precision (with is MachinePrecision) but I briefly checked that the solutions looks much different if you take something like WorkingPrecision -> 8. So probably, Borelli might have been using lower precision (the age of the book might tell you). $\endgroup$ – Virgil May 9 '15 at 20:22
8
$\begingroup$
eqns = {
   x1'[t] == -10*x1[t] + 10*x2[t],
   x2'[t] == a1*x1[t] - x2[t] - x1[t]*x3[t],
   x3'[t] == -8*x3[t]/3 + x1[t]*x2[t],
   x1[0] == 1, x2[0] == 1, x3[0] == 1};

Plot[
 Evaluate[x1[t] /. Table[
    NDSolve[eqns, {x1[t], x2[t], x3[t]}, {t, 0, 10},
      MaxStepFraction -> 1/10000][[1]],
    {a1, {50, 55}}]],
 {t, 0, 10},
 PlotStyle -> {ColorData[106, 1],
   Directive[ColorData[106, 2], AbsoluteDashing[{2, 3}]]},
 Frame -> {{True, False}, {True, False}},
 FrameLabel -> {"t", "x1"},
 Axes -> False,
 PlotLegends -> {"a1=50", "a1=55"},
 PlotPoints -> 100]

enter image description here

$\endgroup$

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.