An ordinary differential equation

This page shows a part of a tutorial by SciML (https://docs.sciml.ai/DiffEqDocs/stable/examples/classical_physics/).

In the tutorial, a simple harmonic oscillator is shown. It is given by

$$\ddot{x} + \omega^2 x = 0$$

with the analytical solution

$$\begin{eqnarray} x(t) &= A \cos (\omega t - \phi) \\ v(t) &= - A \omega \sin (\omega t - \phi) \end{eqnarray}$$

where

$$A = \sqrt{c_1 + c_2} \: \: \: \: \: \text{and} \: \: \: \: \: \tan \phi = \frac{c_2}{c_1}$$

with $c_1$, $c_2$ constants determined by the initial conditions such that $c_1$ is the initial position and $\omega c_2$ is the initial velocity.

Instead of transforming this to a system of ODEs to solve with ODEProblem, we can use SecondOrderODEProblem as follows.

using OrdinaryDiffEq, Plots

With parameter:

ω = 1;

and initial conditions:

x₀ = [0.0];

dx₀ = [π / 2];

tspan = (0.0, 2π);

ϕ = atan((dx₀[1] / ω) / x₀[1]);

A = √(x₀[1]^2 + dx₀[1]^2);

We define the problem as follows:

function harmonicoscillator(ddu, du, u, ω, t)
    ddu .= -ω^2 * u
end;

And pass it to the solvers:

prob = SecondOrderODEProblem(harmonicoscillator, dx₀, x₀, tspan, ω);

sol = solve(prob, DPRKN6())

retcode: Success
Interpolation: specialized 6th order interpolation
t: 22-element Vector{Float64}:
 0.0
 0.0006362147454811361
 0.006998362200292497
 0.037700516496439546
 0.10910842448395996
 0.2238656477827624
 0.37593069379812605
 ⋮
 4.034792459379082
 4.529739575401122
 5.069939559332106
 5.704759067328683
 6.250544837784967
 6.283185307179586
u: 22-element Vector{RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}:
 ([1.5707963267948966], [0.0])
 ([1.570796008889919], [0.0009993637178359062])
 ([1.5707578604483294], [0.010992911903844608])
 ([1.5696801498665618], [0.05920580535072108])
 ([1.561455709804332], [0.1710472642047122])
 ([1.5315995565525569], [0.34871750550843816])
 ([1.4611018566595004], [0.576699631382785])
 ⋮
 ([-0.9847669721810477], [-1.2237790077953246])
 ([-0.28531209538082897], [-1.5446677796331536])
 ([0.5497487836445709], [-1.4714542945233335])
 ([1.315265332819728], [-0.8587641916275773])
 ([1.569958816339949], [-0.05126123019995285])
 ([1.5707954668650472], [1.1681569871696004e-6])

let
    plot(
        sol,
        vars=[2, 1],
        linewidth=2,
        title="Simple Harmonic Oscillator",
        xaxis="Time",
        yaxis="Elongation",
        label=["x" "dx"],
    )
    plot!(t -> A * cos(ω * t - ϕ), lw=3, ls=:dash, label="Analytical Solution x")
    plot!(t -> -A * ω * sin(ω * t - ϕ), lw=3, ls=:dash, label="Analytical Solution dx")
end

Built with Julia 1.12.4 and

OrdinaryDiffEq 6.108.0
Plots 1.41.5

To run this tutorial locally, download this file and open it with Pluto.jl.

Dynamics

An ordinary differential equation