Message Boards Message Boards

0
|
7796 Views
|
12 Replies
|
1 Total Likes
View groups...
Share
Share this post:

Finding a specific solution to an ODE system

Posted 9 years ago

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!

Attachments:
POSTED BY: Mike Jones
12 Replies
Posted 9 years ago

Finding a specific solution to an ODE system

POSTED BY: Mike Jones
Posted 9 years ago

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

POSTED BY: Mike Jones
Posted 9 years ago

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

POSTED BY: David Keith
Posted 9 years ago

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

POSTED BY: Mike Jones
Posted 9 years ago

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?

POSTED BY: Mike Jones
Posted 9 years ago

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}]

enter image description here

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

POSTED BY: David Keith
Posted 9 years ago

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
POSTED BY: David Keith
Posted 9 years ago

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?

POSTED BY: Mike Jones

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}]
POSTED BY: Sean Clarke
Posted 9 years ago

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

enter image description here

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

POSTED BY: Mike Jones

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: enter image description here

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

POSTED BY: Sean Clarke
Posted 9 years ago

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!

POSTED BY: Mike Jones
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract