Finding a specific solution to an ODE system

[michaeljones2711 (Mike Jones)]

Hey guys,

I’m in my last year of my physics degree, and part of my dissertation project involves mathematica, modelling the orbital behaviour of a spacecraft using Ion propulsion.

I have 2 coupled differential equations being solved, r’‘(t) and o’'(t), in an NdSolve function (o is pheta in my case, modelling the change in angular momentum). It solves this with no problem, and plots the results fine. I now need to find the value of t and o at which r=5.3. It seems simple enough on paper, but i’m not sure how to translate that into my code.

Attached is a simplified version of my code. Any help would be greatly appreciated.

Thanks a lot!

[seanc (Sean Clarke)]

You are looking for a value of t between -10 and 10 such that:

First[r[t] /. TotalThrust] == 5.3

The “First” is there to get rid of the list. See the previously linked article on how to use the output of NDSolve for more info.

Usually you can take equations like the above and use NSolve or FindRoot on them to find a value of t that solves them.

But, It doesn’t look like that value exists according to your graph! It does look like r[t] goes up to 5.3. NDSolve is set in your code so that it runs from t=-10 to t=10. If I increase that boundary I still don’t see r[t] going anywhere near 5.3.

Here’s what it looks like when I run NDSolve from -20 to 20 and plot it:

Why do you believe that graph should go up to 5.3? Around what value of t do you think this would happen?

[michaeljones2711 (Mike Jones)]

Hello again!

Sorry about the disagreement with the graph, in my currently used model, different initial conditions are set such that r does eventually reach the desired value.

What is the syntax of the Nsolve and Findroot when solving a specific equation like First[r[t] /. TotalThrust] == 5.3?

I’ve attempted with NSolve [First[r[t] /. TotalThrust] == 5.3, t] and it returns a value of 16.22 which is wrong when looking at my graph, and returns the error : Inverse functions are being used by NSolve, so some solutions may not be found; use Reduce for complete solution information. >>

Thanks!

[michaeljones2711 (Mike Jones)]

Well actually, scratch that, 15.22 is a correct solution, however it is not the first solution on the graph, by which i mean

Is there any way to specify the first solution at around t=7

[seanc (Sean Clarke)]

Is there any way to specify the first solution at around t=7

Use FindRoot. FindRoot requires that you give it a starting value to guess from. So you can try something like

FindRoot[First[r[t] /. TotalThrust] == 5.3, {x,7}]

[michaeljones2711 (Mike Jones)]

Trying that, it returns the error FindRoot::nlnum: The function value {-5.3+<<1>>[t]} is not a list of numbers with dimensions {1} at {x} = {7.}. >>, i looked up the syntax and help documents but i couldn’t get it to give me anything except that error. I assume its to do with the First function as it is treating it like a list?

[dlkeith (David Keith)]

It would be easier to respond if you would post complete code or a new notebook. However, if you are not interested in the values for r after it first passes a limit, you could add WhenEvent to the NDSolve equations. In the code below I used your intial notebook and set rLimit to 2.5. WhenEvent stops integration when r hits 2.5, and records the time in tStop. So r is then r[tStop]. (Note that, for some reason, tStop needs to be defined prior to executing NDSolve or the value will not be saved.)

In[1]:= L = 2 Pi;
M = 1;
x = 35 Degree;

In[18]:= rLimit = 2.5; tStop = 0;

