How to extract state vector after NDSolve for subsequent re-use?

Question

I'm using code like this:

DriveMotion = {
    r'[t] == Sin[γ[t]] v[t],
    v'[t] == -Sin[γ[t]] μ/r[t]^2 + Aion,
    γ'[t] == Cos[γ[t]]/r[t] (v[t] - μ/(v[t] r[t])),
    ϕ'[t] == Cos[γ[t]]/r[t] v[t]
    };
Init1 = {
    r[0] == Rinit,
    v[0] == Vcirc,
    γ[0] == 0, ϕ[0] == 0
    };

Trajectory1 = NDSolve[Join[DriveMotion, Init1], {r, v, γ, ϕ}, {t, 0, Tdrive}, 
   MaxSteps -> 5000000, Method -> "Automatic"];

to calculate the trajectory of a hypothetical ion-drive spacecraft. Now I want to extract the end state, i.e. after Tdrive seconds. Currently my approach is like this:

State = {r[t], v[t], γ[t], ϕ[t]};
Trajectory1EOD = Evaluate[State /. t -> Tdrive /. Trajectory1][[1]]
Init2 = {
    r[0] == Trajectory1EOD[[1]],
    v[0] == Trajectory1EOD[[2]],
    γ[0] == Trajectory1EOD[[3]],
    ϕ[0] == Trajectory1EOD[[4]]
    };

This works OK for a state vector with only four elements. But is there a way to accomplish this without copying each element of the state vector individually? Doing so with dozens of elements would be error-prone, at best! (This is my first question here, BTW. Apologies in advance if it is somehow not properly expressed.)

Answer 1

Maybe tutorial/NDSolveStateData can help. It lets you allow NDSolve to keep track of all the information the OP is asking about.

The first step is to call NDSolve`ProcessEquation instead of NDSolve; the same arguments may be given, but note that an interval of integration is optional. Since the OP did not include a value for Tdrive, I left it out. I also gave dummy values to the other parameters. An NDSolve`StateData data structure is returned.

DriveMotion = {r'[t] == Sin[γ[t]] v[t], 
   v'[t] == -Sin[γ[t]] μ/r[t]^2 + Aion, γ'[t] == Cos[γ[t]]/r[t] (v[t] - μ/(v[t] r[t])),
   ϕ'[t] == Cos[γ[t]]/r[t] v[t]};
Init1 = {r[0] == Rinit, 
   v[0] == Vcirc, γ[0] == 0, ϕ[0] == 0};

Block[{Rinit = 1, Vcirc = 1, Aion = 1, μ = 1},
 {state} = 
  NDSolve`ProcessEquations[{DriveMotion, Init1},
   {r, v, γ, ϕ}, {t}, MaxSteps -> 5000000, 
   Method -> "Automatic"]
 ]
(*  {NDSolve`StateData["<" 0. ">"]}  *)

To integrate up to t == 2, we call NDSolve`Iterate on the state with the stopping time.

NDSolve`Iterate[state, 2]

It alters the state. We can get the values of the state variables by calling NDSolve`ProcessSolutions with the argument "Forward" (use "Backward" if the value t == 2 is before the time of the initial condition in Init1).

NDSolve`ProcessSolutions[state, "Forward"]
(*
  {r[2.] -> 3.19208, Derivative[1][r][2.] -> 2.43266,
   v[2.] -> 2.6337, Derivative[1][v][2.] -> 0.90935,
   γ[2.] -> 1.17754, Derivative[1][γ][2.] -> 0.301888,
   ϕ[2.] -> 1.85649, Derivative[1][ϕ][2.] -> 0.316168}
*)

Now we can continue to integrate, say up to t == 4:

NDSolve`Iterate[state, 4]

To get the complete solution, we can use NDSolve`ProcessSolutions with no extra argument:

sol = NDSolve`ProcessSolutions[state]

The variables used are stored in state and you can get at the pieces without having to type each one out:

vars = state@"Variables"
(*  {t, {}, {r, v, γ, ϕ},
     {Derivative[1][r], Derivative[1][v], Derivative[1][γ], Derivative[1][ϕ]},
     {}, {}, {}, {}}  *)

One can extract the dependent variables:

funcs = NDSolve`SolutionDataComponent[vars, "DependentVariables"]
(*  {r, v, γ, ϕ}  *)

Apply the solution function to the time t == 2. and compare with the values where we stopped the integration above.

Through[funcs[2.]] /. sol
(*  {3.19208, 2.6337, 1.17754, 1.85649}  *)

More arguments and options are described in the tutorial linked at the beginning.

Addendum -- Reinitialization

Another glance at the OP's question makes me think that the OP might be interested in reinitializing the ODE at t == 0 with the values of the state variables obtained at t == 2. In the OP's example, which is autonomous, it makes no difference whether we start over at t == 0 and continue for another 2 or continue to t == 4 as above. If the coefficients depended on t, then it would make a difference. In either case, one can start over with the last state as the initial state as follows.

Suppose therefore we have just completed the step NDSolve`Iterate[state, 2] above and that funcs is a list of the dependent variables as above. One can use NDSolve`Reinitialize with new initial conditions to get a newstate data structure.

{newstate} = NDSolve`Reinitialize[state, 
 Through[funcs[0]] == Through[funcs[2.]] /. NDSolve`ProcessSolutions[state, "Forward"]];

Now we compare advancing the newstate from t == 0 to 2 and the old state from t == 2 to 4.

