Imports
/- Copyright (c) 2026 Nicola Bernini. All rights reserved. Released under Apache 2.0 license as described in the file LICENSE. Authors: Nicola Bernini, Florian Wiesner -/ module public import Physlib.ClassicalMechanics.HarmonicOscillator.Basic

The damped harmonic oscillator

i. Overview

The damped harmonic oscillator is a classical mechanical system consisting of a mass m under a restoring force - k x and a damping force - γ ẋ, where k is the spring constant, γ is the damping coefficient, x is the position, and is the velocity.

The equation of motion for the damped harmonic oscillator is:

m ẍ + γ ẋ + k x = 0

Depending on the relationship between the damping coefficient and the natural frequency, the system exhibits three different behaviors:

    Underdamped (γ^2 < 4 * m * k) : oscillatory motion with exponentially decaying amplitude.

    Critically damped (γ^2 = 4 * m * k) : fastest return to equilibrium without oscillation.

    Overdamped (4 * m * k < γ^2) : slow return to equilibrium without oscillation.

In this file, the position and velocity both have type EuclideanSpace ℝ (Fin 1). This coordinate model is useful for a first formalization, but it works only because the one-dimensional configuration space and its tangent space are both isomorphic to one-dimensional Euclidean space. A more geometric formalization should represent the configuration space and its tangent bundle directly.

ii. Key results

The key results in the study of the classical damped harmonic oscillator are the following:

In the Basic module:

    DampedHarmonicOscillator contains the input data to the problem.

    EquationOfMotion defines the damped oscillator equation m ẍ + γ ẋ + k x = 0.

    energy_dissipation_rate computes the rate at which damping removes mechanical energy.

    IsUnderdamped, IsCriticallyDamped, and IsOverdamped define the three damping regimes from the discriminant γ^2 - 4 * m * k.

    angularFrequency selects the real frequency parameter from the damping regime.

    toUndamped_equationOfMotion relates the damped and undamped equations of motion when the damping coefficient is zero.

    lagrangian defines the Caldirola–Kanai lagrangian exp (γ/m * t) * (T - V).

    gradLagrangian defines the variational gradient of the corresponding action, and equationOfMotion_iff_gradLagrangian_zero shows that its vanishing is equivalent to the equation of motion.

In the Solution module:

    InitialConditions contains the initial position and velocity.

    trajectory gives the explicit solution selected from the damping regime.

iii. Table of contents

    A. The input data

    B. The equation of motion and energy dissipation

      B.1. The equation of motion

      B.2. Energy dissipation

    C. Newton's second law

      C.1. The force

      C.2. Equation of motion if and only if Newton's second law

    D. Damping regimes

    E. To undamped oscillator

    F. The Caldirola–Kanai lagrangian and the equation of motion

      F.1. The lagrangian

        F.1.1. Equalities for the lagrangian

        F.1.2. Smoothness of the lagrangian

        F.1.3. Gradients of the lagrangian

      F.2. The variational gradient of the action

        F.2.1. The variational gradient and Newton's second law

      F.3. Equation of motion iff the variational gradient vanishes

iv. References

References for the damped harmonic oscillator include:

    Landau & Lifshitz, Mechanics, page 76, section 25.

    Goldstein, Classical Mechanics, Chapter 2.

References for the Caldirola–Kanai lagrangian include:

    Caldirola, Nuovo Cimento 18 (1941) 393.

    Kanai, Progress of Theoretical Physics 3 (1948) 440.

@[expose] public sectionTODO "Create a new file for the geometric model which properly models the position as a configuration space and velocity as its tangent space, see the HarmonicOscillator file."TODO "Define and prove properties of the quality factor Q."TODO "Define and prove properties of the relaxation time τ."

A. The input data

We start by defining a structure containing the input data of the damped harmonic oscillator. The mass m and spring constant k are inherited from HarmonicOscillator; this file adds the damping coefficient γ.

The classical damped harmonic oscillator is specified by a mass m, a spring constant k, and a damping coefficient γ.

The mass and spring constant are inherited from HarmonicOscillator and are positive. The damping coefficient is assumed to be nonnegative.

The damping coefficient of the oscillator.

The damping coefficient is nonnegative.

@[ext] structure DampedHarmonicOscillator extends HarmonicOscillator where γ : γ_nonneg : 0 γ

The mass/spring nonzero lemmas, the natural angular frequency, and the undamped energy API are inherited from HarmonicOscillator.

B. The equation of motion and energy dissipation