In[19]:= TotalThrust = 
  NDSolve[{r''[t] == - (4*Pi^2 / r[t]^2) + ((o'[t]^2)*r[t]) + 
      HeavisideTheta[5 - t]*0.51*Cos[x], 
    o''[t] == (-2 r'[t]*o'[t])/r[t] + 
      HeavisideTheta[5 - t]*0.51*Sin[x], o[0] == 0, o'[0] == 2*Pi, 
    r[0] == 1, r'[0] == 0, 
    WhenEvent[r[t] >= rLimit, {"StopIntegration", tStop = t}]}, {o [
     t], r[t]}, {t, -10, 10}];

In[21]:= tStop

Out[21]= 5.27189

[michaeljones2711 (Mike Jones)]

That works really well, but i prefer how dynamic the FindRoot method is, as i will be looking into the values after the limit. Any idea’s about the error it returned when i tried to use it?

[dlkeith (David Keith)]

Hi Mike,

Let’s not stop the integration – we’ll just accumulate a list of the positive crossings:

In[157]:= L = 2 Pi;
M = 1;
x = 35 Degree;

In[160]:= rLimit = 2.5; tList = {};

In[161]:= (* notice I used NDSolveValue to get o[t] and r[t] *)
{fo, fr} = 
  NDSolveValue[{r''[t] == -(4*Pi^2/r[t]^2) + ((o'[t]^2)*r[t]) + 
      HeavisideTheta[5 - t]*0.51*Cos[x], 
    o''[t] == (-2 r'[t]*o'[t])/r[t] + 
      HeavisideTheta[5 - t]*0.51*Sin[x], o[0] == 0, o'[0] == 2*Pi, 
    r[0] == 1, r'[0] == 0, 
    WhenEvent[r[t] >= rLimit, AppendTo[tList, {t, r[t]}]]}, {o[t], 
    r[t]}, {t, -10, 10}];

In[162]:= (* here are the positive crossings *)
tList

Out[162]= {{5.27189, 2.5}, {8.97694, 2.5}}

In[163]:= (* The first is of course at t = *)
tList[[1, 1]]

Out[163]= 5.27189

In[164]:= (* We can see them on the plot *)
Plot[fr, {t, -10, 10}, 
 Epilog -> {Red, PointSize[Large], Point /@ tList}]

Regarding my earlier comment, I think tStop, or in this case tList, needs to be defined before its use in WhenEvent because that provides it with global scope. Otherwise it has scope local to NDSolve and expires when NDSolve finishes.

Best Regards,
David

[michaeljones2711 (Mike Jones)]

This actually works fantastically! how do i manipulate the code to instead give me the value of o at our given rlimit?
This was my attempt, and it failed horribly!

rLimit = 5.2; tList = {}; oList = {};

{fo, fr} = 
  NDSolveValue[{r''[t] == - (4*Pi^2 / r[t]^2) + ((o'[t]^2)*r[t]) + 
      HeavisideTheta[5 - t]*(0.51*Cos[x]/((1 - A) + A (1 - B (t)))), 
    o''[t] == (-2 r'[t]*o'[t])/r[t] + 
      HeavisideTheta[
        5 - t]*(( 0.51*Sin[x])/(r[t]*((1 - A) + A (1 - B (t))))), 
    o[0] == 0, o'[0] == 2*Pi, r[0] == 1, r'[0] == 0, 
    WhenEvent[
     r[t] >= rLimit, {AppendTo[tList, {t, r[t]}], 
      AppendTo[oList, {o, r[o]}]}]}, {o[t], r[t]}, {t, -10, 20}];

Regards

Mike

[dlkeith (David Keith)]

In the code I posted, NDSolveValue creates interpolating functions for both o and r. Actually, because I asked it for o[t] and r[t], rather than o and r, the argument t is included in the output. If you execute the NDSolveValue without the semicolon you will see two interpolating functions, each expressed as a function of t. They are assigned to the variables fo and fr.

The first time at which rLimit is crossed is tList[[1,1]], so the value of o at that time is just fo/.t->tList[[1,1]] (Look up Rule and Replace.)

In[10]:= fo /. t -> tList[[1, 1]]

Out[10]= 20.0634

You can of course plot fo[t] just the way I plotted fr[t]. You could even put them both in the same plot.

Best,
David

[michaeljones2711 (Mike Jones)]

Thanks David, that worked perfectly! My final question is now that i have our value for t* at r=5.2, how do i find o’(t*) (the first derivative of o at our value of t)?

Cheers!
Mike

[michaeljones2711 (Mike Jones)]

Finding a specific solution to an ODE system