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

Questions
, ,
1

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.

Query.nb (92.4 KB)
Solution Matlab.pdf (238 KB)

2

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.

3

Thanks Gianluca. It works.

4

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

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

6

Thanks Tim. I realized the inversed sign.

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

7

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 *)
8

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