B.1. The equation of motion

Solving the equation of motion for the acceleration: along a solution the second derivative is -(k/m) x - (γ/m) ẋ.

S:DampedHarmonicOscillatorz:Time EuclideanSpace (Fin 1)hEOM:S.EquationOfMotion zt:Timehm:S.m 0hsum:S.m ∂ₜ (∂ₜ z) t + (S.γ ∂ₜ z t + S.k z t) = 0hma:S.m ∂ₜ (∂ₜ z) t = -(S.k z t) - S.γ ∂ₜ z thkey:∂ₜ (∂ₜ z) t = S.m⁻¹ S.m ∂ₜ (∂ₜ z) tS.m⁻¹ (-(S.k z t) - S.γ ∂ₜ z t) = -(S.m⁻¹ * S.k) z t + -(S.m⁻¹ * S.γ) ∂ₜ z t All goals completed! 🐙

B.2. Energy dissipation

The damped oscillator inherits the mechanical energy from the undamped harmonic oscillator. Along a solution of the damped equation of motion, that energy decreases at a rate proportional to -γ ‖ẋ‖^2.

Along a smooth solution of the damped equation of motion, the derivative of the mechanical energy is -γ ‖ẋ‖^2.

S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)t:Timeh1:S.EquationOfMotion xₜhx:ContDiff xₜhforce:S.m ∂ₜ (∂ₜ xₜ) t + S.k xₜ t = -S.γ ∂ₜ xₜ t(fun t => ∂ₜ xₜ t, S.m ∂ₜ (∂ₜ xₜ) t + S.k xₜ t⟫_) t = -S.γ * ∂ₜ xₜ t, ∂ₜ xₜ t⟫_ All goals completed! 🐙

If 0 < γ and the velocity is nonzero at a time, the mechanical energy is strictly decreasing at that time.

S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)t:Timeh1:S.EquationOfMotion xₜhx:ContDiff xₜhdx:∂ₜ xₜ t 0:0 < S.γ-(S.γ * ∂ₜ xₜ t, ∂ₜ xₜ t⟫_) < 0 All goals completed! 🐙

C. Newton's second law

We define the force of the damped oscillator, and show that the equation of motion is equivalent to Newton's second law.

C.1. The force

We define the force of the damped oscillator as - k x - γ v.

C.2. Equation of motion if and only if Newton's second law

We show that the equation of motion is equivalent to Newton's second law.

lemma equationOfMotion_iff_newtons_2nd_law (xₜ : Time EuclideanSpace (Fin 1)) : S.EquationOfMotion xₜ ( t : Time, S.m ∂ₜ (∂ₜ xₜ) t = force S xₜ t) := S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)S.EquationOfMotion xₜ (t : Time), S.m ∂ₜ (∂ₜ xₜ) t = S.force xₜ t S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)(∀ (t : Time), S.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t + S.k xₜ t = 0) (t : Time), S.m ∂ₜ (∂ₜ xₜ) t = -S.k xₜ t - S.γ ∂ₜ xₜ t S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)(∀ (t : Time), S.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t + S.k xₜ t = 0) (t : Time), S.m ∂ₜ (∂ₜ xₜ) t = -S.k xₜ t - S.γ ∂ₜ xₜ tS:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)(∀ (t : Time), S.m ∂ₜ (∂ₜ xₜ) t = -S.k xₜ t - S.γ ∂ₜ xₜ t) (t : Time), S.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t + S.k xₜ t = 0 S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)(∀ (t : Time), S.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t + S.k xₜ t = 0) (t : Time), S.m ∂ₜ (∂ₜ xₜ) t = -S.k xₜ t - S.γ ∂ₜ xₜ t S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)h: (t : Time), S.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t + S.k xₜ t = 0t:TimeS.m ∂ₜ (∂ₜ xₜ) t = -S.k xₜ t - S.γ ∂ₜ xₜ t linear_combination (norm := All goals completed! 🐙) h t S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)(∀ (t : Time), S.m ∂ₜ (∂ₜ xₜ) t = -S.k xₜ t - S.γ ∂ₜ xₜ t) (t : Time), S.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t + S.k xₜ t = 0 S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)h: (t : Time), S.m ∂ₜ (∂ₜ xₜ) t = -S.k xₜ t - S.γ ∂ₜ xₜ tt:TimeS.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t + S.k xₜ t = 0 linear_combination (norm := All goals completed! 🐙) h t

D. Damping regimes

