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

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.)

• Welcome to Mathematica.SE! I suggest that: 1) You 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. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey May 10 '15 at 18:46
• Not quite sure what you want, but does this work for you: MapIndexed[#1 == Trajectory1EOD[#2[[1]]] &, State /. t -> 0] – Sjoerd C. de Vries May 10 '15 at 19:21
• Have you taken a look at NDSolveStateData? – J. M.'s ennui May 10 '15 at 20:04
• @J. M. Do you not think it is a duplicate of mathematica.stackexchange.com/questions/5556/…? There are a few extra requests in this Q, but the main question seems the same. – Michael E2 May 10 '15 at 20:52
• @Michael, kind of on the fence here, as the other question is asking about PDEs; both problems (PDE and vector-valued ODE) just turn out to be tractable with the built-in continuation framework. – J. M.'s ennui May 10 '15 at 20:57

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 NDSolveProcessEquation 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 NDSolveStateData 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} =
NDSolveProcessEquations[{DriveMotion, Init1},
{r, v, γ, ϕ}, {t}, MaxSteps -> 5000000,
Method -> "Automatic"]
]
(*  {NDSolveStateData["<" 0. ">"]}  *)


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

NDSolveIterate[state, 2]


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

NDSolveProcessSolutions[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:

NDSolveIterate[state, 4]


To get the complete solution, we can use NDSolveProcessSolutions with no extra argument:

sol = NDSolveProcessSolutions[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 = NDSolveSolutionDataComponent[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.

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 NDSolveIterate[state, 2] above and that funcs is a list of the dependent variables as above. One can use NDSolveReinitialize with new initial conditions to get a newstate data structure.

{newstate} = NDSolveReinitialize[state,
Through[funcs[0]] == Through[funcs[2.]] /. NDSolveProcessSolutions[state, "Forward"]];


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

NDSolveIterate[newstate, 2.]

NDSolveIterate[state, 4.]


Compare values:

NDSolveProcessSolutions[state, "Forward"]
NDSolveProcessSolutions[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?

• Clearly as written elsewhere "the designers of NDSolve[] have precisely anticipated" my needs! ;-) I will soon be a happy user of the technique shown here (and at reference.wolfram.com/language/tutorial/NDSolveStateData.html). Thanks! – sdsds May 10 '15 at 23:17
• @sdsds You're welcome. Best, – Michael E2 May 10 '15 at 23:42

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} =
NDSolveProcessEquations[{DriveMotion, DriveInit}, {r,
v, γ, ϕ}, {t}, MaxSteps -> 5000000];
MotionVariables =
NDSolveSolutionDataComponent[Drive@"Variables",
"DependentVariables"];
NDSolveIterate[Drive, Tdrive]
DriveTrajectory = NDSolveProcessSolutions[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;
NDSolveProcessSolutions[Drive, "Forward"];
{Coast} =
NDSolveProcessEquations[{CoastMotion,
ReleaseHold[CoastInit] /. t -> 0}, {r, v, γ, ϕ}, {t},
MaxSteps -> 5000000];
NDSolveIterate[Coast, Tcoast]
CoastTrajectory = NDSolveProcessSolutions[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]];
`