The Duffing oscillator

For the harmonic oscillator, the motion is almost magically simple. Confined to a quadratic potential $V(x)=ax^2$ , a displaced particle oscillates as a sine or cosine forever. The system possesses a single, intrinsic clock: its period is strictly independent of amplitude.

But if we add even the simplest nonlinear correction, $\displaystyle V(x)=ax^2+bx^4,$ the story changes in a deep way. The motion is still periodic. Energy is still conserved. The particle still oscillates between two turning points. But the waveform is no longer an ordinary cosine. It becomes a Jacobi elliptic function.

Consider the one-dimensional Hamiltonian

$\displaystyle H(x,p)=\frac{p^2}{2m}+ax^2+bx^4,\qquad a>0,\quad b>0.$

This potential is confining, even, and strictly convex, guaranteeing a stable equilibrium at the origin. The corresponding equation of motion defines the undamped, undriven hardening Duffing oscillator:

$\displaystyle \ddot{x}+\frac{2a}{m}x+\frac{4b}{m}x^3=0.$

Remarkably, despite the nonlinearity, this system remains exactly integrable. The trigonometric cosine of the harmonic oscillator is replaced by a Jacobi elliptic cosine.

Let $\displaystyle A>0$ denote the amplitude (the positive turning point). Then kinetic enrgy vanishes, $\displaystyle \dot{x}=0$ at $x=A,$ so the energy is $\displaystyle E=aA^2+bA^4.$ Energy conservation gives $\displaystyle E=\frac{m}{2}\dot{x}^2+ax^2+bx^4.$ Hence

$\displaystyle \dot{x}^2=\frac{2}{m}\left(aA^2+bA^4-ax^2-bx^4\right).$

Factor the quartic expression: $\displaystyle aA^2+bA^4-ax^2-bx^4=(A^2-x^2)(a+bA^2+bx^2).$ Thus

$\displaystyle \dot{x}^2=\frac{2}{m}(A^2-x^2)(a+bA^2+bx^2).$

This is the essential first integral. The time coordinate is therefore obtained from

$\displaystyle t=\sqrt{\frac{m}{2}}\int^x\frac{dq}{\sqrt{(A^2-q^2)(a+bA^2+bq^2)}}.$

The square root contains a quartic polynomial in $\displaystyle q$ . Thus the inverse function is elliptic.

We seek a solution of the form $\displaystyle x(t)=A~{\text{cn}}\left(\Omega(t-t_0),k\right).$ in terms of the Jacobi elliptic cosine function.

The useful identity is $\displaystyle \frac{d^2}{du^2}{\text{cn}}(u,k)=(2k^2-1){\text{cn}}(u,k)-2k^2{\text{cn}}^3(u,k).$

Let $\displaystyle u=\Omega(t-t_0).$ Then $\displaystyle x=A~{\text{cn}}(u,k),$ and so

$\displaystyle \ddot{x}=\Omega^2(2k^2-1)x-\frac{2k^2\Omega^2}{A^2}x^3.$

Matching this with $\ddot{x}=-\frac{2a}{m}x-\frac{4b}{m}x^3$ gives two algebraic equations: $\Omega^2(2k^2-1)=-\frac{2a}{m},$ and $\frac{2k^2\Omega^2}{A^2}=\frac{4b}{m}.$ From the second equation, $k^2\Omega^2=\frac{2bA^2}{m}.$ From the first equation, $\Omega^2-2k^2\Omega^2=\frac{2a}{m}.$ Substituting the previous relation gives $\Omega^2-\frac{4bA^2}{m}=\frac{2a}{m}.$ Therefore $\Omega^2=\frac{2(a+2bA^2)}{m}.$ Then $k^2=\frac{bA^2}{a+2bA^2}.$

Thus the exact solution is

$\displaystyle x(t)=A~{\text{cn}}\left(\sqrt{\frac{2(a+2bA^2)}{m}}(t-t_0),k\right)$

with $\displaystyle k^2=\frac{bA^2}{a+2bA^2}.$

The amplitude is related to energy by $E=aA^2+bA^4.$ Equivalently, $A^2=\frac{-a+\sqrt{a^2+4bE}}{2b}.$

The Jacobi function $\displaystyle {\text{cn}}(u,k)$ has real period $\displaystyle 4K(k),$ where

$\displaystyle K(k)=\int_0^{\pi/2}\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}$

is the complete elliptic integral of the first kind. Since $\displaystyle u=\Omega(t-t_0),$ the physical period is $\displaystyle T(A)=\frac{4K(k)}{\Omega}.$

Hence $\displaystyle T(A)=\frac{4K(k)}{\sqrt{2(a+2bA^2)/m}},\qquad k^2=\frac{bA^2}{a+2bA^2}.$ The nonlinear frequency is $\displaystyle \omega(A)=\frac{2\pi}{T(A)}=\frac{\pi\Omega}{2K(k)}.$ Thus

$\displaystyle \omega(A)=\frac{\pi}{2K(k)}\sqrt{\frac{2(a+2bA^2)}{m}}.$

This is the main physical observable. Unlike the harmonic oscillator, the frequency is not constant. It depends on amplitude, hence on energy.

The amplitude is often less intrinsic than the energy. From $\displaystyle E=aA^2+bA^4,$ we get $\displaystyle A^2=\frac{-a+\sqrt{a^2+4bE}}{2b}.$ Then $\displaystyle a+2bA^2=\sqrt{a^2+4bE}.$ Therefore $\displaystyle \Omega(E)=\sqrt{\frac{2\sqrt{a^2+4bE}}{m}}.$ Also, $\displaystyle bA^2=\frac{-a+\sqrt{a^2+4bE}}{2},$ so

$\displaystyle k^2(E)=\frac{bA^2}{a+2bA^2}=\frac{-a+\sqrt{a^2+4bE}}{2\sqrt{a^2+4bE}}.$

Thus $\displaystyle k^2(E)=\frac{1}{2}\left(1-\frac{a}{\sqrt{a^2+4bE}}\right).$

The exact energy-frequency relation is therefore

$\displaystyle \omega(E)=\frac{\pi}{2K(k(E))}\sqrt{\frac{2\sqrt{a^2+4bE}}{m}}.$

This formula exhibits the interpolation between the quadratic and quartic regimes.

Small Amplitudes

At small amplitude, $\displaystyle bA^2\ll a.$

Then $\displaystyle k^2=\frac{bA^2}{a+2bA^2}=\frac{bA^2}{a}+O(A^4).$

Also, $\displaystyle \Omega=\sqrt{\frac{2a}{m}}\left(1+\frac{bA^2}{a}+O(A^4)\right).$

The elliptic integral has expansion

$\displaystyle K(k)=\frac{\pi}{2}\left(1+\frac{k^2}{4}+\frac{9k^4}{64}+O(k^6)\right).$

Therefore $\displaystyle \frac{\pi}{2K(k)}=1-\frac{k^2}{4}+O(k^4).$

Hence $\displaystyle \omega(A)=\sqrt{\frac{2a}{m}}\left(1+\frac{bA^2}{a}\right)\left(1-\frac{bA^2}{4a}\right)+O(A^4).$

Thus $\displaystyle \omega(A)=\sqrt{\frac{2a}{m}}\left(1+\frac{3bA^2}{4a}+O(A^4)\right).$

The leading correction is positive. The oscillator hardens.

In terms of the usual Duffing notation, $\displaystyle \ddot{x}+\omega_0^2x+\alpha x^3=0,$ we have $\displaystyle \omega_0^2=\frac{2a}{m},\quad \alpha=\frac{4b}{m}.$ Then the small-amplitude frequency shift becomes

$\displaystyle \omega(A)=\omega_0\left(1+\frac{3\alpha A^2}{8\omega_0^2}+O(A^4)\right),$

which is the standard Duffing result.

Harmonic Oscillator Limit

Let $\displaystyle b\to 0.$ Then $\displaystyle k\to 0,$ and $\displaystyle {\text{cn}}(u,0)=\cos u.$ Also, $\displaystyle \Omega\to \sqrt{\frac{2a}{m}}.$

Thus the solution reduces to

$\displaystyle x(t)=A\cos\left(\sqrt{\frac{2a}{m}}(t-t_0)\right).$

The period becomes $\displaystyle T_0=2\pi\sqrt{\frac{m}{2a}}.$ This confirms that the elliptic solution is not a different type of motion. It is the exact nonlinear continuation of the harmonic oscillator.

Quartic Limit

Now set $\displaystyle a=0.$ Then $\displaystyle V(x)=bx^4,$ and $\displaystyle E=bA^4.$ The modulus becomes

$\displaystyle k^2=\frac{1}{2}.$

The frequency scale becomes $\displaystyle \Omega=2A\sqrt{\frac{b}{m}}.$

Hence $\displaystyle x(t)=A~{\text{cn}}\left(2A\sqrt{\frac{b}{m}}(t-t_0),\frac{1}{\sqrt{2}}\right).$

The period is $\displaystyle T=\frac{4K(1/\sqrt{2})}{2A\sqrt{b/m}}.$ Therefore $\displaystyle T=\frac{2K(1/\sqrt{2})}{A}\sqrt{\frac{m}{b}}.$

Since $\displaystyle A=\left(\frac{E}{b}\right)^{1/4},$ we have $\displaystyle T\sim E^{-1/4},\qquad \omega\sim E^{1/4}.$

This scaling can also be obtained without elliptic functions. In the pure quartic problem, dimensional analysis gives $\displaystyle [E]=[bA^4],$ so $\displaystyle A\sim \left(\frac{E}{b}\right)^{1/4}.$ The acceleration scale is $\displaystyle \ddot{x}\sim \frac{bA^3}{m},$ so the time scale satisfies $\displaystyle \frac{A}{T^2}\sim \frac{bA^3}{m}.$

Thus $\displaystyle T\sim \sqrt{\frac{m}{bA^2}}\sim \sqrt{\frac{m}{b}}\frac{1}{A}\sim m^{1/2}b^{-1/4}E^{-1/4}.$ The elliptic calculation supplies the exact numerical constant $\displaystyle 2K(1/\sqrt{2}).$

There is only one essential dimensionless nonlinearity parameter.

Let $\displaystyle \omega_0=\sqrt{\frac{2a}{m}},\quad \tau=\omega_0 t,\quad x=Ay.$

The equation $\displaystyle \ddot{x}+\frac{2a}{m}x+\frac{4b}{m}x^3=0$ becomes

$\displaystyle \frac{d^2y}{d\tau^2}+y+\lambda y^3=0,$

where $\displaystyle \lambda=\frac{2bA^2}{a}.$

Thus the shape of the motion is governed by $\displaystyle \lambda.$

In these variables, $\displaystyle k^2=\frac{bA^2}{a+2bA^2}.$ Since $\displaystyle bA^2=\frac{\lambda a}{2},$ we obtain

$\displaystyle k^2=\frac{\lambda}{2(1+\lambda)}.$

Thus $\displaystyle \lambda\ll 1\quad\Longrightarrow\quad k^2\ll 1,$ and the motion is nearly sinusoidal.

Meanwhile, $\displaystyle \lambda\to\infty\quad\Longrightarrow\quad k^2\to\frac{1}{2},$ which is the pure quartic limit.

This is a useful structural fact: for the stable positive quartic oscillator, the elliptic modulus never approaches $\displaystyle 1$ . It ranges only over $\displaystyle 0\leq k^2<\frac{1}{2}.$ So the elliptic behavior is significant, but it is not the separatrix-type behavior associated with $\displaystyle k\to 1$ .

