Message Boards Message Boards

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

Posted 5 years ago

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.

POSTED BY: Meir Lewkowicz
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.

POSTED BY: Gianluca Gorni

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 BY: Tim Laska
Posted 5 years ago

Thanks Gianluca. It works.

POSTED BY: Meir Lewkowicz

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 BY: Tim Laska

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

POSTED BY: Gianluca Gorni
Posted 5 years 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.

POSTED BY: Meir Lewkowicz
Posted 5 years ago

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

POSTED BY: Meir Lewkowicz
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