This section as well others related to it gives an introduction to nonlinear motion that is observed in many applications. Linear model is valid only in the very small amplitude approximations that a one-dimensional system trapped in the neighborhood of a point of stable equlibrium, which leads to a simple harmonic oscillator. Injection of energy into such a system increases the amplitude of its oscillations, causing the particle to begin to explore the completely different phenomenon that we observe for a simple idealized spring. Therefore, this section can be considered as an introduction to the theory of nonlinear oscillations and related topics.

To make exposition relatively simple and understandable at undergraduate level, we mostly work with some examples utilizing software. The reader should anticipate that numerical methods will play more prominent role in our qualitative analysis.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part III of the course APMA0340
Introduction to Linear Algebra with Mathematica

Spring-mass system

    Consider a simple system consisting of a spring and attached to it a single mass m. Suppose that the mass is put into motion my stretching it and applying some impulse. Assuming its one dimensional motion, we choose rectangular coordinates so that the origin x = 0 corresponds to its stable equilibrium position (where the object would stay at rest if it was released at rest). Once released, the restoring force causes the mass to move back toward its stable equilibrium position, where the net force on it is zero. However, by the time the mass gets there, it gains momentum and continues to move through, producing the opposite deformation. It is then forced to the opposite through equilibrium, and the process is repeated until dissipative forces (e.g., friction) dampen the motion.

Due to various sources of nonlinearities, micro/nano-electro-mechanical-system (MEMS/NEMS) resonators present highly nonlinear behaviors including softening- or hardening-type frequency responses, bistability, chaos, etc. The general Duffing equation with quadratic and cubic nonlinearities serves as a characterizing model for a wide class of MEMS/NEMS resonators as well as lots of other engineering and physical systems.

The simplest oscillations occur when the restoring force is directly proportional to displacement. The name that was given to this relationship between force and displacement is Hooke’s law: \( f(x) = k\, x . \) Here, f is the restoring force, x(t) is the displacement from equilibrium or deformation at time t, and k is a constant related to the difficulty in deforming the system (often called the spring constant or force constant), with units of newtons per meter (N/m). For example, k is directly related to Young’s modulus when we stretch a string. Then Newton's second law yields

\[ m\,\ddot{x} = -f(x) , \]
where the minus sign indicates that the restoring force acts in the direction opposite to the displacement, and dot stands for the derivative with respect to time variable t. However, the restoring force is rearly Hookean, and its determination could be very complicated. For simplicity, we start with undampted case when f(x) does not depend on the velocity.

Multiplying the main equation \( m\,\ddot{x} + f(x) = 0 \) by the velocity \( \dot{x} = {\text d}x/{\text d}t , \) we obtain

\[ \frac{m}{2}\,\ddot{x}\,\dot{x} + \dot{x}\, f(x) =0 . \]
Upon integration, we get
\[ \frac{m}{2}\,\dot{x}^2 + \Pi (x) = E, \qquad \Pi (x) = \int_0^t \dot{x}\, f(x) \,{\text d}t , \]
where E is a constant, representing the total energy of the system. We immediately recognize that Π is the potential energy and the first term \( \mbox{K} = \frac{m}{2}\,\dot{x}^2 \) represents the kinetic energy. The potential energy function of a one-dimensional oscillator about its stable equilibrium point x = 0 can usually be expanded in a Taylor series:
\[ \Pi (x) = \Pi_0 + \Pi^{(1)}_0 x + \frac{1}{2!}\,\Pi^{(2)}_0 x^2 + \frac{1}{3!}\,\Pi^{(3)}_0 x^3 + \frac{1}{4!}\,\Pi^{(4)}_0 x^4 + \cdots , \]
where Π0 is the value of Π(x) at the equilibrium point where Π(x) = 0. This constant term Π0 can be set equal to zero without loss of generality because the potential energy is only defined to within an additive constant. Furthermore, the coefficient of the first order term vanishes because at the point of stable equilibrium, the potential energy function has local minimum. Consequently, the power series expansion starts with a leading term of even power; so at least we have
\[ \Pi (x) = \frac{1}{2!}\,\Pi^{(2)}_0 x^2 + \frac{1}{3!}\,\Pi^{(3)}_0 x^3 + \frac{1}{4!}\,\Pi^{(4)}_0 x^4 + \cdots , \]
and the behavior of the oscillator depends on the values of the remaining coefficients in the series. Since the potential energy function has a local minimum at x = 0, the coefficient of the second order term cannot be negative: Π0(2) > 0. Any spring--mass model has the property that more energy injected into the system leads to increased amplitude.

The above expansion cannot be always true because the potential energy can be represented by a piecewise continuous function for which the Taylor series expansion is not valid. Moreover, even when such an expansion is possible, there are cases where the coefficient of the quadratic term vanishes and the series starts with a term higher than the second order.

Let us return to the string-mass problem when the potential energy retains the Taylor series expansion. Then such system can eb modeled by the following initial value problem:

\[ m\, \ddot{x} + k\,x + a_2 x^2 + a_3 x^3 + \cdots + a_n x^n =0, \qquad x(0) = x_0 , \quad \dot{x}(0) = v_0 . \]
The coresponding total energy is
\[ \frac{m}{2}\, \dot{x}^2 + \frac{k}{2}\, x^2 + \frac{a_2}{3}\, x^3 + \cdots + \frac{a_n}{n+1}\, x^{n+1} = E \equiv \frac{m}{2}\, v_0^2 + \frac{k}{2}\, x_0^2 + \frac{a_2}{3}\, x_0^3 + \cdots + \frac{a_n}{n+1}\, x^{n+1}_0 , \]
where E is a constant. Setting \( y = \dot{x} , \) the contours of the energy function
\[ E(x,y) = \frac{m}{2}\, v_0^2 + \frac{k}{2}\, x_0^2 + \frac{a_2}{3}\, x_0^3 + \cdots + \frac{a_n}{n+1}\, x^{n+1}_0 , \]
represent trajectories in the phase plane. More precisely, a trajectory beginning at (x0, v0) has an energy level of
\[ E = \frac{m}{2}\, v_0^2 + \frac{k}{2}\, x_0^2 + \frac{a_2}{3}\, x_0^3 + \cdots + \frac{a_n}{n+1}\, x^{n+1}_0 , \]
and the trajectory can be viewed as the contour
\[ E(x,y) = E . \]
As is usual, to study the differential equation modelled spring-mass system, we normalize the mass m = 1 and the leading coefficient k = 1 by setting
\[ \tau = t\,\sqrt{k/m} \]
and returning back to the time variable t instead of τ. Next we rewrite the given edifferential equation of the second order as the system of first order equations:
\[ \left. \begin{array}{l} \dot{x} = y , \\ \dot{y} = -k\,x- a_2 x^2 - \cdots - a_n x^n . \end{array} \right\} \]
The equilibrium solutions of the above system of differential equations are precisely the critical values of the energy function E(x, y):
\[ \frac{\partial E}{\partial x} (x,y) =0, \qquad \frac{\partial E}{\partial y} (x,y) =0 . \]
It is easy to see that
\[ \frac{\partial^2 E}{\partial x^2} = k + 2\,a_2 x + \cdots + n\,a_n x^{n-1} , \quad \frac{\partial^2 E}{\partial x\,\partial y} = 0, \quad \frac{\partial^2 E}{\partial y^2} = 1, \]
so that the expression of second derivatives becomes
\[ D = \frac{\partial^2 E}{\partial x^2}\, \frac{\partial^2 E}{\partial y^2} - \frac{\partial^2 E}{\partial x\,\partial y} = \frac{\partial^2 E}{\partial x^2} \left( x_0, y_0 \right) , \]
which we denote by \( D = \frac{\partial^2 E}{\partial x^2} \left( x_0, y_0 \right) . \) Thus, the classification of the critical value depends solely upon the sigh of D:
if D > 0 and Exx = D > 0, then E has a local minimum;
if D < 0, then E has a a saddle point at (x0, y0).
If D = 0, the test is inconclusive.

The above system of differential equations always thas the critical point at the origin. Since its Jacobian

\[ {\bf J} (0,0) = \begin{bmatrix} 0&1 \\ -k& 0 \end{bmatrix} \]
has two pure imaginary eigenvalues
\[ \lambda = \pm \sqrt{-k} = \pm {\bf j}\,\sqrt{k} , \]
the linerized system has a center at the origin. It is also a center for the nonlinear system because the origin is a local minimum for the energy function.

We can reduce the spring-mass differential equation \( \ddot{x} + x + a_2 x^2 + a_3 x^3 + \cdots + a_n x^n =0 \) to a separable differential equation of the first order using substitution: \( \dot{x} = 1/y . \) This yields