NDSolve`Iterate[newstate, 2.]

NDSolve`Iterate[state, 4.]

Compare values:

NDSolve`ProcessSolutions[state, "Forward"]
NDSolve`ProcessSolutions[newstate, "Forward"]
(*
{r[4.`] -> 10.195149962394567`, Derivative[1][r][4.`] -> 4.532956349402`,
 v[4.`] -> 4.567362769921008`,  Derivative[1][v][4.`] -> 0.9904516398189146`,
 γ[4.`] -> 1.4479746578634392`, Derivative[1][γ][4.`] -> 0.054627031974180444`,
 ϕ[4.`] -> 2.1329726485832348`, Derivative[1][ϕ][4.`] -> 0.05488509738922278`}

{r[2.`] -> 10.195149985177293`, Derivative[1][r][2.`] -> 4.5329563552874`,
 v[2.`] -> 4.567362772976815`,  Derivative[1][v][2.`] -> 0.9904516398569216`,
 γ[2.`] -> 1.4479746616356868`, Derivative[1][γ][2.`] -> 0.054627029636341445`,
 ϕ[2.`] -> 2.1329726759733605`, Derivative[1][ϕ][2.`] -> 0.05488509503939425`}
*)

A note on funcs

I usually make my own list of funcs at the very beginning: At some point you type the list {r, v, γ, ϕ} of functions to be solved for. Why not store it in a variable for later use?

Answer 2

Here I present my first working solution that uses the components of NDSolv. (I hope posting a new answer is the right way to do this.) The goal is to use the "ion drive" for awhile, and then "coast" from where that leaves us up to the vicinity of the Moon.

This first block initializes some constants:

G = 6.672*10^-11;
M = 5.97219*10^24;(*Mass of Earth*)
μ = G M;
Rearth = 6.378*10^6;
Rmoon = 384.4*10^6;
Rinit = Rearth + 400000; (*Radius of Earth plus height of orbit*)

Vcirc = Sqrt[μ/Rinit];(*Circular orbital velocity*)

SecondsFromDays = 60*60*24*days; 
Tdrive = Evaluate[SecondsFromDays /. days -> 75];
Tcoast = Evaluate[SecondsFromDays /. days -> 30];

The next block establishes the equations of motion. There's some repetition here I would prefer to avoid. Sorry if I'm just being dense, but neither of the two ways I tried yielded what was required. Hints on this would be greatly appreciated!(See the comment at the bottom.)

StateFunctions = {r[t], v[t], γ[t], ϕ[t]};
Motion = {r'[t], 
    v'[t], γ'[t], ϕ'[
     t]} == {Sin[γ[t]] v[t], -Sin[γ[t]] μ/r[t]^2 + 
     Aion, Cos[γ[t]]/r[t] (v[t] - μ/(v[t] r[t])), 
    Cos[γ[t]]/r[t] v[t]};

Here we generate the trajectory while under drive:

DriveMotion = Motion /. Aion -> 0.001;
DriveInit = StateFunctions == {Rinit, Vcirc, 0, 0} /. t -> 0;
{Drive} = 
  NDSolve`ProcessEquations[{DriveMotion, DriveInit}, {r, 
    v, γ, ϕ}, {t}, MaxSteps -> 5000000];
MotionVariables = 
  NDSolve`SolutionDataComponent[Drive@"Variables", 
   "DependentVariables"];
NDSolve`Iterate[Drive, Tdrive]
DriveTrajectory = NDSolve`ProcessSolutions[Drive];

Next we snag the end state of the drive trajectory and make it the initial state of the coast trajectory. I simply couldn't figure out how to do this without 'Hold' and 'ReleaseHold'. Again, hints would be appreciated! We then calculate the coast trajectory.

CoastMotion = Motion /. Aion -> 0;
CoastInit = 
  Hold[StateFunctions] == Through[MotionVariables[Tdrive]] /. 
   DriveTrajectory;
NDSolve`ProcessSolutions[Drive, "Forward"];
{Coast} = 
  NDSolve`ProcessEquations[{CoastMotion, 
    ReleaseHold[CoastInit] /. t -> 0}, {r, v, γ, ϕ}, {t}, 
   MaxSteps -> 5000000];
NDSolve`Iterate[Coast, Tcoast]
CoastTrajectory = NDSolve`ProcessSolutions[Coast];

Finally we plot these, along with the orbit of the Moon:

OrbitOfMoonXY = {Rmoon Cos[2 π θ], 
   Rmoon Sin[2 π θ]};
Cartesian = {r[t] Cos[ϕ[t]], r[t] Sin[ ϕ[t]]};
ParametricPlot[{
  Evaluate[OrbitOfMoonXY /. θ -> Range],
  Evaluate[Cartesian /. DriveTrajectory /. t -> Range*(Tdrive)],
  Evaluate[Cartesian /. CoastTrajectory /. t -> Range*(Tcoast)]
  }, {Range, 0, 1}, PlotLabel -> "Trajectory", PlotRange -> All]

Tada! We see the green spiral of the drive trajectory leading seamlessly onto the red ellipse of the coast trajectory! I am sure this can all be done more simply/tersely, and as I plan to use this example as the basis for further simulations I welcome any suggestions about how to improve this. Thanks again to all those who helped me get this far!

Comment -- This can be accomplished with:

StateVariables = {r, v, γ, ϕ};
StateDerivatives = Through[Map[Derivative[1], StateVariables][t]];