Effect of anharmonicity on an oscillator

Single harmonic oscillator has been ubiquitous in both classical as well as quantum mechanics. The harmonic oscillator is governed by the potential:

$U(x) = \dfrac{1}{2}kx^2 $

It is quite easy to see that the equations of motion for this case becomes:

$\dot(x) =p$

$\dot(p) =-x$

Any deviation in potential of the first equation causes the oscillator to forgo its "harmonicity" and thereby, enters anharmonic nature. Anharmonicity plays an important role in several toy models of thermal conduction like Fermi-Pasta-Ulam and $\phi^4$ chain. The focus of this short write-up is to show how the numerical solution changes when the system is no longer harmonic. For simplicity we assume that the anharmonic potential is of the form:

$U(x) = \dfrac{1}{2}x^2 + \dfrac{1}{3}x^3 + \dfrac{1}{4}x^4 +...+ \dfrac{1}{2n+1}x^{2n+1} $

 Note, we terminate at the index 2n+1 because otherwise the solutions diverge. The logic of such divergence is easy to understand: The force is always repulsive and further the oscillator moves, greater is the repulsive force. The MATLAB codes are attached below for the interested readers. The four resulting solutions are given below:

It is interesting to point out that the differential equations evidently become stiff under the effect of higher order nonlinear terms. The simplest way to account for it is to use the ode23s differential equation solver in MATLAB. The ode45 solver that has been used to plot these results utilize the adaptive Runge-Kutta algorithm (Runge-Kutta-Fehlberg), and can be easily implemented in a C / Fortran code. However, implementing the C / Fortran analogue of ode23s is not very trivial.

I will try to upload the results with ode23s solver for a better comparison of the effect of nonlinearity. Until then, you can use the ode45 solver codes attached below.

The MATLAB codes are given below:

main file
function file