GitHub

An ordinary differential equation

This page shows a part of a tutorial by SciML (https://tutorials.sciml.ai/html/models/01-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())
timestampvalue1value2
10.01.57080.0
20.0006362151.57080.000999364
30.006998361.570760.0109929
40.03770051.569680.0592058
50.1091081.561460.171047
60.2238661.53160.348718
70.3759311.46110.5767
80.5656371.326140.841874
90.7978091.096851.12442
101.074870.7474551.38156
...
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.9.1 and

OrdinaryDiffEq 6.53.1
Plots 1.38.16

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