1

I have this Diff Eq:

$$2g''+g'^{2}+ag+a-b=0 $$

We can manipulate this into an equivalent relation:

$$g'^{2}=ce^{-g}-ag+b $$

...which makes c a conserved quantity:

$$c=\left(g'^{2}+ag-b\right)e^{g}$$

If $sign(a)=-sign(c) $ and $|c|<|b-a|e^{\frac{b-a}{a}} $ then g will be periodic. This can be verified by plotting $y=ce^{-g} $ and $y=ag-b $ and noting exactly two points of intersection at which $g'=0 $ in Eq(2)

So my question is, what's the best way to calculate the period of g?

Obviously I could iteratively calculate g(x) from initial conditions, but I doubt that's the best approach because of the accumulating numerical errors. Also, it's calculating a lot of intermediate information that I don't need.

Is there any clever way to numerically calculate the period directly?

Jerry Guern
  • 2,764
  • If you have a conserved quantity, then you can do symplectic numerics and thereby obtain the period without accumulated error. (There is still error of course, but it does not worsen from one period to the next.) – Ian Oct 04 '15 at 00:18
  • @Ian Thanks. Could you post a reference link? Googling that term got me to symplectic integrators but I don't know if that what you meant and didn't see how it would get me a period of a diff eq. – Jerry Guern Oct 04 '15 at 00:29
  • It's the same notion. It gets you the period, or rather an approximation of the period, because the numerical solution will also be periodic (but with a slightly different period). – Ian Oct 04 '15 at 00:54
  • First, convert to a system: introduce $q=g,p=g'$, so that $p'=\frac{1}{2} \left ( -p^2 - aq - a + b \right )$ and $q' = p$. Then attempt to see whether there exists $H(p,q)$ such that $p'=-\frac{\partial H}{\partial q}$ and $q'=\frac{\partial H}{\partial p}$. In this case it looks like this actually doesn't happen, which would mean that your conserved quantity is more complicated than it often is in classical mechanics. But for a periodic trajectory there must be some conserved quantity; but without knowing it you can't set up an appropriate conservative integrator. – Ian Oct 04 '15 at 14:35
  • OK, so $G(g,g')=e^g g'^2 + a g e^g - b e^g$ is constant. But it is not in Hamiltonian form in the "natural" coordinates $(g,g')$. Instead we should choose generalized coordinates in which to formulate Hamilton's equations. Taking $g'$ as a coordinate seems like a good idea: we get $\frac{\partial G}{\partial g'}=2 e^g g'$. So if $g'$ is a momentum coordinate, then the position coordinate has time derivative $2 e^g g'$, so the position coordinate is (up to a constant that we can choose to be zero) $2 e^g$. – Ian Oct 04 '15 at 15:23
  • So your generalized position is $2 e^g$ and your generalized momentum is $g'$. So your Hamiltonian in these coordinates is $H(p,q)=\frac{p^2 q}{2} + \frac{a q \ln(q)}{2} - \frac{bq}{2}$. So unless I've made an error somewhere (quite likely, as I didn't do this very carefully), you can solve your problem by using a symplectic integrator on the Hamilton's equations arising from this Hamiltonian, which read $q'=pq,p'=-\left ( \frac{p^2}{2} + \frac{a \ln(q) + a}{2} - \frac{b}{2} \right )$. – Ian Oct 04 '15 at 15:23
  • If you want to find your $g$ over time, you find $q$ and take $g=\ln(q/2)$. But you can find the period in the generalized coordinates. – Ian Oct 04 '15 at 15:31
  • I did make at least one error: $\ln(q)$ should be $\ln(q/2)$, which also changes one of the derivatives slightly. – Ian Oct 04 '15 at 15:52
  • No, that much is fine: $g'$ only appears in the $e^g g'^2$ term, which is where $p$ appears in my $H$. And no, $g=\ln(q/2)$, because my development took $q=2e^g,p=g'$. Perhaps it was unclear because I got $q$ by starting with $q'=2e^g g'$ and then integrating. – Ian Oct 04 '15 at 15:57
  • @Ian It looks like this whole thing is going to come out equivalent to just using my first two equations to eliminate g', getting a diff eq in just g, and using that to numerically calculate g(x) until it returns to initial conditions. – Jerry Guern Oct 04 '15 at 16:17
  • Technically, yes, except that in the form I've suggested, you have a well-developed theory for how to cook up a numerical method whose trajectories are truly periodic. Without this structure one would expect the numerical method to not have any fully conserved quantity, which would make it slightly aperiodic, which can be a problem. – Ian Oct 04 '15 at 16:18

1 Answers1

1

As already noted, this equation preserves the quantity $e^g(g'^2+ag-b)$. So it is worth trying to formulate the problem as Hamilton's equations for this Hamiltonian in some system of coordinates. The coordinates cannot be $(p,q)=(g',g)$ as they would often be, so we have to try something else. If we try to keep $p=g'$, we have $\frac{\partial H}{\partial p}=2e^g g'$ which must be the same as $q'$ for Hamilton's equations. Thus we can take $q=2e^g$. We can now reformulate the equation as the Hamilton's equations for the Hamiltonian

$$H(p,q)=\frac{q}{2} \left ( p^2+a\ln(q/2)-b \right ).$$

These equations are

$$q'=\frac{\partial H}{\partial p}=pq \\ p'=-\frac{\partial H}{\partial q}=-\frac{1}{2} \left ( p^2+a \ln(q/2)-b \right ) - \frac{a}{2}.$$

This will be periodic provided $q$ never reaches zero or $+\infty$ (These would correspond to $g$ blowing up in finite time.) One can estimate the period by using a symplectic integrator on these equations. The numerical solution from this integrator will preserve a perturbed Hamiltonian, and therefore should give a good estimate of the period when the original system is periodic.

Ian
  • 104,572
  • Thank you! One final question: Is it best to combine these into $p'' = -pp' -ap/2$ and numerically generate one full period of p(x)? Or is it better to just use your final two equations to iterate on p(x) and q(x) simultaneously? – Jerry Guern Oct 04 '15 at 16:50
  • Ideally you would maintain the Hamilton's equations structure. However, I have just noticed that most of the famous symplectic integrators are for separable Hamiltonians, i.e. those of the form $H(p,q)=F(p)+G(q)$ (or possibly with an explicit time dependence in either term). The Hamiltonian I wrote above is not separable. So you will either need to choose a system of coordinates in which the Hamiltonian becomes separable, or else use a symplectic integrator for non-separable Hamiltonians. – Ian Oct 04 '15 at 16:58
  • I'm a bit skeptical that the former is even possible; perhaps it would be beneficial to ask another question about the latter. – Ian Oct 04 '15 at 17:00
  • I suspect the former is always possible though often impractically difficult and messy. The latter seems pretty straightforward even if the Hamiltonian wasn't separable. Thanks for all of your help on this. – Jerry Guern Oct 04 '15 at 17:04
  • Feel like trying the bounty problem now? :-) http://math.stackexchange.com/questions/1455526/how-to-solve-an-integral-with-a-gaussian-mixture-denominator – Jerry Guern Oct 04 '15 at 17:05