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

Posted 6 months ago
1029 Views
|
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. Attachments:
7 Replies
Sort By:
Posted 6 months ago
 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 6 months ago
 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.
Posted 6 months ago
 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 *)