Jim's Orbit Problem

Copied from Dormand p163-165 #12

As two second order equations!

g[{x_, y_, vx_, vy_}]  :=  
      Module[
{r1, 
          r2}, 
{
r1 = 
            Sqrt[(((x + mm))^2) + (y^2)]; 
r2 = 
            Sqrt[(((x - ee ))^2) + (y^2)]; 
          2 vy  +  x  -  ee (x + mm)/(r1^3) -  
            mm  (x - ee)/(r2^3), 
-2 vx  +  
            y  -  ee y/(r1^3) -  
            mm  y/(r2^3)
}
]

FastG =  
      Compile[
{x, y, vx, 
          vy}, 
{r1 = 
            Sqrt[(((x + mm))^2) + (y^2)]; 
r2 = 
            Sqrt[(((x - ee ))^2) + (y^2)]; 
          2 vy  +  x  -  ee (x + mm)/(r1^3) -  
            mm  (x - ee)/(r2^3), 
-2 vx  +  
            y  -  ee y/(r1^3) -  
            mm  y/(r2^3)
}
]

FastG[2.1, 2.3, 3.1, 4.5]
g[{2.1, 2.3, 3.1, 4.5}]

Timing[ Table[
          g[{Random[], Random[], Random[], 
              Random[]}], {(10^4)}];]
Timing[ Table[
          FastG[Random[], Random[], Random[], 
            Random[]], {(10^4)}];]

A  =  Table[ Random[], {(10^7)}];

Timing[
Table[
A[[ 
            Random[Integer, {1, 
                (10^7)}]]], 
{(10^4)}];]

mm =  0.012277471;
ee  =  1 -  mm;
Temp  = 
        NDSolve[
{r1 = 
              Sqrt[(((x[t] + mm))^2) + y[t]^2]; 
r2 = 
              Sqrt[(((x[t] - ee ))^2) + 
                  y[t]^2]; 
x''[t]  ==  
              2 y'[t]  +  x[t]  -  ee (x[t] + mm)/(r1^3) -  
                mm  (x[t] - ee)/(r2^3), 
y''[
                t]  ==  -2 x'[t]  +  y[t]  -  
                ee y[t]/(r1^3) -  mm  y[t]/(r2^3), 
x[
                0] == 1.15, 
y[0] ==  
              0, 
x'[0] == 
              0, 
y'[0] == 
              0.0086882909
}, 
{x[t], 
            y[t]}, 
{t, 0,  
            29.4602}, 
MaxSteps -> 
            (10^6)];

{xSol[t_], ySol[t_]} = {x[t], y[t]} /. 
          Temp[[1]];

ParametricPlot[ {xSol[t], ySol[t]}, {t, 0,  
          29.4602}, 
AspectRatio -> 
          Automatic, 
PlotRange -> 
          All
];

As four first order equations!

Clear[y, Derivative]

mm =  0.012277471;

    ee  =  1 -  mm;

    Temp  = 
        NDSolve[
{r1 = 
              Sqrt[(((x[t] + mm))^2) + y[t]^2]; 
r2 = 
              Sqrt[(((x[t] - ee ))^2) + 
                  y[t]^2]; 
vx'[t]  ==  
              2 vy[t]  +  x[t]  -  ee (x[t] + mm)/(r1^3) -  
                mm  (x[t] - ee)/(r2^3), 
vy'[
                t]  ==  -2 vx[t]  +  y[t]  -  
                ee y[t]/(r1^3) -  
                mm  y[t]/(r2^3), 
x'[t] == 
              vx[t], 
y'[t] == 
              vy[t], 
vx[0] == 
              0, 
vy[0] == 
              0.0086882909, 
x[0] == 
              1.15, 
y[0] ==  
              0
}, 
{vx[t], vy[t], x[t], 
            y[t]}, 
{t, 0,  
            29.4602}, 
MaxSteps -> 
            (10^6)];


    {xSol[t_], ySol[t_]} = {x[t], y[t]} /. 
          Temp[[1]];


    Fred = 
        ParametricPlot[ {xSol[t], ySol[t]}, {t, 0,  
            29.4602}, 
AspectRatio -> 
            Automatic
];

Table[
Show[
Fred, 
Graphics[
{Hue[0], PointSize[0.05], 
            Point[{xSol[t], 
                ySol[t]}]
}
]
], 
{t, 0,  24.6, 
        0.1}
]

ySol[t]

ts  = xSol[t][[0, 3, 
          1]];
hs  =  
        Drop[RotateLeft[ts] - ts, -1];
\(\(ListPlot[
        Log[10, hs], 
PlotRange -> 
          All, 
AxesLabel -> {"n", \*""(Subscript[Log, 10])h""}
];\)\)

Here is the step sizes!

\(\(ListPlot[
        Log[10, Take[hs, 100]], 
PlotRange -> 
          All, 
AxesLabel -> {"n", \*""(Subscript[Log, 10])h""}
];\)\)

Robertson's Stiff Problem

Copied from Dormand p198

Clear[α, β, γ, y1, y2, y3, t, y1Sol, y2Sol, 
      y3Sol]

    
{α, β, γ} =  {0.04, 1.0×(10^4), 
        3.0×(10^7)}

    
Temp  = 
        NDSolve[
{y1'[
                t]  ==  -α y1[t] +  β y2[t] y3[
                    t], 
y2'[
                t] ==  α y1[t]  -  β y2[t] y3[
                    t] - γ y2[t]^2, 
y3'[
                t] == γ y2[t]^2, 
y1[
                0] == 1, 
y2[0] == 
              0, 
y3[0] ==  
              0
}, 
{y1[t], y2[t], 
            y3[t]}, 
{t, 0,  
            500}, 
Method -> 
            Automatic
  
  
, 
          MaxSteps -> 
            (10^6) ];

{y1Sol[t_], y2Sol[t_], y3Sol[t_]} = {y1[t], y2[t], y3[t]} /. 
        Temp[[1]]

Needs["Graphics`Graphics`"]
LogLinearListPlot[ 
Table[{(E^t), 
            y2Sol[E^t]}, {t, -12,  Log[500], 
            (10^(-1))}], 
PlotRange -> 
          All
];

ts  = y2Sol[t][[0, 3, 
          1]];
ListPlot[Log[10, Drop[ts - RotateRight[ts], 1]]]