The following section takes up the task of numerically approximating solutions to differential equations. `Julia` has a huge set of state-of-the-art tools for this task starting with the [DifferentialEquations](https://github.com/SciML/DifferentialEquations.jl) package. We don't use that package in this section, focusing on simpler methods and implementations for pedagogical purposes, but any further exploration should utilize the tools provided therein. A brief introduction to the package follows in an upcoming [section](./differential_equations.html).
----
Consider the differential equation:
```math
y'(x) = y(x) \cdot x, \quad y(1)=1,
```
which can be solved with `SymPy`:
```julia;
@syms x, y, u()
D = Differential(x)
x0, y0 = 1, 1
F(y,x) = y*x
dsolve(D(u)(x) - F(u(x), x))
```
With the given initial condition, the solution becomes:
```julia;
out = dsolve(D(u)(x) - F(u(x),x), u(x), ics=Dict(u(x0) => y0))
We see that the answer found from using a polynomial series matches that of Euler's method for a bit, but as time evolves, the approximate solution given by Euler's method more closely tracks the slope field.
was posed by Johann Bernoulli in 1696. It asked for the curve between
two points for which an object will fall faster along that curve than
any other. For an example, a bead sliding on a wire will take a certain amount of time to get from point $A$ to point $B$, the time depending on the shape of the wire. Which shape will take the least amount of time?
```julia; hold=true; echo=false
imgfile = "figures/bead-game.jpg"
caption = """
A child's bead game. What shape wire will produce the shortest time for a bed to slide from a top to the bottom?
"""
ImageFile(:ODEs, imgfile, caption)
```
Restrict our attention to the $x$-$y$ plane, and consider a path,
between the point $(0,A)$ and $(B,0)$. Let $y(x)$ be the distance from
knew the straight line was not the curve, but incorrectly thought the
answer was a part of a circle.
```julia; hold=true; echo=false
imgfile = "figures/galileo.gif"
caption = """
As early as 1638, Galileo showed that an object falling along `AC` and then `CB` will fall faster than one traveling along `AB`, where `C` is on the arc of a circle.
From the [History of Math Archive](http://www-history.mcs.st-and.ac.uk/HistTopics/Brachistochrone.html).
"""
ImageFile(:ODEs, imgfile, caption)
```
This simulation also suggests that a curved path is better than the shorter straight one:
MS = [brach(f, 1/100, 0, 1, 0, 1/100, 100) for f in fs]
function make_brach_graph(n)
p = plot(xlim=(0,1), ylim=(-1/3, 1), legend=false)
for (i,f) in enumerate(fs)
plot!(f, 0, 1)
U = MS[i]
x = min(1.0, U[n,1])
scatter!(p, [x], [f(x)])
end
p
end
n = 4
anim = @animate for i=[1,5,10,15,20,25,30,35,40,45,50,55,60]
make_brach_graph(i)
end
imgfile = tempname() * ".gif"
gif(anim, imgfile, fps = 1)
caption = """
The race is on. An illustration of beads falling along a path, as can be seen, some paths are faster than others. The fastest path would follow a cycloid. See [Bensky and Moelter](https://pdfs.semanticscholar.org/66c1/4d8da6f2f5f2b93faf4deb77aafc7febb43a.pdf) for details on simulating a bead on a wire.
"""
ImageFile(imgfile, caption)
```
Now, the natural question is which path is best? The solution can be
[reduced](http://mathworld.wolfram.com/BrachistochroneProblem.html) to
solving this equation for a positive $c$:
```math
1 + (y'(x))^2 = \frac{c}{y}, \quad c > 0.
```
Reexpressing, this becomes:
```math
\frac{dy}{dx} = \sqrt{\frac{C-y}{y}}.
```
This is a separable equation and can be solved, but even `SymPy` has
trouble with this integral. However, the result has been known to be a piece of a cycloid since the insightful
Jacob Bernoulli used an analogy from light bending to approach the problem. The answer is best described parametrically
through:
```math
x(u) = c\cdot u - \frac{c}{2}\sin(2u), \quad y(u) = \frac{c}{2}( 1- \cos(2u)), \quad 0 \leq u \leq U.
```
The values of $U$ and $c$ must satisfy $(x(U), y(U)) = (B, A)$.
Rather than pursue this, we will solve it numerically for a fixed
value of $C$ over a fixed interval to see the shape.
The equation can be written in terms of $y'=F(y,x)$, where
```math
F(y,x) = \sqrt{\frac{c-y}{y}}.
```
But as $y_0 = 0$, we immediately would have a problem with the first step, as there would be division by $0$.
This says that for the optimal solution, the bead picks up speed by first sliding straight down before heading off towards $B$. That's great for the physics, but runs roughshod over our Euler method, as the first step has an infinity.
For this, we can try the *backwards Euler* method which uses the slope at $(x_{n+1}, y_{n+1})$, rather than $(x_n, y_n)$. The update step becomes:
```math
y_{n+1} = y_n + h \cdot F(y_{n+1}, x_{n+1}).
```
Seems innocuous, but the value we are trying to find, $y_{n+1}$, is
now on both sides of the equation, so is only *implicitly* defined. In
this code, we use the `find_zero` function from the `Roots` package. The
caveat is, this function needs a good initial guess, and the one we
use below need not be widely applicable.
```julia;
function back_euler(F, x0, xn, y0, n)
h = (xn - x0)/n
xs = zeros(n+1)
ys = zeros(n+1)
xs[1] = x0
ys[1] = y0
for i in 1:n
xs[i + 1] = xs[i] + h
## solve y[i+1] = y[i] + h * F(y[i+1], x[i+1])
ys[i + 1] = find_zero(y -> ys[i] + h * F(y, xs[i + 1]) - y, ys[i]+h)
end
linterp(xs, ys)
end
```
We then have with $C=1$ over the interval $[0,1.2]$ the following:
```julia;
𝐹(y, x; C=1) = sqrt(C/y - 1)
𝑥0, 𝑥n, 𝑦0 = 0, 1.2, 0
cyc = back_euler(𝐹, 𝑥0, 𝑥n, 𝑦0, 50)
plot(x -> 1 - cyc(x), 𝑥0, 𝑥n)
```
Remember, $y$ is the displacement from the top, so it is
non-negative. Above we flipped the graph to make it look more like
expectation. In general, the trajectory may actually dip below the
ending point and come back up. The above won't see this, for as
written $dy/dx \geq 0$, which need not be the case, as the defining
equation is in terms of $(dy/dx)^2$, so the derivative could have any
sign.
##### Example: stiff equations
The Euler method is *convergent*, in that as $h$ goes to $0$, the
approximate solution will converge to the actual answer. However, this
does not say that for a fixed size $h$, the approximate value will be
good. For example, consider the differential equation $y'(x) =
-5y$. This has solution $y(x)=y_0 e^{-5x}$. However, if we try the
Euler method to get an answer over $[0,2]$ with $h=0.5$ we don't see
The differential equation describing the simple pendulum is
```math
\theta''(t) = - \frac{g}{l}\sin(\theta(t)).
```
The typical approach to solving for $\theta(t)$ is to use the small-angle approximation that $\sin(x) \approx x$, and then the differential equation simplifies to:
$\theta''(t) = -g/l \cdot \theta(t)$, which is easily solved.
Here we try to get an answer numerically. However, the problem, as stated, is not a first order equation due to the $\theta''(t)$ term. If we let $u(t) = \theta(t)$ and $v(t) = \theta'(t)$, then we get *two* coupled first order equations:
We can try the Euler method here. A simple approach might be this iteration scheme:
```math
\begin{align*}
x_{n+1} &= x_n + h,\\
u_{n+1} &= u_n + h v_n,\\
v_{n+1} &= v_n - h \cdot g/l \cdot \sin(u_n).
\end{align*}
```
Here we need *two* initial conditions: one for the initial value
$u(t_0)$ and the initial value of $u'(t_0)$. We have seen if we start at an angle $a$ and release the bob from rest, so $u'(0)=0$ we get a sinusoidal answer to the linearized model. What happens here? We let $a=1$, $L=5$ and $g=9.8$:
We write a function to solve this starting from $(x_0, y_0)$ and ending at $x_n$:
```julia;
function euler2(x0, xn, y0, yp0, n; g=9.8, l = 5)
xs, us, vs = zeros(n+1), zeros(n+1), zeros(n+1)
xs[1], us[1], vs[1] = x0, y0, yp0
h = (xn - x0)/n
for i = 1:n
xs[i+1] = xs[i] + h
us[i+1] = us[i] + h * vs[i]
vs[i+1] = vs[i] + h * (-g / l) * sin(us[i])
end
linterp(xs, us)
end
```
Let's take $a = \pi/4$ as the initial angle, then the approximate
solution should be $\pi/4\cos(\sqrt{g/l}x)$ with period $T =
2\pi\sqrt{l/g}$. We try first to plot then over 4 periods:
```julia;
𝗅, 𝗀 = 5, 9.8
𝖳 = 2pi * sqrt(𝗅/𝗀)
𝗑0, 𝗑n, 𝗒0, 𝗒p0 = 0, 4𝖳, pi/4, 0
plot(euler2(𝗑0, 𝗑n, 𝗒0, 𝗒p0, 20), 0, 4𝖳)
```
Something looks terribly amiss. The issue is the step size, $h$, is
too large to capture the oscillations. There are basically only $5$
steps to capture a full up and down motion. Instead, we try to get $20$ steps per period
so $n$ must be not $20$, but $4 \cdot 20 \cdot T \approx 360$. To this
graph, we add the approximate one:
```julia;
plot(euler2(𝗑0, 𝗑n, 𝗒0, 𝗒p0, 360), 0, 4𝖳)
plot!(x -> pi/4*cos(sqrt(𝗀/𝗅)*x), 0, 4𝖳)
```
Even now, we still see that something seems amiss, though the issue is
not as dramatic as before. The oscillatory nature of the pendulum is
seen, but in the Euler solution, the amplitude grows, which would
necessarily mean energy is being put into the system. A familiar
instance of a pendulum would be a child on a swing. Without pumping
the legs - putting energy in the system - the height of the swing's
arc will not grow. Though we now have oscillatory motion, this growth
indicates the solution is still not quite right. The issue is likely
due to each step mildly overcorrecting and resulting in an overall
growth. One of the questions pursues this a bit further.
## Questions
##### Question
Use Euler's method with $n=5$ to approximate $u(1)$ where
```math
u'(x) = x - u(x), \quad u(0) = 1
```
```julia; hold=true; echo=false
F(y,x) = x - y
x0, xn, y0 = 0, 1, 1
val = euler(F, x0, xn, y0, 5)(1)
numericq(val)
```
##### Question
Consider the equation
```math
y' = x \cdot \sin(y), \quad y(0) = 1.
```
Use Euler's method with $n=50$ to find the value of $y(5)$.
```julia; hold=true; echo=false
F(y, x) = x * sin(y)
x0, xn, y0 = 0, 5, 1
n = 50
u = euler(F, x0, xn, y0, n)
numericq(u(xn))
```
##### Question
Consider the ordinary differential equation
```math
\frac{dy}{dx} = 1 - 2\frac{y}{x}, \quad y(1) = 0.
```
Use Euler's method to solve for $y(2)$ when $n=50$.
```julia; hold=true; echo=false
F(y, x) = 1 - 2y/x
x0, xn, y0 = 1, 2, 0
n = 50
u = euler(F, x0, xn, y0, n)
numericq(u(xn))
```
##### Question
Consider the ordinary differential equation
```math
\frac{dy}{dx} = \frac{y \cdot \log(y)}{x}, \quad y(2) = e.
```
Use Euler's method to solve for $y(3)$ when $n=25$.
```julia; hold=true; echo=false
F(y, x) = y*log(y)/x
x0, xn, y0 = 2, 3, exp(1)
n = 25
u = euler(F, x0, xn, y0, n)
numericq(u(xn))
```
##### Question
Consider the first-order non-linear ODE
```math
y' = y \cdot (1-2x), \quad y(0) = 1.
```
Use Euler's method with $n=50$ to approximate the solution $y$ over $[0,2]$.
What is the value at $x=1/2$?
```julia; hold=true; echo=false
F(y, x) = y * (1-2x)
x0, xn, y0 = 0, 2, 1
n = 50
u = euler(F, x0, xn, y0, n)
numericq(u(1/2))
```
What is the value at $x=3/2$?
```julia; hold=true; echo=false
F(y, x) = y * (1-2x)
x0, xn, y0 = 0, 2, 1
n = 50
u = euler(F, x0, xn, y0, n)
numericq(u(3/2))
```
##### Question: The pendulum revisited.
The issue with the pendulum's solution growing in amplitude can be
addressed using a modification to the Euler method attributed to
[Cromer](http://astro.physics.ncsu.edu/urca/course_files/Lesson14/index.html). The
fix is to replace the term `sin(us[i])` in the line `vs[i+1] = vs[i] + h * (-g / l) *
sin(us[i])` of the `euler2` function with `sin(us[i+1])`, which uses the updated angular
velocity in the ``2``nd step in place of the value before the step.
Modify the `euler2` function to implement the Euler-Cromer method. What do you see?
```julia; hold=true; echo=false
choices = [
"The same as before - the amplitude grows",
"The solution is identical to that of the approximation found by linearization of the sine term",
"The solution has a constant amplitude, but its period is slightly *shorter* than that of the approximate solution found by linearization",
"The solution has a constant amplitude, but its period is slightly *longer* than that of the approximate solution found by linearization"]
