# Solving ODEs

## Boundary value problem

We solve the Airy equation

\[\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.*