The sign of the discriminant γ^2 - 4 * m * k separates the underdamped, critically damped, and overdamped regimes. We also define the decay rate and the regime-selected real frequency that appears in the explicit solution formulas.

The system is underdamped when γ² < 4mk.

def IsUnderdamped : Prop := S.discriminant < 0

The system is critically damped when γ² = 4mk.

def IsCriticallyDamped : Prop := S.discriminant = 0

The system is overdamped when 4mk < γ².

def IsOverdamped : Prop := 0 < S.discriminant

Every damped oscillator is underdamped, critically damped, or overdamped.

lemma isUnderdamped_or_isCriticallyDamped_or_isOverdamped : S.IsUnderdamped S.IsCriticallyDamped S.IsOverdamped := lt_trichotomy S.discriminant 0

The system is undamped when γ = 0.

def IsUndamped : Prop := S.γ = 0

The relationship between the discriminant, decay rate, and natural angular frequency.

S:DampedHarmonicOscillatorS.γ ^ 2 - 4 * S.m * S.k = 4 * S.m ^ 2 * ((S.γ / (2 * S.m)) ^ 2 - S.k / S.m) S:DampedHarmonicOscillator(S.γ ^ 2 - 4 * S.m * S.k) * 2 ^ 2 = 4 * (S.γ ^ 2 - S.m * S.k * 2 ^ 2) All goals completed! 🐙

The decay rate is nonnegative.

lemma decayRate_nonneg : 0 S.decayRate := S:DampedHarmonicOscillator0 S.decayRate exact div_nonneg S.γ_nonneg (S:DampedHarmonicOscillator0 2 * S.m All goals completed! 🐙)

An undamped oscillator lies in the underdamped regime.

S:DampedHarmonicOscillator:S.γ = 00 ^ 2 - 4 * S.m * S.k < 0 All goals completed! 🐙

An underdamped system has decay rate less than the natural frequency.

S:DampedHarmonicOscillatorhS:4 * S.m ^ 2 * (S.decayRate ^ 2 - S.ω ^ 2) < 0hsq:S.decayRate ^ 2 < S.ω ^ 2S.decayRate < S.ω All goals completed! 🐙

A critically damped system has decay rate equal to the natural frequency.

S:DampedHarmonicOscillatorhS:4 * S.m ^ 2 * (S.decayRate ^ 2 - S.ω ^ 2) = 0hsq:S.decayRate ^ 2 = S.ω ^ 2S.ω = S.decayRate All goals completed! 🐙

The damping coefficient is twice mass times the decay rate.

S:DampedHarmonicOscillatorS.γ = 2 * S.m * (S.γ / (2 * S.m)) All goals completed! 🐙

The spring constant is m * ω^2.

S:DampedHarmonicOscillatorS.k = S.m * (S.k / S.m) All goals completed! 🐙

In the critically damped regime, k = m * decayRate^2.

All goals completed! 🐙

An overdamped system has decay rate greater than the natural frequency.

S:DampedHarmonicOscillatorhS:0 < 4 * S.m ^ 2 * (S.decayRate ^ 2 - S.ω ^ 2)hsq:S.ω ^ 2 < S.decayRate ^ 2S.ω < S.decayRate All goals completed! 🐙

In the underdamped regime, the selected frequency uses the oscillation frequency.

lemma angularFrequency_eq_underdamped (hS : S.IsUnderdamped) : S.angularFrequency = sqrt (- S.discriminant) / (2 * S.m) := S:DampedHarmonicOscillatorhS:S.IsUnderdampedS.angularFrequency = (-S.discriminant) / (2 * S.m) All goals completed! 🐙

In the critically damped regime, the selected frequency is zero.

In the overdamped regime, the selected frequency uses the real split rate.

S:DampedHarmonicOscillatorhS:0 < S.discriminantS.angularFrequency = S.discriminant / (2 * S.m) All goals completed! 🐙

In the underdamped regime, the selected angular frequency squares to ω^2 - decayRate^2.

S:DampedHarmonicOscillatorhS:S.IsUnderdamped-(4 * S.m ^ 2 * (S.decayRate ^ 2 - S.ω ^ 2)) / (2 * S.m) ^ 2 = S.ω ^ 2 - S.decayRate ^ 2 S:DampedHarmonicOscillatorhS:S.IsUnderdamped-(4 * (S.decayRate ^ 2 - S.ω ^ 2)) = 2 ^ 2 * (S.ω ^ 2 - S.decayRate ^ 2) All goals completed! 🐙 S:DampedHarmonicOscillatorhS:S.IsUnderdamped0 -S.discriminant All goals completed! 🐙

The selected angular frequency is positive in the underdamped regime.

S:DampedHarmonicOscillatorhS:S.IsUnderdamped0 < (-S.discriminant) / (2 * S.m) exact div_pos (sqrt_pos.mpr (neg_pos.mpr hS)) (S:DampedHarmonicOscillatorhS:S.IsUnderdamped0 < 2 * S.m All goals completed! 🐙)

The selected angular frequency is nonzero in the underdamped regime.

lemma angularFrequency_ne_zero_of_underdamped (hS : S.IsUnderdamped) : S.angularFrequency 0 := (S.angularFrequency_pos_of_underdamped hS).ne'

In the overdamped regime, the selected angular frequency squares to decayRate^2 - ω^2.

S:DampedHarmonicOscillatorhS:S.IsOverdamped4 * S.m ^ 2 * (S.decayRate ^ 2 - S.ω ^ 2) / (2 * S.m) ^ 2 = S.decayRate ^ 2 - S.ω ^ 2 S:DampedHarmonicOscillatorhS:S.IsOverdamped4 * (S.decayRate ^ 2 - S.ω ^ 2) = (S.decayRate ^ 2 - S.ω ^ 2) * 2 ^ 2 All goals completed! 🐙 S:DampedHarmonicOscillatorhS:S.IsOverdamped0 S.discriminant All goals completed! 🐙

The selected angular frequency is positive in the overdamped regime.

S:DampedHarmonicOscillatorhS:S.IsOverdamped0 < S.discriminant / (2 * S.m) exact div_pos (sqrt_pos.mpr hS) (S:DampedHarmonicOscillatorhS:S.IsOverdamped0 < 2 * S.m All goals completed! 🐙)

The selected angular frequency is nonzero in the overdamped regime.

lemma angularFrequency_ne_zero_of_overdamped (hS : S.IsOverdamped) : S.angularFrequency 0 := (S.angularFrequency_pos_of_overdamped hS).ne'

E. To undamped oscillator

We show that the damped harmonic oscillator reduces to the undamped harmonic oscillator when the damping coefficient is zero. The underlying mass and spring data are already inherited from HarmonicOscillator; the proof argument records that this conversion is being used only in the zero-damping case.

We also show that the equations of motion are equivalent in this case.

Convert a damped oscillator to its underlying undamped oscillator when γ = 0.

set_option linter.unusedVariables false in@[nolint unusedArguments] def toUndamped (S : DampedHarmonicOscillator) (_hS : S.IsUndamped) : HarmonicOscillator := S.toHarmonicOscillator

When γ = 0, the damped equation of motion is equivalent to the equation of motion for the corresponding undamped harmonic oscillator.

All goals completed! 🐙

F. The Caldirola–Kanai lagrangian and the equation of motion

The damped harmonic oscillator is not conservative, so the undamped lagrangian T - V does not reproduce the damped equation of motion. Instead we use the Caldirola–Kanai lagrangian, which multiplies the undamped lagrangian by the time-dependent factor exp (γ/m * t):

$$L(t, x, v) = e^{(\gamma/m) t}\left(\frac{1}{2} m |v|^2 - \frac{1}{2} k |x|^2\right)$$

Setting the variational gradient of the corresponding action equal to zero recovers exactly the damped equation of motion m ẍ + γ ẋ + k x = 0.

F.1. The lagrangian

We define the Caldirola–Kanai lagrangian as the lagrangian of the underlying undamped harmonic oscillator multiplied by the exponential factor exp (γ/m * t).

F.1.1. Equalities for the lagrangian

We prove some simple equalities for the lagrangian, in particular that when applied to a trajectory it is the exponential factor times the kinetic energy minus the potential energy, and that it reduces to the undamped lagrangian when the damping coefficient is zero.

All goals completed! 🐙All goals completed! 🐙

When the damping coefficient is zero, the Caldirola–Kanai lagrangian is the lagrangian of the underlying undamped harmonic oscillator.

S:DampedHarmonicOscillatorhS:S.IsUndamped:S.γ = 0S.lagrangian = S.lagrangian S:DampedHarmonicOscillatorhS:S.IsUndamped:S.γ = 0t:Timex:EuclideanSpace (Fin 1)v:EuclideanSpace (Fin 1)S.lagrangian t x v = S.lagrangian t x v All goals completed! 🐙
F.1.2. Smoothness of the lagrangian

The lagrangian is smooth in all its arguments.

@[fun_prop] lemma contDiff_lagrangian (n : WithTop ℕ∞) : ContDiff n S.lagrangian := S:DampedHarmonicOscillatorn:ℕ∞ωContDiff n S.lagrangian S:DampedHarmonicOscillatorn:ℕ∞ωContDiff n fun p => rexp (S.γ / S.m * p.1.val) * S.lagrangian p All goals completed! 🐙
F.1.3. Gradients of the lagrangian

We now show results related to the gradients of the lagrangian with respect to the position and velocity. They follow from the corresponding gradients of the undamped lagrangian, using that the gradient scales with the constant exp (γ/m * t).

private lemma gradient_const_mul {f : EuclideanSpace (Fin 1) } {x : EuclideanSpace (Fin 1)} (c : ) (hf : DifferentiableAt f x) : gradient (fun y => c * f y) x = c gradient f x := f:EuclideanSpace (Fin 1) x:EuclideanSpace (Fin 1)c:hf:DifferentiableAt f xgradient (fun y => c * f y) x = c gradient f x All goals completed! 🐙S:DampedHarmonicOscillatort:Timex:EuclideanSpace (Fin 1)v:EuclideanSpace (Fin 1)hf:DifferentiableAt (fun y => S.lagrangian t y v) xrexp (S.γ / S.m * t.val) -S.k x = -(rexp (S.γ / S.m * t.val) * S.k) x All goals completed! 🐙All goals completed! 🐙

F.2. The variational gradient of the action

We now write down the variational gradient of the action of the damped harmonic oscillator, for a trajectory $x(t)$ this is equal to

$$t\mapsto \left.\frac{\partial L(t, q, \dot x (t))}{\partial q}\right|{q = x(t)} - \frac{d}{dt} \left.\frac{\partial L(t, x(t), v)}{\partial v}\right|{v = \dot x (t)}$$

Setting this equal to zero corresponds to the Euler-Lagrange equations, and thereby the equation of motion.

All goals completed! 🐙
F.2.1. The variational gradient and Newton's second law

We simplify the variational gradient of the action to the exponential factor times the difference between the force and mass times acceleration.

a:y:Time EuclideanSpace (Fin 1)hy:Differentiable yt:Timerexp (a * t.val) ∂ₜ y t + (rexp (a * t.val) * (a + t.val * (fderiv (fun t' => a) t) 1)) y t = rexp (a * t.val) (∂ₜ y t + a y t) All goals completed! 🐙

The variational gradient of the Caldirola–Kanai action is the exponential factor times the difference of the force and mass times acceleration appearing in Newton's second law.

S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)hx:ContDiff xₜhdx:Differentiable (∂ₜ xₜ)t:Timeh2:∂ₜ (fun t' => gradient (fun x => S.lagrangian t' (xₜ t') x) (∂ₜ xₜ t')) t = rexp (S.γ / S.m * t.val) (S.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t)-(rexp (S.γ / S.m * t.val) * S.k) xₜ t - rexp (S.γ / S.m * t.val) (S.m ∂ₜ (∂ₜ xₜ) t + S.γ ∂ₜ xₜ t) = rexp (S.γ / S.m * t.val) (-S.k xₜ t - S.γ ∂ₜ xₜ t - S.m ∂ₜ (∂ₜ xₜ) t) All goals completed! 🐙

F.3. Equation of motion iff the variational gradient vanishes

The equation of motion of the damped harmonic oscillator holds if and only if the variational gradient of the Caldirola–Kanai action vanishes, since the exponential factor is never zero.

S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)hx:ContDiff xₜ(∀ (t : Time), S.m ∂ₜ (∂ₜ xₜ) t = S.force xₜ t) (x : Time), rexp (S.γ / S.m * x.val) (S.force xₜ x - S.m ∂ₜ (∂ₜ xₜ) x) = 0 x S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)hx:ContDiff xₜt:TimeS.m ∂ₜ (∂ₜ xₜ) t = S.force xₜ t rexp (S.γ / S.m * t.val) (S.force xₜ t - S.m ∂ₜ (∂ₜ xₜ) t) = 0 t S:DampedHarmonicOscillatorxₜ:Time EuclideanSpace (Fin 1)hx:ContDiff xₜt:TimeS.m ∂ₜ (∂ₜ xₜ) t = S.force xₜ t S.force xₜ t = S.m ∂ₜ (∂ₜ xₜ) t All goals completed! 🐙