Phase Space

Energy conservation gives the phase curve $\displaystyle \frac{m}{2}\dot{x}^2+ax^2+bx^4=E.$ For $\displaystyle b=0$ , this is an ellipse in the $\displaystyle (x,\dot{x})$ plane: $\displaystyle \frac{m}{2}\dot{x}^2+ax^2=E.$

For $\displaystyle b>0$ , the orbit is no longer an ellipse. The turning points are still $\displaystyle x=\pm A,$ but the velocity is

$\displaystyle \dot{x}=\pm \sqrt{\frac{2}{m}(A^2-x^2)(a+bA^2+bx^2)}.$

The extra factor $\displaystyle a+bA^2+bx^2$ weights the velocity differently across the orbit. Near the origin, $\displaystyle \dot{x}^2\big|_{x=0}=\frac{2}{m}A^2(a+bA^2).$ For large amplitude, this is dominated by $\displaystyle \dot{x}^2\big|_{x=0}\sim \frac{2b}{m}A^4.$ So the maximum speed scales as $\displaystyle v_{\max}\sim A^2\sqrt{\frac{2b}{m}}$ in the quartic regime, rather than linearly in $\displaystyle A$ as in the harmonic regime. The phase-space orbit is therefore compressed in time: large-amplitude trajectories are traversed faster.

Action Variables

For a one-dimensional bound system, the action is

$\displaystyle J(E)=\frac{1}{2\pi}\oint pdx.$

Here $\displaystyle p=m\dot{x},$ so $\displaystyle p(x)=\sqrt{2m(E-ax^2-bx^4)}.$

Therefore $\displaystyle J(E)=\frac{2}{\pi}\int_0^A \sqrt{2m(E-ax^2-bx^4)}~dx.$

The frequency is related to the action by $\displaystyle \omega(E)=\frac{dE}{dJ}.$ Equivalently, $\displaystyle \frac{dJ}{dE}=\frac{1}{\omega(E)}.$

This gives another conceptual explanation for why elliptic integrals appear: the action integral contains the square root of a quartic polynomial.

The period can also be written as $\displaystyle T(E)=\oint \frac{dx}{\dot{x}}.$

Using the turning point symmetry, $\displaystyle T(E)=4\int_0^A \frac{dx}{\sqrt{\frac{2}{m}(E-ax^2-bx^4)}}.$

Substituting the factorization,

$\displaystyle T(E)=4\sqrt{\frac{m}{2}}\int_0^A\frac{dx}{\sqrt{(A^2-x^2)(a+bA^2+bx^2)}}.$

Let $\displaystyle x=A\sin\theta.$ Then $\displaystyle dx=A\cos\theta,d\theta,$ and $\displaystyle A^2-x^2=A^2\cos^2\theta.$ Thus

$\displaystyle T(E)=4\sqrt{\frac{m}{2}}\int_0^{\pi/2}\frac{d\theta}{\sqrt{a+bA^2+bA^2\sin^2\theta}}.$

Factor $\displaystyle a+2bA^2$ from the denominator: $\displaystyle a+bA^2+bA^2\sin^2\theta=(a+2bA^2)\left(1-\frac{bA^2}{a+2bA^2}\cos^2\theta\right).$

Hence $\displaystyle T(E)=4\sqrt{\frac{m}{2(a+2bA^2)}}\int_0^{\pi/2}\frac{d\theta}{\sqrt{1-k^2\cos^2\theta}},$

where $\displaystyle k^2=\frac{bA^2}{a+2bA^2}.$

Since replacing $\displaystyle \cos^2\theta$ by $\displaystyle \sin^2\theta$ is just the change $\displaystyle \theta\mapsto \pi/2-\theta$ , the integral is $\displaystyle K(k)$ . Therefore