\[ y' = y^3 \left( x + a_2 x^2 + a_3 x^3 + \cdots + a_n x^n \right) \qquad \Longrightarrow \qquad \frac{{\text d}y}{y^3} = {\text d} x \left( x + a_2 x^2 + a_3 x^3 + \cdots + a_n x^n \right) . \]
Integrating both sides, we obtain
\[ \frac{1}{y} = \dot{x} = \pm \sqrt{C - \frac{x^2}{2} + \frac{a_2}{3}\, x^3 - \cdots - \frac{a_n}{n+1}\, x^{n+1}} , \]
where C is an arbitrary constant. FRom energy conservation \( \frac{1}{2}\,m\dot{x}^2 + \Pi (x) = E , \) it follows that the time dt required for the particle to move from x to x + dx can be described as
\[ {\text d}t = \frac{{\text d}x}{\sqrt{\frac{2}{m} \left[ E - \Pi (x) \right]}} \qquad\Longrightarrow \qquad \mbox{time of transition from }x_0 \ \mbox{ to } x = \int_{x_0}^x \frac{{\text d}y}{\sqrt{\frac{2}{m} \left[ E - \Pi (y) \right]}} . \]

Example: We consider the pendulum equation and its two approximations:

\begin{align*} \Pi (\theta ) &= mg\ell \left[ 1 - \cos \theta \right] : \qquad &\mbox{exact simple pendulum} \\ &\approx mg\ell \left[ \frac{1}{2}\,\theta^2 \right] : \qquad &\mbox{(harmonic) approximation} \\ &\approx mg\ell \left[ \frac{1}{2}\,\theta^2 - \frac{1}{24}\,\theta^4 \right] : \qquad &\mbox{harmonic approximation + quartic correction} \end{align*}
Using Mathematica, we plot some trajectories. We start with cubic term
\[ E(p,x) = p^2 + x^2 - \frac{1}{3}\, x^3 , \qquad p = \dot{x} , \]
and then go to quadric term:
\[ E(p,x) = p^2 + x^2 - \frac{1}{4}\, x^4 , \qquad p = \dot{x} , \]
q4 = ContourPlot[{p^2 + x^2 - x^3/3 == 4}, {x, -15, 10}, {p, -6, 6}, AspectRatio -> Automatic, Axes -> True, Frame -> False, ClippingStyle -> Automatic, ContourLabels -> True]
q3 = ContourPlot[{p^2 + x^2 - x^3/3 == 3}, {x, -15, 10}, {p, -6, 6}, AspectRatio -> Automatic, Axes -> True, Frame -> False, ClippingStyle -> Automatic, ContourLabels -> True]
q2 = ContourPlot[{p^2 + x^2 - x^3/3 == 2}, {x, -15, 10}, {p, -6, 6}, AspectRatio -> Automatic, Axes -> True, Frame -> False, ClippingStyle -> Automatic, ContourLabels -> True]
q1 = ContourPlot[{p^2 + x^2 - x^3/3 == 1}, {x, -15, 10}, {p, -6, 6}, AspectRatio -> Automatic, Axes -> True, Frame -> False, ClippingStyle -> Automatic, ContourLabels -> True]
q43 = ContourPlot[{p^2 + x^2 - x^3/3 == 4/3}, {x, -15, 10}, {p, -6, 6}, AspectRatio -> Automatic, Axes -> True, Frame -> False, ClippingStyle -> Automatic, ContourLabels -> True, ColorFunction -> Hue]
q05 = ContourPlot[{p^2 + x^2 - x^3/3 == 0.5}, {x, -15, 10}, {p, -6, 6}, AspectRatio -> Automatic, Axes -> True, Frame -> False, ClippingStyle -> Automatic, ContourLabels -> True]
Show[q05, q1, q43, q2, q3, q4]
Phase portrait of cubic potential energy


  1. Fay, T.H. and Joubert, S.V., Energy and nonsymmetric nonlinear spring, International Journal of Mathematical Education in Science and Technology, 1999, Vol. 30, No. 6, pp. 889--902;
  2. Kovacic, I. and Brennan, M.J., The Duffing equation, nonlinear oscillators and their behavior 2011, Joun Wiley & Sons.
  3. Liao, S. Beyond perturbation introduction to homotopy analysis method, 2003, CRC Press.
  4. Main, I.G., Vibrations and waves in physics 1993, Cambridge University Press
  5. Olabode, D.L., Lamboni, B., Orou, J.B.C., Active Control of Chaotic Oscillations in Nonlinear Chemical Dynamics, JAMP. 2019, Vol. 7, No. 3, doi: 10.4236/jamp.2019.73040
  6. Oliveira, A.R.E., History of Krylov-Bogoliubov-Mitropolsky Methods of Nonlinear Oscillations, doi: 10.4236/ahs.2017.61003
  7. Sarafian, H., Dynamic Dipole-Dipole Magnetic Interaction and Damped Nonlinear Oscillations, JEMAA, 2009, Vol. 1, No. 4, doi: 10.4236/jemaa.2009.14030
  8. Sarafian, H., Nonlinear Oscillations of a Magneto Static Spring-Mass, JEMAA, 2011, Vol. 3, No. 5, doi: 10.4236/jemaa.2011.35022


Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions