# Using perturbation theory to solve classical anharmonic oscillations

A point mass $m$ is allowed to move in the $x-y$ plane.

Given $k$, $a$ to be the spring constant and springs' natural length respectively.

The springs are parallel to the $x$ axis.

The Lagrangian of the system is given by:

$$L= \frac{m}{2} (\dot{x}^2+\dot{y}^2) - \frac{k}{2}\left( a\sqrt{\left( \frac{x}{a} - 1\right)^2 + \frac{y^2}{a^2}} -a\right)^2 - \frac{k}{2}\left( a\sqrt{\left( \frac{x}{a} + 1\right)^2 + \frac{y^2}{a^2}} -a\right)^2$$

Assuming small oscillation approximation for the above where $x \ll a$ and $y \ll a$

Using $1+\frac{x}{2}-\frac{x^2}{8}+\frac{x^3}{16}-\frac{5 x^4}{128}+O\left(x^5\right)$

We get

$$L \approx \frac{m}{2} (\dot{x}^2+\dot{y}^2)+ \frac{k}{a^2}x^2y^2 - kx^2-\frac{k}{4a^2} y^4$$

How do I solve this in order to plot $x$ vs $y$ motions (Lissajous figures), understanding the motion, analyze the equations of motion etc.

An interesting observation in this problem is that if we start where both springs are extended by $x_0$ at equilibrium, namely the springs are placed in $(-a-x_0, 0)$ and $(+a+x_0, 0)$a we get an harmonic motion in both $x$ and $y$ directions for small $(x,y)$ displacements the Lagrangian

$$L= \frac{m}{2} (\dot{x}^2+\dot{y}^2) - \frac{k}{2}\left( a\sqrt{\left( \frac{x+x_0}{a} - 1\right)^2 + \frac{y^2}{a^2}} -a\right)^2 - \frac{k}{2}\left( a\sqrt{\left( \frac{x+x_0}{a} + 1\right)^2 + \frac{y^2}{a^2}} -a\right)^2$$

is approximated as:

$$L_{oscillations} \approx \frac{m}{2} \left( \dot{x}^2+\dot{y}^2\right) - kx_0^2 - \frac{1}{x_0 + a} \left((k x_0 + a k) x^2 + k x_0 y^2\right)$$

There will almost certainly not be a closed-form solution for these equations. You could try going to higher order in the perturbation series, i.e., assume that our solutions can be expressed as a smooth power series in a small parameter $\epsilon$: $$x(t) = \epsilon x^{(1)} + \epsilon^2 x^{(2)} + \dots \\ y(t) = \epsilon y^{(1)} + \epsilon^2 y^{(2)} + \dots$$ Plugging these equations into the Euler-Lagrange equations then yields a set of equations that could in principle be solved order by order in $\epsilon$. (Conventional "small oscillation" theory just corresponds to doing this to $\mathcal{O}(\epsilon^1)$.)

The problem here is that if you do this for this system, the equation for $y^{(1)}$ will end up being $\ddot{y}^{(1)} = 0$, and so $y^{(1)} = At+B$. (This arises from the fact that there is no $y^2$ term in the expanded potential.) This is true enough for very short periods of time, but we know that in the long run the ball will not head off to $\pm \infty$. So this technique is of limited utility here; in particular, since $y^{(1)}$ will eventually get "large", the power series ansatz I wrote above will not necessarily converge to the "true" solution. (This technique is more useful when you don't have these so-called secular perturbations—i.e., parts of the solution that grow in time rather than oscillating–gumming up the works.)

Honestly, I would just use your favorite piece of software to get numerical approximations to some representative solutions, and plot them. Since you're going to be using numerical methods anyway, you might as well use the full potential rather than the approximation you have above. I suspect that the trajectory will in general not form a nice closed Lissajous-esque figure, but will rather fill out some particular region of space. Here, for example, is the trajectory for $x(0) = y(0) = 0, \dot{x}(0) = 1, \dot{y}(0) = 1/2$ (in units where $k = a = m = 1$): The blue curve is the actual trajectory for $0 \leq t \leq 200$; the red curve is the region of the plane that is allowed by energy conservation ($U(x,y) \leq E = \frac{1}{2}(\dot{x}_0^2 + \dot{y}_0^2)$.) I had suspected that the trajectory would fill out the red curve, but it does not appear to do so even when I run the simulation for fairly long intervals. I find this interesting, but have no clue how one would go about proving it.

Per the OP's request, the Mathematica code I used is below:

U[x_, y_] = (Sqrt[(x - 1)^2 + y^2] - 1)^2 + (Sqrt[(x + 1)^2 + y^2] - 1)^2
vx0 = 1; vy0 = 1/2;
eqns = {x''[t] == -D[U[a, b], a], y''[t] == -D[U[a, b], b],
x == y == 0, x' == vx0, y' == vy0} /. {a -> x[t],
b -> y[t]};
tmax = 200;
NDSolve[eqns, {x[t], y[t]}, {t, 0, tmax}];
trajectory =  ParametricPlot[{x[t], y[t]} /. First[%], {t, 0, tmax},  PlotPoints -> 200];
enbound = ContourPlot[U[x, y] == 1/2 (vx0^2 + vy0^2), {x, -3, 3}, {y, -3, 3}, ContourStyle -> Red];
Show[trajectory, enbound, PlotRange -> All]

• $y^{(n)}$ here is not a derivative. Rather, we are assuming that $y(t)$ can be expressed as a power series in some small parameter $\epsilon$. The coefficients of $\epsilon^n$ in this series are the functions $y^{(n)}(t)$. You might want to read up on perturbation theory if you've never done these kinds of higher-order calculations before. – Michael Seifert Oct 20 '16 at 14:43
• Michael, thank you for your explanation all is delightful. Please consider the intro on the perturbation ansatz and expansion. – 0x90 Oct 21 '16 at 17:21