Message Boards Message Boards


[✓] Why NDSolve give a run-away solution for a hamiltonian system?

Posted 6 months ago
7 Replies
4 Total Likes

I have a Hamiltonian system: r'=f(r,p), p'=g(r,p}. I ask for r(t), p(t). From energy conservation I know the limits for r. The solution for r(t) I get is correct for the beginning t's, but then runs off slowly to very large values. I have a solution from Matlab a colleague sent, which shows the correct oscillating behaviour in the correct limits for r(t). I checked my input, seems all fine. Any ideas? Thanks I attach the program and the solution from Matlab.

7 Replies

Check your Hamiltonian equations. If I use these:

{p'[t] == (Simplify[D[H[p, r], r], r + \[Epsilon] > 0] /. {r -> r[t], 
     p -> p[t]}),
 r'[t] == -(Simplify[D[H[p, r], p], r + \[Epsilon] > 0] /. {r -> r[t],
       p -> p[t]})}

I get the same plots as Matlab.

Nice elegant answer, but when I plotted your solution, the plot was coming out in reverse to Matlab. I needed to reverse the signs of the r and p equations. Do you think the signs should be reversed?

trq = NDSolve[{p'[
     t] == -(Simplify[D[H[p, r], r], 
        r + \[Epsilon] > 0] /. {r -> r[t], p -> p[t]}), 
   r'[t] == (Simplify[D[H[p, r], p], 
       r + \[Epsilon] > 0] /. {r -> r[t], p -> p[t]}), r[0] == 1.2, 
   p[0] == .1}, {r, p}, {t, 18}]
Plot[Evaluate[{r[t], p[t]} /. trq], {t, 0, 18}, PlotRange -> All, 
 AxesOrigin -> {0, 0}]
Posted 6 months ago

Thanks Gianluca. It works.

Yes, it's either that or the time-reversed.

Posted 6 months ago

Thanks Tim. I realized the inversed sign.

Frankly, I still do not understand why Gianluca's code works, and my 'straightforward' one, does not.

You are welcome Meir.

There appears to be a discrepancy between your momentum equation and Gianluca's. We can check the equivalence by comparing the Right Hand Sides (RHS) of the position and momentum equations your and Gianluca's approach. If they are equivalent, FullSimplify should evaluate to "True". As the following code demonstrates, this does not apply to the momentum equation.

v[r_] := Sqrt[-(Q^2/(eps + r)^2) + (2*M)/(eps + r)]/mp;
H[p_, r_] := 
  Sqrt[m^2 + p^2 + 
     p^4/mp^2] - (p*Sqrt[-(Q^2/(eps + r)^2) + (2*M)/(eps + r)])/mp;
meirposition = (2*p[t] + (4*p[t]^3)/mp^2)/(2*
      Sqrt[m^2 + p[t]^2 + p[t]^4/mp^2]) - 
   Sqrt[-(Q^2/(eps + r[t])^2) + (2*M)/(eps + r[t])]/mp;
gianlucaposition = (Simplify[D[H[p, r], p], 
     r + eps > 0] /. {r -> r[t], p -> p[t]});
meirmomentum = -((p[t]*(-(Q^2/(eps + r[t])^2) + M/(eps + r[t])))/(mp*
       r[t]*Sqrt[-(Q^2/(eps + r[t])^2) + (2*M)/(eps + r[t])]));
gianlucamomentum = -(Simplify[D[H[p, r], r], 
      r + eps > 0] /. {r -> r[t], p -> p[t]});
FullSimplify[meirposition == gianlucaposition, 
 r[t] + eps > 0] (* True *)
FullSimplify[meirmomentum == gianlucamomentum, 
 r[t] + eps > 
  0] (* (eps*p[t]*(eps*M-Q^2+M*r[t]))/(mp*r[t]*Sqrt[2*eps*M-Q^2+2*M*r[\
t]])\[Equal]0 *)
Posted 6 months ago

Thanks Tim. I found the the missing epsilon. Thanks to all offering their help. I appreciate it.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract