# The Bridge to Nonlinearity
## Introduction: When Linear Theory Breaks Down
In Module 1, we studied the harmonic oscillator where the restoring force is proportional to displacement: $F = -ky$. This linearity gave us exact sinusoidal solutions, constant frequency independent of amplitude, and the principle of superposition. But real physical systems rarely maintain perfect linearity, especially for larger displacements.
In this module, we explore what happens when we include the next term in the Taylor series expansion of the restoring force—leading us from the familiar pendulum to the celebrated Duffing oscillator.
## From Pendulum to Duffing: A Natural Progression
### The Exact Pendulum Equation
Consider a mass $m$ attached to a rigid, weightless rod of length $L$, free to rotate in a vertical plane about a frictionless pivot. The equation of motion is:
$mL\ddot{\theta} + \gamma L\dot{\theta} + mg\sin(\theta) = f(t)$
where:
- $\theta$ = angular displacement from vertical (positive counterclockwise)
- $\gamma$ = Stokes' drag coefficient
- $g$ = gravitational acceleration
- $f(t)$ = external forcing
![[Teaching/MATH310/Lecture Notes/Nonlinear Oscillators/Oscillators - Bridge to Nonlinearity/pendulum_diagram.png]]
*Figure: Forces acting on a pendulum showing the gravitational components and angular displacement*
The nonlinearity enters through $\sin(\theta)$—a transcendental function that makes exact analytical solutions impossible for large amplitudes.
### The Hierarchy of Approximations
To understand this nonlinearity, we expand $\sin(\theta)$ in a [Taylor series](https://en.wikipedia.org/wiki/Taylor_series):
$\sin(\theta) = \theta - \frac{\theta^3}{3!} + \frac{\theta^5}{5!} - \cdots$
This leads to a hierarchy of approximations:
#### 1. Linear (Small Angle) Approximation
For $|\theta| \ll 1$, keeping only the first term:
$\ddot{\theta} + \omega^2\theta = 0, \quad \text{where } \omega^2 = \frac{g}{L}$
This is our familiar harmonic oscillator with frequency $\omega$ independent of amplitude.
#### 2. Cubic (Duffing) Approximation
Keeping the first two terms:
$\ddot{\theta} + \omega^2\left(\theta - \frac{\theta^3}{6}\right) = 0$
Rescaling by defining $y = \theta$ and rearranging:
$m\ddot{y} + k_1 y + k_3 y^3 = 0$
where $k_1 = mg/L$ and $k_3 = -mg/(6L)$.
This is the **Duffing equation**—the simplest nonlinear oscillator that captures essential features missing from linear theory.
## The Duffing Oscillator: A Universal Model
### General Form
The Duffing oscillator represents any system with a cubic nonlinearity:
$m\ddot{y} + \gamma\dot{y} + k_1 y + k_3 y^3 = f(t)$
The force law:
$F_s(y) = -k_1 y - k_3 y^3$
represents the next odd term (preserving symmetry) in the Taylor expansion of a general [nonlinear spring](https://en.wikipedia.org/wiki/Nonlinear_system) force.
### Classification: Soft vs. Hard Springs
The sign of $k_3$ fundamentally changes the system's behavior:
#### Hard Spring ($k_3 > 0$)
- Restoring force **increases** faster than linear for large displacements
- Spring "hardens" or "supports" the linear term
- Frequency **increases** with amplitude
- Physical example: Stretched elastic materials approaching their limit
#### Soft Spring ($k_3 < 0$)
- Restoring force **increases** slower than linear for large displacements
- Cubic term "retards" the linear term
- Frequency **decreases** with amplitude
- Physical example: Pendulum for moderate angles
The pendulum naturally gives a soft spring with $k_3 = -mg/(6L) < 0$.
## Energy and Phase Space Structure
### The Nonlinear Hamiltonian
For the conservative Duffing oscillator ($\gamma = 0$, $f(t) = 0$):
$H(y, \dot{y}) = \frac{m\dot{y}^2}{2} + \frac{k_1 y^2}{2} + \frac{k_3 y^4}{4}$
This energy function consists of:
- Kinetic energy: $T = \frac{m\dot{y}^2}{2}$
- Potential energy: $V(y) = \frac{k_1 y^2}{2} + \frac{k_3 y^4}{4}$
### Numerical Demonstration: Energy Surfaces
The Hamiltonian structure reveals how energy determines the allowed phase space trajectories.
#### Implementation
## Python
Step 1. Define the energy functions
```python
import numpy as np
import matplotlib.pyplot as plt
class DuffingEnergy:
def __init__(self, k1=1.0, k3=0.1, m=1.0):
self.k1 = k1
self.k3 = k3
self.m = m
def potential(self, y):
"""Potential energy V(y)"""
return 0.5*self.k1*y**2 + 0.25*self.k3*y**4
def hamiltonian(self, y, v):
"""Total energy H(y,v) = T + V"""
return 0.5*self.m*v**2 + self.potential(y)
```
Step 2. Visualize phase space contours
```python
# Create mesh
y = np.linspace(-2.5, 2.5, 100)
v = np.linspace(-2, 2, 100)
Y, V = np.meshgrid(y, v)
# Calculate energy surface
duff = DuffingEnergy(k1=-1.0, k3=0.5) # Double well
H = duff.hamiltonian(Y, V)
# Plot contours
plt.contour(Y, V, H, levels=15)
plt.xlabel('Position y')
plt.ylabel('Velocity v')
plt.title('Phase Space Energy Contours')
plt.grid(True)
plt.show()
```
## R
Step 1. Define the energy functions
```r
# Create Duffing energy calculator
create_duffing_energy <- function(k1 = 1.0, k3 = 0.1, m = 1.0) {
list(k1 = k1, k3 = k3, m = m)
}
# Potential energy V(y)
potential <- function(duff, y) {
0.5 * duff$k1 * y^2 + 0.25 * duff$k3 * y^4
}
# Total energy H(y,v) = T + V
hamiltonian <- function(duff, y, v) {
0.5 * duff$m * v^2 + potential(duff, y)
}
```
Step 2. Visualize phase space contours
```r
# Create mesh
y <- seq(-2.5, 2.5, length.out = 100)
v <- seq(-2, 2, length.out = 100)
mesh <- expand.grid(y = y, v = v)
# Calculate energy surface
duff <- create_duffing_energy(k1 = -1.0, k3 = 0.5) # Double well
H <- hamiltonian(duff, mesh$y, mesh$v)
H_matrix <- matrix(H, nrow = length(y), ncol = length(v))
# Plot contours
contour(y, v, H_matrix, nlevels = 15,
xlab = "Position y", ylab = "Velocity v",
main = "Phase Space Energy Contours")
grid()
```
## MATLAB
Step 1. Define the energy functions
```matlab
function duff = create_duffing_energy(k1, k3, m)
if nargin < 1, k1 = 1.0; end
if nargin < 2, k3 = 0.1; end
if nargin < 3, m = 1.0; end
duff.k1 = k1;
duff.k3 = k3;
duff.m = m;
end
function V = potential(duff, y)
V = 0.5 * duff.k1 * y.^2 + 0.25 * duff.k3 * y.^4;
end
function H = hamiltonian(duff, y, v)
H = 0.5 * duff.m * v.^2 + potential(duff, y);
end
```
Step 2. Visualize phase space contours
```matlab
% Create mesh
y = linspace(-2.5, 2.5, 100);
v = linspace(-2, 2, 100);
[Y, V] = meshgrid(y, v);
% Calculate energy surface
duff = create_duffing_energy(-1.0, 0.5, 1.0); % Double well
H = hamiltonian(duff, Y, V);
% Plot contours
contour(Y, V, H, 15);
xlabel('Position y');
ylabel('Velocity v');
title('Phase Space Energy Contours');
grid on;
```
#### Results
![[Teaching/MATH310/Lecture Notes/Nonlinear Oscillators/Oscillators - Bridge to Nonlinearity/phase_space_verification.png]]
*Figure: Phase space structures for different spring types, showing how nonlinearity affects the shape of trajectories.*
The energy surfaces reveal distinct behaviors:
- **Single well (k₁ > 0)**: Closed orbits at all energies
- **Double well (k₁ < 0, k₃ > 0)**: Separatrix divides phase space
- **Energy conservation**: Trajectories follow constant energy contours
**Complete code**: `Module2_Bridge_to_Nonlinearity/codes/duffing_energy_surfaces_FAST.py`
#### Exploration Exercises
1. **3D visualization**: Plot the Hamiltonian as a 3D surface.
2. **Poincaré sections**: Sample the phase space at regular intervals.
3. **Action-angle variables**: Transform to canonical coordinates.
### Potential Energy Landscapes
The potential $V(y)$ determines the system's qualitative behavior:
#### Case 1: $k_1 > 0$ (Single Well)
- Potential has a single minimum at $y = 0$
- All trajectories are bounded oscillations
- Phase portrait shows nested closed curves (centers)
#### Case 2: $k_1 < 0, k_3 > 0$ (Double Well)
- Potential has two minima at $y = \pm\sqrt{-k_1/k_3}$
- Local maximum at $y = 0$ (unstable equilibrium)
- Phase portrait shows:
- Small oscillations in each well
- Large oscillations over both wells
- Separatrix dividing these regions
### Phase Space Trajectories
From energy conservation, trajectories satisfy:
$H(y, \dot{y}) = E = \text{constant}$
Solving for velocity:
$\dot{y} = \pm\sqrt{\frac{2(E - V(y))}{m}}$
This gives the phase space trajectory implicitly. The ± sign indicates clockwise motion in phase space.
### Critical Points and Stability
Equilibrium points occur where $\dot{y} = 0$ and $\ddot{y} = 0$:
$k_1 y + k_3 y^3 = 0$
Solutions:
- $y = 0$ (always an equilibrium)
- $y = \pm\sqrt{-k_1/k_3}$ (exist only if $k_1$ and $k_3$ have opposite signs)
Stability determined by the linearized system:
- **Center**: Stable equilibrium with purely imaginary eigenvalues (oscillatory)
- **Saddle**: Unstable equilibrium with real eigenvalues of opposite sign
- **Separatrix**: Special trajectory connecting saddle points, dividing phase space into regions with qualitatively different motion
![[Teaching/MATH310/Lecture Notes/Nonlinear Oscillators/Oscillators - Bridge to Nonlinearity/phase_portraits_comparison.png]]
*Figure: Phase portrait showing centers, saddle points, and the separatrix in a nonlinear system*
## Amplitude-Dependent Phenomena
### Nonlinear Frequency Shift
Unlike the harmonic oscillator, the Duffing oscillator's frequency depends on amplitude. We can derive this through perturbation analysis.
#### Perturbation Expansion Method
Starting with the normalized Duffing equation where the potential is:
$V(q) = \frac{1}{2}m\omega_0^2 q^2 + \frac{1}{4}\varepsilon m q^4$
We expand the solution in powers of the perturbation parameter $\varepsilon$:
$q(t) = q_0(t) + \varepsilon q_1(t) + \varepsilon^2 q_2(t) + \cdots$
This leads to a hierarchy of linear equations:
- Order $\varepsilon^0$: $\ddot{q}_0 + \omega_0^2 q_0 = 0$
- Order $\varepsilon^1$: $\ddot{q}_1 + \omega_0^2 q_1 + q_0^3 = 0$
#### The Resonance Problem
For initial conditions $q_0(0) = a$ and $\dot{q}_0(0) = 0$, we have $q_0(t) = a\cos(\omega_0 t)$.
Substituting into the $\varepsilon^1$ equation:
$\ddot{q}_1 + \omega_0^2 q_1 = -a^3\cos^3(\omega_0 t) = -\frac{a^3}{4}[3\cos(\omega_0 t) + \cos(3\omega_0 t)]$
The dangerous term $3\cos(\omega_0 t)$ resonates with the natural frequency, leading to a solution containing:
$q_1(t) \sim -\frac{3a^3}{8\omega_0} t \sin(\omega_0 t)$
This **secular term** grows linearly with time, violating energy conservation!
#### Frequency Renormalization
To resolve this, we recognize that the frequency itself must be modified by the nonlinearity:
$\omega(\varepsilon) = \omega_0 + \varepsilon\omega_1 + \varepsilon^2\omega_2 + \cdots$
By choosing:
$\omega_1 = \frac{3a^2}{8\omega_0}$
we suppress the resonance, yielding the **amplitude-dependent frequency**:
$\boxed{\omega(a) = \omega_0 + \varepsilon\frac{3a^2}{8\omega_0} + O(\varepsilon^2)}$
### Numerical Demonstration: Amplitude-Frequency Relationship
The frequency shift can be verified numerically by measuring oscillation periods at different amplitudes.
#### Implementation Strategy
## Python
Step 1. Define the Duffing oscillator class
```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks
class DuffingOscillator:
def __init__(self, k1=1.0, k3=0.1, m=1.0):
self.k1 = k1
self.k3 = k3
self.m = m
self.omega0 = np.sqrt(k1/m) if k1 > 0 else 0
def dynamics(self, t, state):
"""ODE system for Duffing oscillator"""
y, v = state
return [v, -(self.k1*y + self.k3*y**3)/self.m]
def perturbation_frequency(self, amplitude):
"""Theoretical frequency from perturbation theory"""
if self.omega0 == 0:
return 0
epsilon = self.k3/self.k1
return self.omega0 + epsilon * 3*amplitude**2/(8*self.omega0)
```
Step 2. Measure frequency numerically
```python
def measure_frequency(self, amplitude, t_max=50):
"""Measure frequency using peak detection"""
# Solve ODE
sol = solve_ivp(self.dynamics, [0, t_max], [amplitude, 0],
method='RK45', rtol=1e-8)
# Find peaks and compute period
peaks, _ = find_peaks(sol.y[0], height=amplitude*0.8)
if len(peaks) > 1:
period = np.mean(np.diff(sol.t[peaks]))
return 2*np.pi/period
return self.omega0
```
Step 3. Compare hard and soft springs
```python
# Test parameters
oscillators = {
'Linear': DuffingOscillator(k1=1.0, k3=0.0),
'Hard Spring': DuffingOscillator(k1=1.0, k3=0.5),
'Soft Spring': DuffingOscillator(k1=1.0, k3=-0.3)
}
amplitudes = np.linspace(0.1, 2.0, 10)
plt.figure(figsize=(10, 6))
for name, osc in oscillators.items():
frequencies = [osc.measure_frequency(a) for a in amplitudes]
plt.plot(amplitudes, frequencies, 'o-', label=name, linewidth=2)
plt.xlabel('Amplitude')
plt.ylabel('Frequency (rad/s)')
plt.title('Amplitude-Frequency Relationship')
plt.legend()
plt.grid(True)
plt.show()
```
## R
Step 1. Define the Duffing oscillator functions
```r
library(deSolve)
library(pracma)
# Create Duffing oscillator
create_duffing <- function(k1 = 1.0, k3 = 0.1, m = 1.0) {
list(
k1 = k1,
k3 = k3,
m = m,
omega0 = if (k1 > 0) sqrt(k1/m) else 0
)
}
# ODE dynamics
duffing_dynamics <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
dy <- v
dv <- -(k1*y + k3*y^3)/m
list(c(dy, dv))
})
}
# Theoretical frequency from perturbation theory
perturbation_frequency <- function(osc, amplitude) {
if (osc$omega0 == 0) return(0)
epsilon <- osc$k3/osc$k1
osc$omega0 + epsilon * 3*amplitude^2/(8*osc$omega0)
}
```
Step 2. Measure frequency numerically
```r
measure_frequency <- function(osc, amplitude, t_max = 50) {
# Initial conditions
y0 <- c(y = amplitude, v = 0)
times <- seq(0, t_max, length.out = 5000)
# Parameters for ODE
parms <- list(k1 = osc$k1, k3 = osc$k3, m = osc$m)
# Solve ODE
sol <- ode(y = y0, times = times, func = duffing_dynamics, parms = parms,
method = "ode45", rtol = 1e-8)
# Find peaks to determine period
peaks <- findpeaks(sol[, "y"], minpeakheight = amplitude * 0.8)
if (!is.null(peaks) && nrow(peaks) > 1) {
peak_indices <- peaks[, 2]
periods <- diff(times[peak_indices])
avg_period <- mean(periods)
frequency <- 2*pi/avg_period
return(frequency)
}
return(osc$omega0)
}
```
Step 3. Compare hard and soft springs
```r
# Create oscillators
oscillators <- list(
Linear = create_duffing(k1 = 1.0, k3 = 0.0),
HardSpring = create_duffing(k1 = 1.0, k3 = 0.5),
SoftSpring = create_duffing(k1 = 1.0, k3 = -0.3)
)
amplitudes <- seq(0.1, 2.0, length.out = 10)
# Plot setup
plot(NULL, xlim = range(amplitudes), ylim = c(0.5, 1.5),
xlab = "Amplitude", ylab = "Frequency (rad/s)",
main = "Amplitude-Frequency Relationship")
colors <- c("black", "red", "blue")
for (i in seq_along(oscillators)) {
osc <- oscillators[[i]]
frequencies <- sapply(amplitudes, function(a) measure_frequency(osc, a))
lines(amplitudes, frequencies, col = colors[i], lwd = 2)
points(amplitudes, frequencies, col = colors[i], pch = 19)
}
legend("topright", names(oscillators), col = colors, lwd = 2, pch = 19)
grid()
```
## MATLAB
Step 1. Define the Duffing oscillator functions
```matlab
function osc = create_duffing(k1, k3, m)
if nargin < 1, k1 = 1.0; end
if nargin < 2, k3 = 0.1; end
if nargin < 3, m = 1.0; end
osc.k1 = k1;
osc.k3 = k3;
osc.m = m;
if k1 > 0
osc.omega0 = sqrt(k1/m);
else
osc.omega0 = 0;
end
end
function dydt = duffing_dynamics(t, state, osc)
y = state(1);
v = state(2);
dydt = zeros(2,1);
dydt(1) = v;
dydt(2) = -(osc.k1*y + osc.k3*y^3)/osc.m;
end
function freq = perturbation_frequency(osc, amplitude)
if osc.omega0 == 0
freq = 0;
return;
end
epsilon = osc.k3/osc.k1;
freq = osc.omega0 + epsilon * 3*amplitude^2/(8*osc.omega0);
end
```
Step 2. Measure frequency numerically
```matlab
function frequency = measure_frequency(osc, amplitude, t_max)
if nargin < 3, t_max = 50; end
% Initial conditions
y0 = [amplitude; 0];
% Solve ODE
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
[t, y] = ode45(@(t, state) duffing_dynamics(t, state, osc), ...
[0 t_max], y0, options);
% Find peaks to determine period
[pks, locs] = findpeaks(y(:,1), 'MinPeakHeight', amplitude*0.8);
if length(locs) > 1
periods = diff(t(locs));
avg_period = mean(periods);
frequency = 2*pi/avg_period;
else
frequency = osc.omega0;
end
end
```
Step 3. Compare hard and soft springs
```matlab
% Create oscillators
oscillators = struct();
oscillators.Linear = create_duffing(1.0, 0.0, 1.0);
oscillators.HardSpring = create_duffing(1.0, 0.5, 1.0);
oscillators.SoftSpring = create_duffing(1.0, -0.3, 1.0);
amplitudes = linspace(0.1, 2.0, 10);
figure;
hold on;
names = fieldnames(oscillators);
colors = {'k', 'r', 'b'};
for i = 1:length(names)
osc = oscillators.(names{i});
frequencies = zeros(size(amplitudes));
for j = 1:length(amplitudes)
frequencies(j) = measure_frequency(osc, amplitudes(j));
end
plot(amplitudes, frequencies, 'o-', 'Color', colors{i}, ...
'LineWidth', 2, 'DisplayName', names{i});
end
xlabel('Amplitude');
ylabel('Frequency (rad/s)');
title('Amplitude-Frequency Relationship');
legend('Location', 'best');
grid on;
hold off;
```
#### Results
![[Teaching/MATH310/Lecture Notes/Nonlinear Oscillators/Oscillators - Bridge to Nonlinearity/frequency_amplitude_FAST.png]]
*Figure: Amplitude-dependent frequency for different spring types. Hard springs (k₃>0) show increasing frequency with amplitude, while soft springs (k₃<0) show decreasing frequency.*
The numerical results confirm the theoretical prediction:
- **Hard spring (k₃ > 0)**: Frequency increases up to 35% at large amplitudes
- **Soft spring (k₃ < 0)**: Frequency decreases up to 32% at large amplitudes
- **Linear spring (k₃ = 0)**: Frequency remains constant
**Complete code**: `Module2_Bridge_to_Nonlinearity/codes/duffing_frequency_amplitude_FAST.py`
#### Exploration Exercises
1. **Parameter study**: Modify k₃ in the code to see how the nonlinearity strength affects the frequency shift.
2. **Phase portraits**: Add phase space plots to visualize how orbits change shape with amplitude.
3. **Energy check**: Verify that energy remains constant during oscillations.
#### Physical Consequences
This result has profound implications:
1. **Frequency increases with amplitude** for hard springs ($\varepsilon > 0$)
2. **Frequency decreases with amplitude** for soft springs ($\varepsilon < 0$)
3. **Sound analogy**: If this oscillator produced sound, its pitch would depend on volume—louder means higher (or lower) pitch!
4. **Harmonic generation**: The solution contains a $3\omega$ component:
$q_1(t) = -\frac{a^3}{32\omega_0^2}[\cos(\omega t) - \cos(3\omega t)]$
This third harmonic is a signature of cubic nonlinearity.
### Numerical Demonstration: Harmonic Generation
The cubic nonlinearity generates odd harmonics (3ω, 5ω, ...) with amplitudes that scale as powers of the fundamental amplitude.
#### Implementation
Step 1. Compute the frequency spectrum
```python
def fast_spectrum(y, dt):
"""Compute FFT spectrum"""
N = len(y)
window = 0.5 - 0.5 * np.cos(2*np.pi*np.arange(N)/N)
y_windowed = y * window
Y = fft(y_windowed)
freqs = fftfreq(N, dt)
# Positive frequencies only
pos_mask = freqs > 0
return freqs[pos_mask], np.abs(Y[pos_mask])
```
Step 2. Analyze harmonic content vs amplitude
```python
osc = DuffingHarmonics(k1=1.0, k3=0.5)
amplitudes = [0.5, 1.0, 1.5]
for amp in amplitudes:
# Solve ODE
sol = solve_ivp(osc.dynamics, [0, 20], [amp, 0],
method='DOP853', rtol=1e-6)
# Compute spectrum
freqs, spectrum = osc.fast_spectrum(sol.y[0], dt)
# Plot normalized to fundamental
plt.semilogy(freqs/osc.omega0, spectrum/max(spectrum), label=f'a={amp}')
```
#### Results
![[Teaching/MATH310/Lecture Notes/Nonlinear Oscillators/Oscillators - Bridge to Nonlinearity/harmonics_FAST.png]]
*Figure: Harmonic spectrum for different amplitudes. Cubic nonlinearity generates odd harmonics with amplitudes that scale as a³, a⁵, etc.*
Key observations:
- **Odd harmonics only**: Cubic nonlinearity preserves symmetry
- **Scaling laws**: 3rd harmonic ∝ a³, 5th harmonic ∝ a⁵
- **Waveform distortion**: Larger amplitudes show more harmonics
**Complete code**: `Module2_Bridge_to_Nonlinearity/codes/duffing_harmonics_FAST.py`
#### Exploration Exercises
1. **Even nonlinearity**: Add a quadratic term (k₂y²) and observe even harmonics.
2. **Forced response**: Drive the oscillator and see harmonic mixing.
3. **Spectral cascade**: Track how energy flows to higher frequencies.
### Exact Solutions via Elliptic Functions
While perturbation theory works for weak nonlinearity, the exact solution of the Duffing equation reveals a profound connection to one of mathematics' most beautiful subjects: elliptic functions.
#### From Duffing to Elliptic Integrals
Starting from energy conservation:
$H = \frac{m\dot{y}^2}{2} + \frac{k_1 y^2}{2} + \frac{k_3 y^4}{4} = E$
Solving for velocity:
$\dot{y} = \pm\sqrt{\frac{2(E - V(y))}{m}} = \pm\sqrt{\frac{2E - k_1 y^2 - \frac{k_3 y^4}{2}}{m}}$
Separating variables and integrating:
$\int \frac{dy}{\sqrt{2E - k_1 y^2 - \frac{k_3 y^4}{2}}} = \pm\sqrt{\frac{1}{m}}\int dt$
With appropriate scaling (choosing $m = 1$, $E = 1/2$, $k_3 = -\kappa^2$, and $k_1 = 1 + \kappa^2$), this becomes:
$\boxed{\int \frac{dy}{\sqrt{(1-y^2)(1-\kappa^2 y^2)}} = \pm \int dt}$
This is the **standard form of an elliptic integral of the first kind**.
#### The Jacobi Elliptic Functions
The solution to this integral involves the Jacobi elliptic functions. For our Duffing oscillator:
$y(t) = \text{cn}(u, \kappa)$
where $u = \omega t$ and $\kappa$ is the **elliptic modulus** (determined by initial conditions).
The Jacobi elliptic functions—$\text{sn}(u, \kappa)$, $\text{cn}(u, \kappa)$, and $\text{dn}(u, \kappa)$—are generalizations of trigonometric functions:
- **When $\kappa = 0$** (linear limit):
- $\text{sn}(u, 0) = \sin(u)$
- $\text{cn}(u, 0) = \cos(u)$
- $\text{dn}(u, 0) = 1$
- **When $\kappa = 1$** (separatrix):
- $\text{sn}(u, 1) = \tanh(u)$
- $\text{cn}(u, 1) = \text{sech}(u)$
- $\text{dn}(u, 1) = \text{sech}(u)$
- **When $0 < \kappa < 1$**: True elliptic behavior—periodic but not sinusoidal
#### The Period Formula
The period of oscillation is given by:
$T = 4K(\kappa)\sqrt{\frac{m}{|k_1 + k_3 a^2|}}$
where $K(\kappa)$ is the **complete elliptic integral of the first kind**:
$K(\kappa) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1 - \kappa^2\sin^2\theta}}$
This shows explicitly how the period depends on amplitude through the modulus $\kappa$.
#### Historical Context and Mathematical Significance
##### Origins in Geometry
Elliptic integrals were first studied in the mid-18th century by mathematicians trying to find the arc length of an ellipse—hence the name "elliptic." Notable contributions came from:
- **Fagnano** (1718): Doubling and bisection of elliptic arcs
- **Euler** (1760s): Systematic study of elliptic integrals
- **Legendre** (1811-1817): Comprehensive treatise and standard forms
##### The Revolution of Inversion
The breakthrough came when **Abel** (1827) and **Jacobi** (1829) independently realized that instead of trying to evaluate the integrals, they should study their inverse functions—the elliptic functions. This was analogous to defining $\sin(x)$ as the inverse of $\int \frac{dy}{\sqrt{1-y^2}}$.
##### Profound Connections
**Jacobi's Remarkable Observation** (noted by Arnold in 1997):
> "Jacobi noted, as mathematics' most fascinating property, that in it one and the same function controls both the representations of a whole number as a sum of four squares and the real movement of a pendulum."
This refers to the appearance of theta functions (related to elliptic functions) in:
1. **Number theory**: Lagrange's four-square theorem
2. **Physics**: Exact pendulum solutions
##### Modern Applications
Elliptic functions and their associated curves have found applications far beyond oscillators:
1. **Cryptography**: Elliptic curve cryptography (ECC) provides security for internet communications
2. **Number Theory**:
- Proof of Fermat's Last Theorem (Wiles, 1995)
- Birch and Swinnerton-Dyer conjecture (Millennium Prize Problem)
3. **String Theory**: Modular forms and elliptic functions appear in partition functions
4. **Engineering**: Optimal filter design, signal processing
5. **General Relativity**: Geodesics in Schwarzschild and Kerr spacetimes
#### The Elliptic Functions as a Bridge
The appearance of elliptic functions in the Duffing oscillator reveals a fundamental truth: **nonlinearity connects elementary physics to advanced mathematics**.
While linear oscillators require only high school trigonometry, the cubic term in the Duffing equation opens a door to:
- 19th-century special function theory
- Modern algebraic geometry
- Contemporary cryptography
- Unsolved problems in mathematics
This exemplifies how even simple nonlinearities lead to mathematical structures of extraordinary richness and continuing relevance.
## Connecting to Physical Systems
### The Pendulum Revisited
For the pendulum with moderate amplitudes, the Duffing approximation:
$\ddot{\theta} + \omega^2\theta - \frac{\omega^2\theta^3}{6} = 0$
captures the essential nonlinearity while remaining analytically tractable. The soft spring character ($k_3 < 0$) correctly predicts that larger swings have longer periods.
### Numerical Demonstration: Pendulum Approximations
We can compare the exact pendulum equation with its linear and Duffing approximations to see the range of validity.
#### Implementation
Step 1. Define the three pendulum models
```python
class ExactPendulum(PendulumModel):
def dynamics(self, t, state):
theta, theta_dot = state
theta_ddot = -(self.g/self.L * np.sin(theta))
return [theta_dot, theta_ddot]
class LinearPendulum(PendulumModel):
def dynamics(self, t, state):
theta, theta_dot = state
theta_ddot = -(self.omega0**2 * theta)
return [theta_dot, theta_ddot]
class DuffingPendulum(PendulumModel):
def dynamics(self, t, state):
theta, theta_dot = state
theta_ddot = -(self.omega0**2 * (theta - theta**3/6))
return [theta_dot, theta_ddot]
```
Step 2. Compare period vs amplitude
```python
theta0_range = np.array([0.1, 0.3, 0.5, 0.8, 1.0, 1.5, 2.0, 2.5])
models = {'Exact': ExactPendulum(),
'Linear': LinearPendulum(),
'Duffing': DuffingPendulum()}
for name, model in models.items():
periods = [model.fast_measure_period(theta0) for theta0 in theta0_range]
plt.plot(np.degrees(theta0_range), periods/T0, 'o-', label=name)
```
#### Results
![[Teaching/MATH310/Lecture Notes/Nonlinear Oscillators/Oscillators - Bridge to Nonlinearity/pendulum_approximations_FAST.png]]
*Figure: Comparison of pendulum approximations. The Duffing (cubic) approximation extends the valid range significantly beyond the linear approximation.*
![[Teaching/MATH310/Lecture Notes/Nonlinear Oscillators/Oscillators - Bridge to Nonlinearity/taylor_series_FAST.png]]
*Figure: Taylor series approximations of sin(θ). The cubic term in the Duffing approximation captures the nonlinearity up to about 60°.*
Key findings:
- **Linear approximation**: Valid for θ < 15° (error < 1%)
- **Duffing approximation**: Valid for θ < 60° (error < 1%)
- **Period increase**: At 90°, period is 18% longer than small-angle prediction
**Complete code**: `Module2_Bridge_to_Nonlinearity/codes/pendulum_approximations_FAST.py`
#### Exploration Exercises
1. **Higher-order terms**: Add the fifth-order term (θ⁵/120) and see the improvement.
2. **Phase portraits**: Compare the phase space trajectories for the three models.
3. **Energy surfaces**: Plot the potential energy for each approximation.
### Beyond Pendulums
The Duffing equation models numerous physical systems:
- **Mechanical**: Beams and plates under large deflections
- **Electrical**: Nonlinear inductors in circuits
- **Optical**: Nonlinear refractive index in intense laser fields
- **Biological**: Nonlinear elasticity in tissues
- **Architectural**: Building sway under wind loads
## Mathematical Structure and Symmetry
### Symmetry Properties
The cubic nonlinearity preserves the system's symmetry:
- If $y(t)$ is a solution, so is $-y(t)$
- The phase portrait is symmetric about both axes
- Only odd powers of $y$ appear in the force law
### Loss of Superposition
A fundamental consequence of nonlinearity: the principle of superposition fails.
- If $y_1(t)$ and $y_2(t)$ are solutions, $y_1(t) + y_2(t)$ is generally **not** a solution
- Linear combinations of normal modes do not yield new normal modes
- Energy does not partition cleanly between modes
### Qualitative Changes with Parameters
As parameters vary, the system undergoes qualitative changes:
- **Bifurcations**: Critical parameter values where the number or stability of equilibria changes
- **Hysteresis**: System response depends on history (increasing vs. decreasing parameter)
- **Jump phenomena**: Sudden transitions between oscillation states
## Summary: The Nonlinear Difference
The transition from linear to [nonlinear oscillators](https://en.wikipedia.org/wiki/Nonlinear_system) reveals fundamentally new phenomena:
1. **Amplitude-dependent frequency**: Period depends on energy
2. **Multiple [equilibria](https://en.wikipedia.org/wiki/Equilibrium_point)**: Possibility of coexisting stable states
3. **[Separatrices](https://en.wikipedia.org/wiki/Separatrix_(dynamical_systems))**: Boundaries between qualitatively different motions
4. **Special functions**: [Elliptic functions](https://en.wikipedia.org/wiki/Elliptic_function) replace trigonometric ones
5. **Rich [phase space](https://en.wikipedia.org/wiki/Phase_space)**: Centers, [saddles](https://en.wikipedia.org/wiki/Saddle_point), and complex trajectory structures
These features, absent in linear systems, make nonlinear oscillators both challenging and fascinating. The Duffing oscillator, as the simplest model exhibiting these behaviors, serves as our gateway to understanding more complex nonlinear dynamics, including chaos—topics we'll explore in subsequent modules.
## Numerical Codes Summary
All demonstration codes are optimized for educational use, running in seconds rather than minutes while maintaining physical accuracy.
### Available Python Codes
1. **`duffing_frequency_amplitude_FAST.py`** - Amplitude-frequency relationships
- Demonstrates frequency shifts for hard/soft springs
- Compares numerical results with perturbation theory
- Runtime: ~5 seconds
2. **`duffing_phase_portraits.py`** - Phase space structures
- Shows single well, double well, and separatrices
- Visualizes different parameter regimes
- Runtime: ~8 seconds
3. **`pendulum_approximations_FAST.py`** - Pendulum approximation comparison
- Exact vs linear vs Duffing approximations
- Period-amplitude relationships
- Runtime: ~6 seconds
4. **`duffing_harmonics_FAST.py`** - Harmonic generation analysis
- FFT spectrum analysis
- Demonstrates odd harmonic generation
- Runtime: ~7 seconds
5. **`duffing_energy_surfaces_FAST.py`** - Energy and Hamiltonian visualization
- 3D energy surfaces
- Phase space contours
- Runtime: ~5 seconds
6. **`elliptic_functions_demo_FAST.py`** - Connection to elliptic functions
- Jacobi elliptic functions
- Transition from periodic to soliton
- Runtime: ~4 seconds
### Running the Codes
```bash
cd Module2_Bridge_to_Nonlinearity/codes
python3 duffing_frequency_amplitude_FAST.py
```
### Suggested Explorations
1. **Parameter Study**: Vary k₃ from -1 to +1 to see the transition from soft to hard spring
2. **Initial Conditions**: Try different initial velocities to explore different energy levels
3. **Damping Effects**: Add small damping (γ = 0.1) to see energy decay
4. **Forcing**: Add periodic forcing to observe resonance and chaos (preview of Module 3)
## Key Takeaways
- The Duffing oscillator naturally emerges from the pendulum equation through systematic approximation
- The sign of the cubic coefficient ($k_3$) determines whether the spring is "soft" or "hard"
- Nonlinearity introduces amplitude-dependent frequencies and complex phase space structures
- Exact solutions require elliptic functions, reflecting the fundamental difference from linear systems
- The loss of superposition makes nonlinear systems qualitatively different, not just quantitatively more complex
### Verified Relationships
Our numerical demonstrations confirm:
- **Hard spring (k₃ > 0)**: Frequency increases with amplitude (up to 35% at a = 1.5)
- **Soft spring (k₃ < 0)**: Frequency decreases with amplitude (down to 32% at a = 1.5)
- **Harmonic generation**: 3rd harmonic amplitude ∝ a³, 5th harmonic ∝ a⁵
- **Theory agreement**: Perturbation theory is accurate to within 2% for moderate amplitudes
The mathematical tools developed here—phase portraits, potential energy landscapes, and perturbation methods—will be essential as we explore driven systems and the route to chaos in Module 3.