Solving ODEs

Boundary value problem

\[\frac{d^2 y}{dx^2}-x y = 0,\]

in an interval a..b, subject to the boundary conditions $y(a)=\mathrm{Ai}(a)$ and $y(b)=\mathrm{Ai}(b)$, where $\mathrm{Ai}$ represents the Airy function of the first kind.

using ApproxFun
using SpecialFunctions

Construct the domain

a, b = -20, 10
d = a..b;

We construct the Airy differential operator $L = d^2/dx^2 - x$:

x = Fun(d);
D = Derivative(d);
L = D^2 - x;

We impose boundary conditions by setting the function at the boundaries equal to the values of the Airy function at these points. First, we construct the Dirichlet boundary condition operator, that evaluates the function at the boundaries

B = Dirichlet(d);

The right hand side for the boundary condition is obtained by evaluating the Airy function We use airyai from SpecialFunctions.jl for this

B_vals = [airyai(a), airyai(b)];

The main step, where we solve the differential equation

u = [B; L] \ [B_vals, 0];

We plot the solution

import Plots
Plots.plot(u, xlabel="x", ylabel="u(x)", legend=false)

Initial value problem

Undamped harmonic oscillator

\[\frac{d^2 u}{dt^2} + 4 u = 0,\]

in an interval 0..2pi, subject to the initial conditions $u(0)=0$ and $u'(0)=2$.

using ApproxFun
using LinearAlgebra

Construct the domain

a, b = 0, 2pi
d = a..b;

We construct the differential operator $L = d^2/dx^2 + 4$:

D = Derivative(d)
L = D^2 + 4;

We have not chosen the space explicitly, and the solve chooses it to be Chebyshev(d) internally

We incorporate the initial conditions using ivp, which is a shortcut for evaluating the function and it's derivative at the left boundary of the domain

A = [L; ivp()];

Initial conditions

u0 = 0;
dtu0 = 2;

We solve the differential equation

u = \(A, [0,u0,dtu0], tolerance=1e-6);

We plot the solution

t = a:0.1:b
Plots.plot(t, u.(t), xlabel="t", label="u(t)")
Plots.plot!(t, sin.(2t), label="Analytical", seriestype=:scatter)

Damped harmonic oscillator

\[\frac{d^2 y}{dt^2} + 2\zeta\omega_0\frac{dy}{dt} + \omega_0^2 y = 0,\]

from $t=0$ to $t=T$, subject to the initial conditions $y(0)=4$ and $y'(0)=0$.

Construct the domain

T = 12
d = 0..T;

Dt = Derivative(d);
ζ = 0.2 # damping ratio
ω0 = 2 # oscillation frequency
L = Dt^2 + 2ζ * ω0 * Dt + ω0^2;

initial conditions

y0 = 4; # initial displacement
dty0 = 3; # initial velocity

The differential operator along with the initial condition evaluation operator

A = [L; ivp()];

We solve the differential equation

y = \(A, [0,y0,dty0], tolerance=1e-6);

Plot the solution

import Plots
Plots.plot(y, xlabel="t", ylabel="y(t)", legend=false)

Increasing Precision

Solving differential equations with high precision types is available. The following calculates $e$ to 300 digits by solving the ODE $u^\prime = u$:

using ApproxFun
u = setprecision(1000) do
    d = BigFloat(0)..BigFloat(1)
    D = Derivative(d)
    [ldirichlet(); D-1] \ [1; 0]
end;

import Plots
Plots.plot(u; legend=false, xlabel="x", ylabel="u(x)")

This page was generated using Literate.jl.