daviddoji/CalculusWithJuliaNotes.jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

initial

Browse Source
This commit is contained in:
jverzani
2022-05-24 13:51:49 -04:00
parent 5c0fd1b6fe
commit 244f492f9e
240 changed files with 65211 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
CwJ/ODEs/Project.toml Normal file
View File

@@ -0,0 +1,10 @@
[deps]
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"

BIN
CwJ/ODEs/cache/euler.cache vendored Normal file
View File

Binary file not shown.

BIN
CwJ/ODEs/cache/odes.cache vendored Normal file
View File

Binary file not shown.

384
CwJ/ODEs/differential_equations.jmd Normal file
View File

@@ -0,0 +1,384 @@
# The `DifferentialEquations` suite
This section uses these add-on packages:
```julia
using OrdinaryDiffEq
using Plots
using ModelingToolkit
```
```julia; echo=false; results="hidden"
using CalculusWithJulia.WeaveSupport
const frontmatter = (
title = "The `DifferentialEquations` suite",
description = "Calculus with Julia: The `DifferentialEquations` suite",
tags = ["CalculusWithJulia", "odes", "the `differentialequations` suite"],
);
fig_size = (600, 400)
nothing
```
----
The
[`DifferentialEquations`](https://github.com/SciML/DifferentialEquations.jl)
suite of packages contains solvers for a wide range of various
differential equations. This section just briefly touches touch on
ordinary differential equations (ODEs), and so relies only on
`OrdinaryDiffEq` part of the suite. For more detail on this type and
many others covered by the suite of packages, there are many other
resources, including the
[documentation](https://diffeq.sciml.ai/stable/) and accompanying
[tutorials](https://github.com/SciML/SciMLTutorials.jl).
## SIR Model
We follow along with an introduction to the SIR model for the spread of disease by [Smith and Moore](https://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-introduction). This model received a workout due to the COVID-19 pandemic.
The basic model breaks a population into three cohorts: The **susceptible** individuals, the **infected** individuals, and the **recovered** individuals. These add to the population size, ``N``, which is fixed, but the cohort sizes vary in time. We name these cohort sizes ``S(t)``, ``I(t)``, and ``R(t)`` and define ``s(t)=S(t)/N``, ``i(t) = I(t)/N`` and ``r(t) = R(t)/N`` to be the respective proportions.
The following *assumptions* are made about these cohorts by Smith and Moore:
> No one is added to the susceptible group, since we are ignoring births and immigration. The only way an individual leaves the susceptible group is by becoming infected.
This implies the rate of change in time of ``S(t)`` depends on the current number of susceptibles, and the amount of interaction with the infected cohorts. The model *assumes* each infected person has ``b`` contacts per day that are sufficient to spread the disease. Not all contacts will be with susceptible people, but if people are assumed to mix within the cohorts, then there will be on average ``b \cdot S(t)/N`` contacts with susceptible people per infected person. As each infected person is modeled identically, the time rate of change of ``S(t)`` is:
```math
\frac{dS}{dt} = - b \cdot \frac{S(t)}{N} \cdot I(t) = -b \cdot s(t) \cdot I(t)
```
It is negative, as no one is added, only taken off. After dividing by
``N``, this can also be expressed as ``s'(t) = -b s(t) i(t)``.
> assume that a fixed fraction ``k`` of the infected group will recover during any given day.
This means the change in time of the recovered depends on ``k`` and the number infected, giving rise to the equation
```math
\frac{dR}{dt} = k \cdot I(t)
```
which can also be expressed in proportions as ``r'(t) = k \cdot i(t)``.
Finally, from ``S(t) + I(T) + R(t) = N`` we have ``S'(T) + I'(t) + R'(t) = 0`` or ``s'(t) + i'(t) + r'(t) = 0``.
Combining, it is possible to express the rate of change of the infected population through:
```math
\frac{di}{dt} = b \cdot s(t) \cdot i(t) - k \cdot i(t)
```
The author's apply this model to flu statistics from Hong Kong where:
```math
\begin{align*}
S(0) &= 7,900,000\\
I(0) &= 10\\
R(0) &= 0\\
\end{align*}
```
In `Julia` we define these, `N` to model the total population, and `u0` to be the proportions.
```julia
S0, I0, R0 = 7_900_000, 10, 0
N = S0 + I0 + R0
u0 = [S0, I0, R0]/N # initial proportions
```
An *estimated* set of values for ``k`` and ``b`` are ``k=1/3``, coming from the average period of infectiousness being estimated at three days and ``b=1/2``, which seems low in normal times, but not for an infected person who may be feeling quite ill and staying at home. (The model for COVID would certainly have a larger ``b`` value).
Okay, the mathematical modeling is done; now we try to solve for the unknown functions using `DifferentialEquations`.
To warm up, if ``b=0`` then ``i'(t) = -k \cdot i(t)`` describes the infected. (There is no circulation of people in this case.) The solution would be achieved through:
```julia; hold=true
k = 1/3
f(u,p,t) = -k * u # solving u′(t) = - k u(t)
time_span = (0.0, 20.0)
prob = ODEProblem(f, I0/N, time_span)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
plot(sol)
```
The `sol` object is a set of numbers with a convenient `plot` method. As may have been expected, this graph shows exponential decay.
A few comments are in order. The problem we want to solve is
```math
\frac{di}{dt} = -k \cdot i(t) = F(i(t), k, t)
```
where ``F`` depends on the current value (``i``), a parameter (``k``), and the time (``t``). We did not utilize ``p`` above for the parameter, as it was easy not to, but could have, and will in the following. The time variable ``t`` does not appear by itself in our equation, so only `f(u, p, t) = -k * u` was used, `u` the generic name for a solution which in this case is ``i``.
The problem we set up needs an initial value (the ``u0``) and a time span to solve over. Here we want time to model real time, so use floating point values.
The plot shows steady decay, as there is no mixing of infected with others.
Adding in the interaction requires a bit more work. We now have what is known as a *system* of equations:
```math
\begin{align*}
\frac{ds}{dt} &= -b \cdot s(t) \cdot i(t)\\
\frac{di}{dt} &= b \cdot s(t) \cdot i(t) - k \cdot i(t)\\
\frac{dr}{dt} &= k \cdot i(t)\\
\end{align*}
```
Systems of equations can be solved in a similar manner as a single ordinary differential equation, though adjustments are made to accommodate the multiple functions.
We use a style that updates values in place, and note that `u` now holds ``3`` different functions at once:
```julia
function sir!(du, u, p, t)
k, b = p
s, i, r = u[1], u[2], u[3]
ds = -b * s * i
di = b * s * i - k * i
dr = k * i
du[1], du[2], du[3] = ds, di, dr
end
```
The notation `du` is suggestive of both the derivative and a small increment. The mathematical formulation follows the derivative, the numeric solution uses a time step and increments the solution over this time step. The `Tsit5()` solver, used here, adaptively chooses a time step, `dt`; were the `Euler` method used, this time step would need to be explicit.
```julia; echo=false
note("""
The `sir!` function has the trailing `!` indicating -- by convention -- it *mutates* its first value, `du`. In this case, through an assignment, as in `du[1]=ds`. This could use some explanation. The *binding* `du` refers to the *container* holding the ``3`` values, whereas `du[1]` refers to the first value in that container. So `du[1]=ds` changes the first value, but not the *binding* of `du` to the container. That is, `du` mutates. This would be quite different were the call `du = [ds,di,dr]` which would create a new *binding* to a new container and not mutate the values in the original container.
""", title="Mutation not re-binding")
```
With the update function defined, the problem is setup and a solution found with in the same manner:
```julia;
p = (k=1/3, b=1/2) # parameters
time_span = (0.0, 150.0) # time span to solve over, 5 months
prob = ODEProblem(sir!, u0, time_span, p)
sol = solve(prob, Tsit5())
plot(sol)
plot!(x -> 0.5, linewidth=2) # mark 50% line
```
The lower graph shows the number of infected at each day over the five-month period displayed. The peak is around 6-7% of the population at any one time. However, over time the recovered part of the population reaches over 50%, meaning more than half the population is modeled as getting sick.
Now we change the parameter ``b`` and observe the difference. We passed in a value `p` holding our two parameters, so we just need to change that and run the model again:
```julia; hold=true
p = (k=1/2, b=2) # change b from 1/2 to 2 -- more daily contact
prob = ODEProblem(sir!, u0, time_span, p)
sol = solve(prob, Tsit5())
plot(sol)
```
The graphs are somewhat similar, but the steady state is reached much more quickly and nearly everyone became infected.
What about if ``k`` were bigger?
```julia; hold=true
p = (k=2/3, b=1/2)
prob = ODEProblem(sir!, u0, time_span, p)
sol = solve(prob, Tsit5())
plot(sol)
```
The graphs show that under these conditions the infections never take off; we have ``i' = (b\cdot s-k)i = k\cdot((b/k) s - 1) i`` which is always negative, since `(b/k)s < 1`, so infections will only decay.
The solution object is indexed by time, then has the `s`, `i`, `r` estimates. We use this structure below to return the estimated proportion of recovered individuals at the end of the time span.
```julia
function recovered(k,b)
prob = ODEProblem(sir!, u0, time_span, (k,b));
sol = solve(prob, Tsit5());
s,i,r = last(sol)
r
end
```
This function makes it easy to see the impact of changing the parameters. For example, fixing ``k=1/3`` we have:
```julia
f(b) = recovered(1/3, b)
plot(f, 0, 2)
```
This very clearly shows the sharp dependence on the value of ``b``; below some level, the proportion of people who are ever infected (the recovered cohort) remains near ``0``; above that level it can climb quickly towards ``1``.
The function `recovered` is of two variables returning a single value. In subsequent sections we will see a few ``3``-dimensional plots that are common for such functions, here we skip ahead and show how to visualize multiple function plots at once using "`z`" values in a graph.
```julia; hold=true
k, ks = 0.1, 0.2:0.1:0.9 # first `k` and then the rest
bs = range(0, 2, length=100)
zs = recovered.(k, bs) # find values for fixed k, each of bs
p = plot(bs, k*one.(bs), zs, legend=false) # k*one.(ks) is [k,k,...,k]
for k in ks
plot!(p, bs, k*one.(bs), recovered.(k, bs))
end
p
```
The 3-dimensional graph with `plotly` can have its viewing angle
adjusted with the mouse. When looking down on the ``x-y`` plane, which
code `b` and `k`, we can see the rapid growth along a line related to
``b/k``.
Smith and Moore point out that ``k`` is roughly the reciprocal of the number of days an individual is sick enough to infect others. This can be estimated during a breakout. However, they go on to note that there is no direct way to observe ``b``, but there is an indirect way.
The ratio ``c = b/k`` is the number of close contacts per day times the number of days infected which is the number of close contacts per infected individual.
This can be estimated from the curves once steady state has been reached (at the end of the pandemic).
```math
\frac{di}{ds} = \frac{di/dt}{ds/dt} = \frac{b \cdot s(t) \cdot i(t) - k \cdot i(t)}{-b \cdot s(t) \cdot i(t)} = -1 + \frac{1}{c \cdot s}
```
This equation does not depend on ``t``; ``s`` is the dependent variable. It could be solved numerically, but in this case affords an algebraic solution: ``i = -s + (1/c) \log(s) + q``, where ``q`` is some constant. The quantity ``q = i + s - (1/c) \log(s)`` does not depend on time, so is the same at time ``t=0`` as it is as ``t \rightarrow \infty``. At ``t=0`` we have ``s(0) \approx 1`` and ``i(0) \approx 0``, whereas ``t \rightarrow \infty``, ``i(t) \rightarrow 0`` and ``s(t)`` goes to the steady state value, which can be estimated. Solving with ``t=0``, we see ``q=0 + 1 - (1/c)\log(1) = 1``. In the limit them ``1 = 0 + s_{\infty} - (1/c)\log(s_\infty)`` or ``c = \log(s_\infty)/(1-s_\infty)``.
## Trajectory with drag
We now solve numerically the problem of a trajectory with a drag force from air resistance.
The general model is:
```math
\begin{align*}
x''(t) &= - W(t,x(t), x'(t), y(t), y'(t)) \cdot x'(t)\\
y''(t) &= -g - W(t,x(t), x'(t), y(t), y'(t)) \cdot y'(t)\\
\end{align*}
```
with initial conditions: ``x(0) = y(0) = 0`` and ``x'(0) = v_0 \cos(\theta), y'(0) = v_0 \sin(\theta)``.
This into an ODE by a standard trick. Here we define our function for updating a step. As can be seen the vector `u` contains both ``\langle x,y \rangle``
and ``\langle x',y' \rangle``
```julia
function xy!(du, u, p, t)
g, γ = p.g, p.k
x, y = u[1], u[2]
x′, y′ = u[3], u[4] # unicode \prime[tab]
W = γ
du[1] = x′
du[2] = y′
du[3] = 0 - W * x′
du[4] = -g - W * y′
end
```
This function ``W`` is just a constant above, but can be easily modified as desired.
```julia; echo=false
note("""
The "standard" trick is to take a second order ODE like ``u''(t)=u`` and turn this into two coupled ODEs by using a new name: ``v=u'(t)`` and then ``v'(t) = u(t)``. In this application, there are ``4`` equations, as we have *both* ``x''`` and ``y''`` being so converted. The first and second components of ``du`` are new variables, the third and fourth show the original equation.
""", title="A second-order ODE is a coupled first-order ODE")
```
The initial conditions are specified through:
```julia
θ = pi/4
v₀ = 200
xy₀ = [0.0, 0.0]
vxy₀ = v₀ * [cos(θ), sin(θ)]
INITIAL = vcat(xy₀, vxy₀)
```
The time span can be computed using an *upper* bound of no drag, for which the classic physics formulas give (when ``y_0=0``) ``(0, 2v_{y0}/g)``
```julia
g = 9.8
TSPAN = (0, 2*vxy₀[2] / g)
```
This allows us to define an `ODEProblem`:
```julia
trajectory_problem = ODEProblem(xy!, INITIAL, TSPAN)
```
When ``\gamma = 0`` there should be no drag and we expect to see a parabola:
```julia; hold=true
ps = (g=9.8, k=0)
SOL = solve(trajectory_problem, Tsit5(); p = ps)
plot(t -> SOL(t)[1], t -> SOL(t)[2], TSPAN...; legend=false)
```
The plot is a parametric plot of the ``x`` and ``y`` parts of the solution over the time span. We can see the expected parabolic shape.
On a *windy* day, the value of ``k`` would be positive. Repeating the above with ``k=1/4`` gives:
```julia; hold=true
ps = (g=9.8, k=1/4)
SOL = solve(trajectory_problem, Tsit5(); p = ps)
plot(t -> SOL(t)[1], t -> SOL(t)[2], TSPAN...; legend=false)
```
We see that the ``y`` values have gone negative. The `DifferentialEquations` package can adjust for that with a *callback* which terminates the problem once ``y`` has gone negative. This can be implemented as follows:
```julia; hold=true
condition(u,t,integrator) = u[2] # called when `u[2]` is negative
affect!(integrator) = terminate!(integrator) # stop the process
cb = ContinuousCallback(condition, affect!)
ps = (g=9.8, k = 1/4)
SOL = solve(trajectory_problem, Tsit5(); p = ps, callback=cb)
plot(t -> SOL(t)[1], t -> SOL(t)[2], TSPAN...; legend=false)
```
Finally, we note that the `ModelingToolkit` package provides symbolic-numeric computing. This allows the equations to be set up symbolically, as in `SymPy` before being passed off to `DifferentialEquations` to solve numerically. The above example with no wind resistance could be translated into the following:
```julia; hold=true
@parameters t γ g
@variables x(t) y(t)
D = Differential(t)
eqs = [D(D(x)) ~ -γ * D(x),
D(D(y)) ~ -g - γ * D(y)]
@named sys = ODESystem(eqs)
sys = ode_order_lowering(sys) # turn 2nd order into 1st
u0 = [D(x) => vxy₀[1],
D(y) => vxy₀[2],
x => 0.0,
y => 0.0]
p = [γ => 0.0,
g => 9.8]
prob = ODEProblem(sys, u0, TSPAN, p, jac=true)
sol = solve(prob,Tsit5())
plot(t -> sol(t)[3], t -> sol(t)[4], TSPAN..., legend=false)
```
The toolkit will automatically generate fast functions and can perform transformations (such as is done by `ode_order_lowering`) before passing along to the numeric solves.

834
CwJ/ODEs/euler.jmd Normal file
View File

@@ -0,0 +1,834 @@
# Euler's method
This section uses these add-on packages:
```julia
using CalculusWithJulia
using Plots
using SymPy
using Roots
```
```julia; echo=false; results="hidden"
using CalculusWithJulia.WeaveSupport
const frontmatter = (
title = "Euler's method",
description = "Calculus with Julia: Euler's method",
tags = ["CalculusWithJulia", "odes", "euler's method"],
);
fig_size = (600, 400)
nothing
```
----
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))
```
Plotting this solution over the slope field
```julia;
p = plot(legend=false)
vectorfieldplot!((x,y) -> [1, F(x,y)], xlims=(0, 2.5), ylims=(0, 10))
plot!(rhs(out), linewidth=5)
```
we see that the vectors that are drawn seem to be tangent to the graph
of the solution. This is no coincidence, the tangent lines to integral
curves are in the direction of the slope field.
What if the graph of the solution were not there, could we use this
fact to *approximately* reconstruct the solution?
That is, if we stitched together pieces of the slope field, would we
get a curve that was close to the actual answer?
```julia; hold=true; echo=false; cache=true
## {{{euler_graph}}}
function make_euler_graph(n)
x, y = symbols("x, y")
F(y,x) = y*x
x0, y0 = 1, 1
h = (2-1)/5
xs = zeros(n+1)
ys = zeros(n+1)
xs[1] = x0 # index is off by 1
ys[1] = y0
for i in 1:n
xs[i + 1] = xs[i] + h
ys[i + 1] = ys[i] + h * F(ys[i], xs[i])
end
p = plot(legend=false)
vectorfieldplot!((x,y) -> [1, F(y,x)], xlims=(1,2), ylims=(0,6))
## Add Euler soln
plot!(p, xs, ys, linewidth=5)
scatter!(p, xs, ys)
## add function
out = dsolve(u'(x) - F(u(x), x), u(x), ics=(u, x0, y0))
plot!(p, rhs(out), x0, xs[end], linewidth=5)
p
end
n = 5
anim = @animate for i=1:n
make_euler_graph(i)
end
imgfile = tempname() * ".gif"
gif(anim, imgfile, fps = 1)
caption = """
Illustration of a function stitching together slope field lines to
approximate the answer to an initial-value problem. The other function drawn is the actual solution.
"""
ImageFile(imgfile, caption)
```
The illustration suggests the answer is yes, let's see. The solution
is drawn over $x$ values $1$ to $2$. Let's try piecing together $5$
pieces between $1$ and $2$ and see what we have.
The slope-field vectors are *scaled* versions of the vector `[1, F(y,x)]`. The `1`
is the part in the direction of the $x$ axis, so here we would like
that to be $0.2$ (which is $(2-1)/5$. So our vectors would be `0.2 *
[1, F(y,x)]`. To allow for generality, we use `h` in place of the
specific value $0.2$.
Then our first pieces would be the line connecting $(x_0,y_0)$ to
```math
\langle x_0, y_0 \rangle + h \cdot \langle 1, F(y_0, x_0) \rangle.
```
The above uses vector notation to add the piece scaled by $h$ to the
starting point. Rather than continue with that notation, we will use
subscripts. Let $x_1$, $y_1$ be the postion of the tip of the
vector. Then we have:
```math
x_1 = x_0 + h, \quad y_1 = y_0 + h F(y_0, x_0).
```
With this notation, it is easy to see what comes next:
```math
x_2 = x_1 + h, \quad y_2 = y_1 + h F(y_1, x_1).
```
We just shifted the indices forward by $1$. But graphically what is
this? It takes the tip of the first part of our "stitched" together
solution, finds the slope filed there (`[1, F(y,x)]`) and then uses
this direction to stitch together one more piece.
Clearly, we can repeat. The $n$th piece will end at:
```math
x_{n+1} = x_n + h, \quad y_{n+1} = y_n + h F(y_n, x_n).
```
For our example, we can do some numerics. We want $h=0.2$ and $5$
pieces, so values of $y$ at $x_0=1, x_1=1.2, x_2=1.4, x_3=1.6,
x_4=1.8,$ and $x_5=2$.
Below we do this in a loop. We have to be a bit careful, as in `Julia`
the vector of zeros we create to store our answers begins indexing at
$1$, and not $0$.
```julia;
n=5
h = (2-1)/n
xs = zeros(n+1)
ys = zeros(n+1)
xs[1] = x0 # index is off by 1
ys[1] = y0
for i in 1:n
xs[i + 1] = xs[i] + h
ys[i + 1] = ys[i] + h * F(ys[i], xs[i])
end
```
So how did we do? Let's look graphically:
```julia;
plot(exp(-1/2)*exp(x^2/2), x0, 2)
plot!(xs, ys)
```
Not bad. We wouldn't expect this to be exact - due to the concavity
of the solution, each step is an underestimate. However, we see it is
an okay approximation and would likely be better with a smaller $h$. A
topic we pursue in just a bit.
Rather than type in the above command each time, we wrap it all up in
a function. The inputs are $n$, $a=x_0$, $b=x_n$, $y_0$, and, most
importantly, $F$. The output is massaged into a function through a
call to `linterp`, rather than two vectors. The `linterp` function we define below just
finds a function that linearly interpolates between the points and is
`NaN` outside of the range of the $x$ values:
```julia;
function linterp(xs, ys)
function(x)
((x < xs[1]) || (x > xs[end])) && return NaN
for i in 1:(length(xs) - 1)
if xs[i] <= x < xs[i+1]
l = (x-xs[i]) / (xs[i+1] - xs[i])
return (1-l) * ys[i] + l * ys[i+1]
end
end
ys[end]
end
end
```
With that, here is our function to find an approximate solution to $y'=F(y,x)$ with initial condition:
```julia;
function 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
ys[i + 1] = ys[i] + h * F(ys[i], xs[i])
end
linterp(xs, ys)
end
```
With `euler`, it becomes easy to explore different values.
For example, we thought the solution would look better with a smaller $h$ (or larger $n$). Instead of $n=5$, let's try $n=50$:
```julia;
u₁₂ = euler(F, 1, 2, 1, 50)
plot(exp(-1/2)*exp(x^2/2), x0, 2)
plot!(u₁₂, x0, 2)
```
It is more work for the computer, but not for us, and clearly a much better approximation to the actual answer is found.
## The Euler method
```julia; hold=true; echo=false
imgfile ="figures/euler.png"
caption = """
Figure from first publication of Euler's method. From [Gander and Wanner](http://www.unige.ch/~gander/Preprints/Ritz.pdf).
"""
ImageFile(:ODEs, imgfile, caption)
```
The name of our function reflects the [mathematician](https://en.wikipedia.org/wiki/Leonhard_Euler) associated with the iteration:
```math
x_{n+1} = x_n + h, \quad y_{n+1} = y_n + h \cdot F(y_n, x_n),
```
to approximate a solution to the first-order, ordinary differential
equation with initial values: $y'(x) = F(y,x)$.
[The Euler method](https://en.wikipedia.org/wiki/Euler_method) uses
linearization. Each "step" is just an approximation of the function
value $y(x_{n+1})$ with the value from the tangent line tangent to the
point $(x_n, y_n)$.
Each step introduces an error. The error in one step is known as the
*local truncation error* and can be shown to be about equal to $1/2
\cdot h^2 \cdot f''(x_{n})$ assuming $y$ has ``3`` or more derivatives.
The total error, or more commonly, *global truncation error*, is the
error between the actual answer and the approximate answer at the end
of the process. It reflects an accumulation of these local errors. This
error is *bounded* by a constant times $h$. Since it gets smaller as
$h$ gets smaller in direct proportion, the Euler method is called
*first order*.
Other, somewhat more complicated, methods have global truncation errors that
involve higher powers of $h$ - that is for the same size $h$, the
error is smaller. In analogy is the fact that Riemann sums have
error that depends on $h$, whereas other methods of approximating the
integral have smaller errors. For example, Simpson's rule had error
related to $h^4$. So, the Euler method may not be employed if there
is concern about total resources (time, computer, ...), it is
important for theoretical purposes in a manner similar to the role of the Riemann
integral.
In the examples, we will see that for many problems the simple Euler
method is satisfactory, but not always so. The task of numerically
solving differential equations is not a one-size-fits-all one. In the
following, a few different modifications are presented to the basic
Euler method, but this just scratches the surface of the topic.
#### Examples
##### Example
Consider the initial value problem $y'(x) = x + y(x)$ with initial
condition $y(0)=1$. This problem can be solved exactly. Here we
approximate over $[0,2]$ using Euler's method.
```julia;
𝑭(y,x) = x + y
𝒙0, 𝒙n, 𝒚0 = 0, 2, 1
𝒇 = euler(𝑭, 𝒙0, 𝒙n, 𝒚0, 25)
𝒇(𝒙n)
```
We graphically compare our approximate answer with the exact one:
```julia;
plot(𝒇, 𝒙0, 𝒙n)
𝒐ut = dsolve(D(u)(x) - 𝑭(u(x),x), u(x), ics = Dict(u(𝒙0) => 𝒚0))
plot(rhs(𝒐ut), 𝒙0, 𝒙n)
plot!(𝒇, 𝒙0, 𝒙n)
```
From the graph it appears our value for `f(xn)` will underestimate the
actual value of the solution slightly.
##### Example
The equation $y'(x) = \sin(x \cdot y)$ is not separable, so need not have an
easy solution. The default method will fail. Looking at the available methods with `sympy.classify_ode(𝐞qn, u(x))` shows a power series method which
can return a power series *approximation* (a Taylor polynomial). Let's
look at comparing an approximate answer given by the Euler method to
that one returned by `SymPy`.
First, the `SymPy` solution:
```julia;
𝐅(y,x) = sin(x*y)
𝐞qn = D(u)(x) - 𝐅(u(x), x)
𝐨ut = dsolve(𝐞qn, hint="1st_power_series")
```
If we assume $y(0) = 1$, we can continue:
```julia;
𝐨ut1 = dsolve(𝐞qn, u(x), ics=Dict(u(0) => 1), hint="1st_power_series")
```
The approximate value given by the Euler method is
```julia;
𝐱0, 𝐱n, 𝐲0 = 0, 2, 1
plot(legend=false)
vectorfieldplot!((x,y) -> [1, 𝐅(y,x)], xlims=(𝐱0, 𝐱n), ylims=(0,5))
plot!(rhs(𝐨ut1).removeO(), linewidth=5)
𝐮 = euler(𝐅, 𝐱0, 𝐱n, 𝐲0, 10)
plot!(𝐮, linewidth=5)
```
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.
##### Example
The
[Brachistochrone problem](http://www.unige.ch/~gander/Preprints/Ritz.pdf)
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
$A$, so $y(0)=0$ and at the end $y$ will be $A$.
[Galileo](http://www-history.mcs.st-and.ac.uk/HistTopics/Brachistochrone.html)
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:
```julia; hold=true; echo=false; cache=true
##{{{brach_graph}}}
function brach(f, x0, vx0, y0, vy0, dt, n)
m = 1
g = 9.8
axs = Float64[0]
ays = Float64[-g]
vxs = Float64[vx0]
vys = Float64[vy0]
xs = Float64[x0]
ys = Float64[y0]
for i in 1:n
x = xs[end]
vx = vxs[end]
ax = -f'(x) * (f''(x) * vx^2 + g) / (1 + f'(x)^2)
ay = f''(x) * vx^2 + f'(x) * ax
push!(axs, ax)
push!(ays, ay)
push!(vxs, vx + ax * dt)
push!(vys, vys[end] + ay * dt)
push!(xs, x + vxs[end] * dt)# + (1/2) * ax * dt^2)
push!(ys, ys[end] + vys[end] * dt)# + (1/2) * ay * dt^2)
end
[xs ys vxs vys axs ays]
end
fs = [x -> 1 - x,
x -> (x-1)^2,
x -> 1 - sqrt(1 - (x-1)^2),
x -> - (x-1)*(x+1),
x -> 3*(x-1)*(x-1/3)
]
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
this:
```julia;
ℱ(y,x) = -5y
𝓍0, 𝓍n, 𝓎0 = 0, 2, 1
𝓊 = euler(ℱ, 𝓍0, 𝓍n, 𝓎0, 4) # n =4 => h = 2/4
vectorfieldplot((x,y) -> [1, ℱ(y,x)], xlims=(0, 2), ylims=(-5, 5))
plot!(x -> y0 * exp(-5x), 0, 2, linewidth=5)
plot!(𝓊, 0, 2, linewidth=5)
```
What we see is that the value of $h$ is too big to capture the decay
scale of the solution. A smaller $h$, can do much better:
```julia;
𝓊₁ = euler(ℱ, 𝓍0, 𝓍n, 𝓎0, 50) # n=50 => h = 2/50
plot(x -> y0 * exp(-5x), 0, 2)
plot!(𝓊₁, 0, 2)
```
This is an example of a
[stiff equation](https://en.wikipedia.org/wiki/Stiff_equation). Such
equations cause explicit methods like the Euler one problems, as small
$h$s are needed to good results.
The implicit, backward Euler method does not have this issue, as we can see here:
```julia;
𝓊₂ = back_euler(ℱ, 𝓍0, 𝓍n, 𝓎0, 4) # n =4 => h = 2/4
vectorfieldplot((x,y) -> [1, ℱ(y,x)], xlims=(0, 2), ylims=(-1, 1))
plot!(x -> y0 * exp(-5x), 0, 2, linewidth=5)
plot!(𝓊₂, 0, 2, linewidth=5)
```
##### Example: The pendulum
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:
```math
v'(t) = -g/l \cdot \sin(u(t)), \quad u'(t) = v(t).
```
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"]
ans = 4
radioq(choices, ans, keep_order=true)
```

BIN
CwJ/ODEs/figures/bead-game.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 7.4 KiB

BIN
CwJ/ODEs/figures/euler.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 91 KiB

BIN
CwJ/ODEs/figures/galileo.gif Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 885 B

BIN
CwJ/ODEs/figures/verrazano-narrows-bridge-anniversary-historic-photos-2.jpeg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 98 KiB

914
CwJ/ODEs/odes.jmd Normal file
View File

@@ -0,0 +1,914 @@
# ODEs
This section uses these add-on packages:
```julia
using CalculusWithJulia
using Plots
using SymPy
```
```julia; echo=false; results="hidden"
using CalculusWithJulia.WeaveSupport
const frontmatter = (
title = "ODEs",
description = "Calculus with Julia: ODEs",
tags = ["CalculusWithJulia", "odes", "odes"],
);
nothing
```
----
Some relationships are easiest to describe in terms of rates or derivatives. For example:
* Knowing the speed of a car and how long it has been driving can
summarize the car's location.
* One of Newton's famous laws, $F=ma$, describes the force on an
object of mass $m$ in terms of the acceleration. The acceleration
is the derivative of velocity, which in turn is the derivative of
position. So if we know the rates of change of $v(t)$ or $x(t)$, we
can differentiate to find $F$.
* Newton's law of [cooling](http://tinyurl.com/z4lmetp). This
describes the temperature change in an object due to a difference in
temperature with the object's surroundings. The formula being,
$T'(t) = -r \left(T(t) - T_a \right)$, where $T(t)$ is temperature at time $t$
and $T_a$ the ambient temperature.
* [Hooke's law](http://tinyurl.com/kbz7r8l) relates force on an object
to the position on the object, through $F = k x$. This is
appropriate for many systems involving springs. Combined with
Newton's law $F=ma$, this leads to an equation that $x$ must
satisfy: $m x''(t) = k x(t)$.
## Motion with constant acceleration
Let's consider the case of constant acceleration. This describes how nearby objects fall to earth, as the force due to gravity is assumed to be a constant, so the acceleration is the constant force divided by the constant mass.
With constant acceleration, what is the velocity?
As mentioned, we have $dv/dt = a$ for any velocity function $v(t)$, but in this case, the right hand side is assumed to be constant. How does this restrict the possible functions, $v(t)$, that the velocity can be?
Here we can integrate to find that any answer must look like the following for some constant of integration:
```math
v(t) = \int \frac{dv}{dt} dt = \int a dt = at + C.
```
If we are given the velocity at a fixed time, say $v(t_0) = v_0$, then we can use the definite integral to get:
```math
v(t) - v(t_0) = \int_{t_0}^t a dt = at - a t_0.
```
Solving, gives:
```math
v(t) = v_0 + a (t - t_0).
```
This expresses the velocity at time $t$ in terms of the initial velocity, the constant acceleration and the time duration.
A natural question might be, is this the *only* possible answer? There are a few useful ways to think about this.
First, suppose there were another, say $u(t)$. Then define $w(t)$ to be the difference: $w(t) = v(t) - u(t)$. We would have that $w'(t) = v'(t) - u'(t) = a - a = 0$. But from the mean value theorem, a function whose derivative is *continuously* $0$, will necessarily be a constant. So at most, $v$ and $u$ will differ by a constant, but if both are equal at $t_0$, they will be equal for all $t$.
Second, since the derivative of any solution is a continuous function, it is true by the fundamental theorem of calculus that it *must* satisfy the form for the antiderivative. The initial condition makes the answer unique, as the indeterminate $C$ can take only one value.
Summarizing, we have
> If ``v(t)`` satisfies the equation: ``v'(t) = a``, ``v(t_0) = v_0,``
> then the unique solution will be ``v(t) = v_0 + a (t - t_0)``.
Next, what about position? Here we know that the time derivative of position yields the velocity, so we should have that the unknown position function satisfies this equation and initial condition:
```math
x'(t) = v(t) = v_0 + a (t - t_0), \quad x(t_0) = x_0.
```
Again, we can integrate to get an answer for any value $t$:
```math
x(t) - x(t_0) = \int_{t_0}^t \frac{dv}{dt} dt = (v_0t + \frac{1}{2}a t^2 - at_0 t) |_{t_0}^t =
(v_0 - at_0)(t - t_0) + \frac{1}{2} a (t^2 - t_0^2).
```
There are three constants: the initial value for the independent variable, $t_0$, and the two initial values for the velocity and position, $v_0, x_0$. Assuming $t_0 = 0$, we can simplify the above to get a formula familiar from introductory physics:
```math
x(t) = x_0 + v_0 t + \frac{1}{2} at^2.
```
Again, the mean value theorem can show that with the initial value specified this is the only possible solution.
## First-order initial-value problems
The two problems just looked at can be summarized by the following. We are looking for solutions to an equation of the form (taking $y$ and $x$ as the variables, in place of $x$ and $t$):
```math
y'(x) = f(x), \quad y(x_0) = y_0.
```
This is called an *ordinary differential equation* (ODE), as it is an equation involving the ordinary derivative of an unknown function, $y$.
This is called a first-order, ordinary differential equation, as there is only the first derivative involved.
This is called an initial-value problem, as the value at the initial point $x_0$ is specified as part of the problem.
#### Examples
Let's look at a few more examples, and then generalize.
##### Example: Newton's law of cooling
Consider the ordinary differential equation given by Newton's law of cooling:
```math
T'(t) = -r (T(t) - T_a), \quad T(0) = T_0
```
This equation is also first order, as it involves just the first derivative, but notice that on the right hand side is the function $T$, not the variable being differentiated against, $t$.
As we have a difference on the right hand side, we rename the variable through $U(t) = T(t) - T_a$. Then, as $U'(t) = T'(t)$, we have the equation:
```math
U'(t) = -r U(t), \quad U(0) = U_0.
```
This shows that the rate of change of $U$ depends on $U$. Large postive values indicate a negative rate of change - a push back towards the origin, and large negative values of $U$ indicate a positive rate of change - again, a push back towards the origin. We shouldn't be surprised to either see a steady decay towards the origin, or oscillations about the origin.
What will we find? This equation is different from the previous two
equations, as the function $U$ appears on both sides. However, we can
rearrange to get:
```math
\frac{dU}{dt}\frac{1}{U(t)} = -r.
```
This suggests integrating both sides, as before. Here we do the "$u$"-substitution $u = U(t)$, so $du = U'(t) dt$:
```math
-rt + C = \int \frac{dU}{dt}\frac{1}{U(t)} dt =
\int \frac{1}{u}du = \log(u).
```
Solving gives: $u = U(t) = e^C e^{-rt}$. Using the initial condition forces $e^C = U(t_0) = T(0) - T_a$ and so our solution in terms of $T(t)$ is:
```math
T(t) - T_a = (T_0 - T_a) e^{-rt}.
```
In words, the initial difference in temperature of the object and the environment exponentially decays to $0$.
That is, as $t > 0$ goes to $\infty$, the right hand will go to $0$ for $r > 0$, so $T(t) \rightarrow T_a$ - the temperature of the object will reach the ambient temperature. The rate of this is largest when the difference between $T(t)$ and $T_a$ is largest, so when objects are cooling the statement "hotter things cool faster" is appropriate.
A graph of the solution for $T_0=200$ and $T_a=72$ and $r=1/2$ is made
as follows. We've added a few line segments from the defining formula,
and see that they are indeed tangent to the solution found for the differential equation.
```julia; hold=true; echo=false
T0, Ta, r = 200, 72, 1/2
f(u, t) = -r*(u - Ta)
v(t) = Ta + (T0 - Ta) * exp(-r*t)
p = plot(v, 0, 6, linewidth=4, legend=false)
[plot!(p, x -> v(a) + f(v(a), a) * (x-a), 0, 6) for a in 1:2:5]
p
```
The above is implicitly assuming that there could be no other
solution, than the one we found. Is that really the case? We will see
that there is a theorem that can answer this, but in this case, the
trick of taking the difference of two equations satisfying the
equation leads to the equation $W'(t) = r W(t), \text{ and } W(0) =
0$. This equation has a general solution of $W(t) = Ce^{rt}$ and the
initial condition forces $C=0$, so $W(t) = 0$, as before. Hence, the
initial-value problem for Newton's law of cooling has a unique
solution.
In general, the equation could be written as (again using $y$ and $x$ as the variables):
```math
y'(x) = g(y), \quad y(x_0) = y_0
```
This is called an *autonomous*, first-order ODE, as the right-hand side does not depend on $x$ (except through ``y(x)``).
Let $F(y) = \int_{y_0}^y du/g(u)$, then a solution to the above is $F(y) = x - x_0$, assuming $1/g(u)$ is integrable.
##### Example: Toricelli's law
[Toricelli's Law](http://tinyurl.com/hxvf3qp) describes the speed a jet of water will leave a vessel through an opening below the surface of the water. The formula is $v=\sqrt{2gh}$, where $h$ is the height of the water above the hole and $g$ the gravitational constant. This arises from equating the kinetic energy gained, $1/2 mv^2$ and potential energy lost, $mgh$, for the exiting water.
An application of Torricelli's law is to describe the volume of water in a tank over time, $V(t)$. Imagine a cylinder of cross sectional area $A$ with a hole of cross sectional diameter $a$ at the bottom, Then $V(t) = A h(t)$, with $h$ giving the height. The change in volume over $\Delta t$ units of time must be given by the value $a v(t) \Delta t$, or
```math
V(t+\Delta t) - V(t) = -a v(t) \Delta t = -a\sqrt{2gh(t)}\Delta t
```
This suggests the following formula, written in terms of $h(t)$ should apply:
```math
A\frac{dh}{dt} = -a \sqrt{2gh(t)}.
```
Rearranging, this gives an equation
```math
\frac{dh}{dt} \frac{1}{\sqrt{h(t)}} = -\frac{a}{A}\sqrt{2g}.
```
Integrating both sides yields:
```math
2\sqrt{h(t)} = -\frac{a}{A}\sqrt{2g} t + C.
```
If $h(0) = h_0 = V(0)/A$, we can solve for $C = 2\sqrt{h_0}$, or
```math
\sqrt{h(t)} = \sqrt{h_0} -\frac{1}{2}\frac{a}{A}\sqrt{2g} t.
```
Setting $h(t)=0$ and solving for $t$ shows that the time to drain the tank would be $(2A)/(a\sqrt{2g})\sqrt{h_0}$.
##### Example
Consider now the equation
```math
y'(x) = y(x)^2, \quad y(x_0) = y_0.
```
This is called a *non-linear* ordinary differential equation, as the $y$ variable on the right hand side presents itself in a non-linear form (it is squared). These equations may have solutions that are not defined for all times.
This particular problem can be solved as before by moving the $y^2$ to the left hand side and integrating to yield:
```math
y(x) = - \frac{1}{C + x},
```
and with the initial condition:
```math
y(x) = \frac{y_0}{1 - y_0(x - x_0)}.
```
This answer can demonstrate *blow-up*. That is, in a finite range for $x$ values, the $y$ value can go to infinity. For example, if the initial conditions are $x_0=0$ and $y_0 = 1$, then $y(x) = 1/(1-x)$ is only defined for $x \geq x_0$ on $[0,1)$, as at $x=1$ there is a vertical asymptote.
## Separable equations
We've seen equations of the form $y'(x) = f(x)$ and $y'(x) = g(y)$ both solved by integrating. The same tricks will work for equations of the form $y'(x) = f(x) \cdot g(y)$. Such equations are called *separable*.
Basically, we equate up to constants
```math
\int \frac{dy}{g(y)} = \int f(x) dx.
```
For example, suppose we have the equation
```math
\frac{dy}{dx} = x \cdot y(x), \quad y(x_0) = y_0.
```
Then we can find a solution, $y(x)$ through:
```math
\int \frac{dy}{y} = \int x dx,
```
or
```math
\log(y) = \frac{x^2}{2} + C
```
Which yields:
```math
y(x) = e^C e^{\frac{1}{2}x^2}.
```
Substituting in $x_0$ yields a value for $C$ in terms of the initial information $y_0$ and $x_0$.
## Symbolic solutions
Differential equations are classified according to their type. Different types have different methods for solution, when a solution exists.
The first-order initial value equations we have seen can be described generally by
```math
\begin{align*}
y'(x) &= F(y,x),\\
y(x_0) &= x_0.
\end{align*}
```
Special cases include:
* *linear* if the function $F$ is linear in $y$;
* *autonomous* if $F(y,x) = G(y)$ (a function of $y$ alone);
* *separable* if $F(y,x) = G(y)H(x)$.
As seen, separable equations are approached by moving the "$y$" terms to one side, the "$x$" terms to the other and integrating. This also applies to autonomous equations then. There are other families of equation types that have exact solutions, and techniques for solution, summarized at this [Wikipedia page](http://tinyurl.com/zywzz4q).
Rather than go over these various families, we demonstrate that `SymPy` can solve many of these equations symbolically.
The `solve` function in `SymPy` solves equations for unknown
*variables*. As a differential equation involves an unknown *function*
there is a different function, `dsolve`. The basic idea is to describe
the differential equation using a symbolic function and then call
`dsolve` to solve the expression.
Symbolic functions are defined by the `@syms` macro (also see `?symbols`) using parentheses to distinguish a function from a variable:
```julia;
@syms x u() # a symbolic variable and a symbolic function
```
We will solve the following, known as the *logistic equation*:
```math
u'(x) = a u(1-u), \quad a > 0
```
Before beginning, we look at the form of the equation. When $u=0$ or
$u=1$ the rate of change is $0$, so we expect the function might be
bounded within that range. If not, when $u$ gets bigger than $1$, then
the slope is negative and when $u$ gets less than $0$, the slope is
positive, so there will at least be a drift back to the range
$[0,1]$. Let's see exactly what happens. We define a parameter,
restricting `a` to be positive:
```julia;
@syms a::positive
```
To specify a derivative of `u` in our equation we can use `diff(u(x),x)` but here, for visual simplicity, use the `Differential` operator, as follows:
```julia;
D = Differential(x)
eqn = D(u)(x) ~ a * u(x) * (1 - u(x)) # use l \Equal[tab] r, Eq(l,r), or just l - r
```
In the above, we evaluate the symbolic function at the variable `x`
through the use of `u(x)` in the expression. The equation above uses `~` to combine the left- and right-hand sides as an equation in `SymPy`. (A unicode equals is also available for this task). This is a shortcut for `Eq(l,r)`, but even just using `l - r` would suffice, as the default assumption for an equation is that it is set to `0`.
The `Differential` operation is borrowed from the `ModelingToolkit` package, which will be introduced later.
To finish, we call `dsolve` to find a solution (if possible):
```julia;
out = dsolve(eqn)
```
This answer - to a first-order equation - has one free constant,
`C_1`, which can be solved for from an initial condition. We can see
that when $a > 0$, as $x$ goes to positive infinity the solution goes
to $1$, and when $x$ goes to negative infinity, the solution goes to $0$
and otherwise is trapped in between, as expected.
The limits are confirmed by investigating the limits of the right-hand:
```julia;
limit(rhs(out), x => oo), limit(rhs(out), x => -oo)
```
We can confirm that the solution is always increasing, hence trapped within ``[0,1]`` by observing that the derivative is positive when `C₁` is positive:
```juila;
diff(rhs(out),x)
```
Suppose that $u(0) = 1/2$. Can we solve for $C_1$ symbolically? We can use `solve`, but first we will need to get the symbol for `C_1`:
```julia;
eq = rhs(out) # just the right hand side
C1 = first(setdiff(free_symbols(eq), (x,a))) # fish out constant, it is not x or a
c1 = solve(eq(x=>0) - 1//2, C1)
```
And we plug in with:
```julia;
eq(C1 => c1[1])
```
That's a lot of work. The `dsolve` function in `SymPy` allows initial conditions to be specified for some equations. In this case, ours is $x_0=0$ and $y_0=1/2$. The extra arguments passed in through a dictionary to the `ics` argument:
```julia;
x0, y0 = 0, Sym(1//2)
dsolve(eqn, u(x), ics=Dict(u(x0) => y0))
```
(The one subtlety is the need to write the rational value as a symbolic expression, as otherwise it will get converted to a floating point value prior to being passed along.)
##### Example: Hooke's law
In the first example, we solved for position, $x(t)$, from an assumption of constant acceleration in two steps. The equation relating the two is a second-order equation: $x''(t) = a$, so two constants are generated. That a second-order equation could be reduced to two first-order equations is not happy circumstance, as it can always be done. Rather than show the technique though, we demonstrate that `SymPy` can also handle some second-order ODEs.
Hooke's law relates the force on an object to its position via $F=ma = -kx$, or $x''(t) = -(k/m)x(t)$.
Suppose $k > 0$. Then we can solve, similar to the above, with:
```julia;
@syms k::positive m::positive
D2 = D ∘ D # takes second derivative through composition
eqnh = D2(u)(x) ~ -(k/m)*u(x)
dsolve(eqnh)
```
Here we find two constants, as anticipated, for we would guess that
two integrations are needed in the solution.
Suppose the spring were started by pulling it down to a bottom and
releasing. The initial position at time $0$ would be $a$, say, and
initial velocity $0$. Here we get the solution specifying initial
conditions on the function and its derivative (expressed through
`u'`):
```julia;
dsolve(eqnh, u(x), ics = Dict(u(0) => -a, D(u)(0) => 0))
```
We get that the motion will follow
$u(x) = -a \cos(\sqrt{k/m}x)$. This is simple oscillatory behavior. As the spring stretches, the force gets large enough to pull it back, and as it compresses the force gets large enough to push it back. The amplitude of this oscillation is $a$ and the period $2\pi/\sqrt{k/m}$. Larger $k$ values mean shorter periods; larger $m$ values mean longer periods.
##### Example: the pendulum
The simple gravity [pendulum](http://tinyurl.com/h8ys6ts) is an idealization of a physical pendulum that models a "bob" with mass $m$ swinging on a massless rod of length $l$ in a frictionless world governed only by the gravitational constant $g$. The motion can be described by this differential equation for the angle, $\theta$, made from the vertical:
```math
\theta''(t) + \frac{g}{l}\sin(\theta(t)) = 0
```
Can this second-order equation be solved by `SymPy`?
```julia;
@syms g::positive l::positive theta()=>"θ"
eqnp = D2(theta)(x) + g/l*sin(theta(x))
```
Trying to do so, can cause `SymPy` to hang or simply give up and repeat its input; no easy answer is forthcoming for this equation.
In general, for the first-order initial value problem characterized by
$y'(x) = F(y,x)$, there are conditions
([Peano](http://tinyurl.com/h663wba) and
[Picard-Lindelof](http://tinyurl.com/3rbde5e)) that can guarantee the
existence (and uniqueness) of equation locally, but there may not be
an accompanying method to actually find it. This particular problem
has a solution, but it can not be written in terms of elementary
functions.
However, as [Huygens](https://en.wikipedia.org/wiki/Christiaan_Huygens) first noted, if the angles involved are small, then we approximate the solution through the linearization $\sin(\theta(t)) \approx \theta(t)$. The resulting equation for an approximate answer is just that of Hooke:
```math
\theta''(t) + \frac{g}{l}\theta(t) = 0
```
Here, the solution is in terms of sines and cosines, with period given by $T = 2\pi/\sqrt{k} = 2\pi\cdot\sqrt{l/g}$. The answer does not depend on the mass, $m$, of the bob nor the amplitude of the motion, provided the small-angle approximation is valid.
If we pull the bob back an angle $a$ and release it then the initial conditions are $\theta(0) = a$ and $\theta'(a) = 0$. This gives the solution:
```julia;
eqnp₁ = D2(u)(x) + g/l * u(x)
dsolve(eqnp₁, u(x), ics=Dict(u(0) => a, D(u)(0) => 0))
```
##### Example: hanging cables
A chain hangs between two supports a distance $L$ apart. What shape
will it take if there are no forces outside of gravity acting on it?
What about if the force is uniform along length of the chain, like a
suspension bridge? How will the shape differ then?
Let $y(x)$ describe the chain at position $x$, with $0 \leq x \leq L$,
say. We consider first the case of the chain with no force save
gravity. Let $w(x)$ be the density of the chain at $x$, taken below to be a constant.
The chain is in equilibrium, so tension, $T(x)$, in the chain will be
in the direction of the derivative. Let $V$ be the vertical component
and $H$ the horizontal component. With only gravity acting on the
chain, the value of $H$ will be a constant. The value of $V$ will vary
with position.
At a point $x$, there is $s(x)$ amount of chain with weight $w \cdot s(x)$. The tension is in the direction of the tangent line, so:
```math
\tan(\theta) = y'(x) = \frac{w s(x)}{H}.
```
In terms of an increment of chain, we have:
```math
\frac{w ds}{H} = d(y'(x)).
```
That is, the ratio of the vertical and horizontal tensions in the increment are in balance with the differential of the derivative.
But $ds = \sqrt{dx^2 + dy^2} = \sqrt{dx^2 + y'(x)^2 dx^2} = \sqrt{1 + y'(x)^2}dx$, so we can simplify to:
```math
\frac{w}{H}\sqrt{1 + y'(x)^2}dx =y''(x)dx.
```
This yields the second-order equation:
```math
y''(x) = \frac{w}{H} \sqrt{1 + y'(x)^2}.
```
We enter this into `Julia`:
```julia;
@syms w::positive H::positive y()
eqnc = D2(y)(x) ~ (w/H) * sqrt(1 + y'(x)^2)
```
Unfortunately, `SymPy` needs a bit of help with this problem, by breaking the problem into
steps.
For the first step we solve for the derivative. Let $u = y'$,
then we have $u'(x) = (w/H)\sqrt{1 + u(x)^2}$:
```julia;
eqnc₁ = subs(eqnc, D(y)(x) => u(x))
```
and can solve via:
```julia;
outc = dsolve(eqnc₁)
```
So $y'(x) = u(x) = \sinh(C_1 + w \cdot x/H)$. This can be solved by direct
integration as there is no $y(x)$ term on the right hand
side.
```julia;
D(y)(x) ~ rhs(outc)
```
We see a simple linear transformation involving the hyperbolic sine. To avoid, `SymPy` struggling with the above equation, and knowing the hyperbolic sine is the derivative of the hyperbolic cosine, we anticipate an answer and verify it:
```julia;
yc = (H/w)*cosh(C1 + w*x/H)
diff(yc, x) == rhs(outc) # == not \Equal[tab]
```
The shape is a hyperbolic cosine, known as the catenary.
```julia; echo=false
imgfile = "figures/verrazano-narrows-bridge-anniversary-historic-photos-2.jpeg"
caption = """
The cables of an unloaded suspension bridge have a different shape than a loaded suspension bridge. As seen, the cables in this [figure](https://www.brownstoner.com/brooklyn-life/verrazano-narrows-bridge-anniversary-historic-photos/) would be modeled by a catenary.
"""
ImageFile(:ODEs, imgfile, caption)
```
----
If the chain has a uniform load -- like a suspension bridge with a deck -- sufficient to make the weight of the chain negligible, then how does the above change? Then the vertical tension comes from $Udx$ and not $w ds$, so the equation becomes instead:
```math
\frac{Udx}{H} = d(y'(x)).
```
This $y''(x) = U/H$, a constant. So it's answer will be a parabola.
##### Example: projectile motion in a medium
The first example describes projectile motion without air resistance. If we use $(x(t), y(t))$ to describe position at time $t$, the functions satisfy:
```math
x''(t) = 0, \quad y''(t) = -g.
```
That is, the $x$ position - where no forces act - has $0$ acceleration, and the $y$ position - where the force of gravity acts - has constant acceleration, $-g$, where $g=9.8m/s^2$ is the gravitational constant. These equations can be solved to give:
```math
x(t) = x_0 + v_0 \cos(\alpha) t, \quad y(t) = y_0 + v_0\sin(\alpha)t - \frac{1}{2}g \cdot t^2.
```
Furthermore, we can solve for $t$ from $x(t)$, to get an equation describing $y(x)$. Here are all the steps:
```julia; hold=true
@syms x0::real y0::real v0::real alpha::real g::real
@syms t x u()
a1 = dsolve(D2(u)(x) ~ 0, u(x), ics=Dict(u(0) => x0, D(u)(0) => v0 * cos(alpha)))
a2 = dsolve(D2(u)(x) ~ -g, u(x), ics=Dict(u(0) => y0, D(u)(0) => v0 * sin(alpha)))
ts = solve(t - rhs(a1), x)[1]
y = simplify(rhs(a2)(t => ts))
sympy.Poly(y, x).coeffs()
```
Though `y` is messy, it can be seen that the answer is a quadratic polynomial in $x$ yielding the familiar
parabolic motion for a trajectory. The output shows the coefficients.
In a resistive medium, there are drag forces at play. If this force is
proportional to the velocity, say, with proportion $\gamma$, then the
equations become:
```math
\begin{align*}
x''(t) &= -\gamma x'(t), & \quad y''(t) &= -\gamma y'(t) -g, \\
x(0) &= x_0, &\quad y(0) &= y_0,\\
x'(0) &= v_0\cos(\alpha),&\quad y'(0) &= v_0 \sin(\alpha).
\end{align*}
```
We now attempt to solve these.
```julia
@syms alpha::real, γ::postive, t::positive, v()
@syms x_0::real y_0::real v_0::real
Dₜ = Differential(t)
eq₁ = Dₜ(Dₜ(u))(t) ~ - γ * Dₜ(u)(t)
eq₂ = Dₜ(Dₜ(v))(t) ~ -g - γ * Dₜ(v)(t)
a₁ = dsolve(eq₁, ics=Dict(u(0) => x_0, Dₜ(u)(0) => v_0 * cos(alpha)))
a₂ = dsolve(eq₂, ics=Dict(v(0) => y_0, Dₜ(v)(0) => v_0 * sin(alpha)))
ts = solve(x - rhs(a₁), t)[1]
yᵣ = rhs(a₂)(t => ts)
```
This gives $y$ as a function of $x$.
There are a lot of symbols. Lets simplify by using constants $x_0=y_0=0$:
```julia;
yᵣ₁ = yᵣ(x_0 => 0, y_0 => 0)
```
What is the trajectory? We see
that the `log` function part will have issues when
$-\gamma x + v_0 \cos(\alpha) = 0$.
If we fix some parameters, we can plot.
```julia;
v₀, γ₀, α = 200, 1/2, pi/4
soln = yᵣ₁(v_0=>v₀, γ=>γ₀, alpha=>α, g=>9.8)
plot(soln, 0, v₀ * cos(α) / γ₀ - 1/10, legend=false)
```
We can see that the resistance makes the path quite non-symmetric.
## Visualizing a first-order initial value problem
The solution, $y(x)$, is known through its derivative. A useful tool to visualize the solution to a first-order differential equation is the [slope field](http://tinyurl.com/jspzfok) (or direction field) plot, which at different values of $(x,y)$, plots a vector with slope given through $y'(x)$.The `vectorfieldplot` of the `CalculusWithJulia` package can be used to produce these.
For example, in a previous example we found a solution to $y'(x) = x\cdot y(x)$, coded as
```julia
F(y, x) = y*x
```
Suppose $x_0=1$ and $y_0=1$. Then a direction field plot is drawn through:
```julia; hold=true
@syms x y
x0, y0 = 1, 1
plot(legend=false)
vectorfieldplot!((x,y) -> [1, F(y,x)], xlims=(x0, 2), ylims=(y0-5, y0+5))
f(x) = y0*exp(-x0^2/2) * exp(x^2/2)
plot!(f, linewidth=5)
```
In general, if the first-order equation is written as $y'(x) = F(y,x)$, then we plot a "function" that takes $(x,y)$ and returns an $x$ value of $1$ and a $y$ value of $F(y,x)$, so the slope is $F(y,x)$.
```julia; echo=false
note(L"""The order of variables in $F(y,x)$ is conventional with the equation $y'(x) = F(y(x),x)$.
""")
```
The plots are also useful for illustrating solutions for different initial conditions:
```julia; hold=true
p = plot(legend=false)
x0, y0 = 1, 1
vectorfieldplot!((x,y) -> [1,F(y,x)], xlims=(x0, 2), ylims=(y0-5, y0+5))
for y0 in -4:4
f(x) = y0*exp(-x0^2/2) * exp(x^2/2)
plot!(f, x0, 2, linewidth=5)
end
p
```
Such solutions are called [integral
curves](https://en.wikipedia.org/wiki/Integral_curve).
These graphs illustrate the fact that the slope field is tangent to the graph of any
integral curve.
## Questions
##### Question
Using `SymPy` to solve the differential equation
```math
u' = \frac{1-x}{u}
```
gives
```julia; hold=true
@syms x u()
dsolve(D(u)(x) - (1-x)/u(x))
```
The two answers track positive and negative solutions. For the initial condition, $u(-1)=1$, we have the second one is appropriate: $u(x) = \sqrt{C_1 - x^2 + 2x}$. At $-1$ this gives: $1 = \sqrt{C_1-3}$, so $C_1 = 4$.
This value is good for what values of $x$?
```julia; hold=true; echo=false
choices = [
"``[-1, \\infty)``",
"``[-1, 4]``",
"``[-1, 0]``",
"``[1-\\sqrt{5}, 1 + \\sqrt{5}]``"]
ans = 4
radioq(choices, ans)
```
##### Question
Suppose $y(x)$ satisfies
```math
y'(x) = y(x)^2, \quad y(1) = 1.
```
What is $y(3/2)$?
```julia; hold=true; echo=false
@syms x u()
out = dsolve(D(u)(x) - u(x)^2, u(x), ics=Dict(u(1) => 1))
val = N(rhs(out(3/2)))
numericq(val)
```
##### Question
Solve the initial value problem
```math
y' = 1 + x^2 + y(x)^2 + x^2 y(x)^2, \quad y(0) = 1.
```
Use your answer to find $y(1)$.
```julia; hold=true; echo=false
eqn = D(u)(x) - (1 + x^2 + u(x)^2 + x^2 * u(x)^2)
out = dsolve(eqn, u(x), ics=Dict(u(0) => 1))
val = N(rhs(out)(1).evalf())
numericq(val)
```
##### Question
A population is modeled by $y(x)$. The rate of population growth is generally proportional to the population ($k y(x)$), but as the population gets large, the rate is curtailed $(1 - y(x)/M)$.
Solve the initial value problem
```math
y'(x) = k\cdot y(x) \cdot (1 - \frac{y(x)}{M}),
```
when $k=1$, $M=100$, and $y(0) = 20$. Find the value of $y(5)$.
```julia; hold=true;echo=false
k, M = 1, 100
eqn = D(u)(x) - k * u(x) * (1 - u(x)/M)
out = dsolve(eqn, u(x), ics=Dict(u(0) => 20))
val = N(rhs(out)(5))
numericq(val)
```
##### Question
Solve the initial value problem
```math
y'(t) = \sin(t) - \frac{y(t)}{t}, \quad y(\pi) = 1
```
Find the value of the solution at $t=2\pi$.
```julia; hold=true; echo=false
eqn = D(u)(x) - (sin(x) - u(x)/x)
out = dsolve(eqn, u(x), ics=Dict(u(PI) => 1))
val = N(rhs(out(2PI)))
numericq(val)
```
##### Question
Suppose $u(x)$ satisfies:
```math
\frac{du}{dx} = e^{-x} \cdot u(x), \quad u(0) = 1.
```
Find $u(5)$ using `SymPy`.
```julia; hold=true; echo=false
eqn = D(u)(x) - exp(-x)*u(x)
out = dsolve(eqn, u(x), ics=Dict(u(0) => 1))
val = N(rhs(out)(5))
numericq(val)
```
##### Question
The differential equation with boundary values
```math
\frac{r^2 \frac{dc}{dr}}{dr} = 0, \quad c(1)=2, c(10)=1,
```
can be solved with `SymPy`. What is the value of $c(5)$?
```julia; hold=true; echo=false
@syms x u()
eqn = diff(x^2*D(u)(x), x)
out = dsolve(eqn, u(x), ics=Dict(u(1)=>2, u(10) => 1)) |> rhs
out(5) # 10/9
choices = ["``10/9``", "``3/2``", "``9/10``", "``8/9``"]
ans = 1
radioq(choices, ans)
```
##### Question
The example with projectile motion in a medium has a parameter
$\gamma$ modeling the effect of air resistance. If `y` is the
answer - as would be the case if the example were copy-and-pasted
in - what can be said about `limit(y, gamma=>0)`?
```julia; hold=true; echo=false
choices = [
"The limit is a quadratic polynomial in `x`, mirroring the first part of that example.",
"The limit does not exist, but the limit to `oo` gives a quadratic polynomial in `x`, mirroring the first part of that example.",
"The limit does not exist -- there is a singularity -- as seen by setting `gamma=0`."
]
ans = 1
radioq(choices, ans)
```

27
CwJ/ODEs/process.jl Normal file
View File

@@ -0,0 +1,27 @@
using WeavePynb
using Mustache
mmd(fname) = mmd_to_html(fname, BRAND_HREF="../toc.html", BRAND_NAME="Calculus with Julia")
## uncomment to generate just .md files
mmd(fname) = mmd_to_md(fname, BRAND_HREF="../toc.html", BRAND_NAME="Calculus with Julia")
fnames = [
"odes",
"euler"
]
function process_file(nm, twice=false)
include("$nm.jl")
mmd_to_md("$nm.mmd")
markdownToHTML("$nm.md")
twice && markdownToHTML("$nm.md")
end
process_files(twice=false) = [process_file(nm, twice) for nm in fnames]
"""
## TODO ODEs
"""

248
CwJ/ODEs/solve.jmd Normal file
View File

@@ -0,0 +1,248 @@
# The problem-algorithm-solve interface
This section uses these add-on packages:
```julia
using Plots
using MonteCarloMeasurements
```
```julia; echo=false; results="hidden"
using CalculusWithJulia.WeaveSupport
const frontmatter = (
title = "The problem-algorithm-solve interface",
description = "Calculus with Julia: The problem-algorithm-solve interface",
tags = ["CalculusWithJulia", "odes", "the problem-algorithm-solve interface"],
);
fig_size = (600, 400)
nothing
```
----
The [DifferentialEquations.jl](https://github.com/SciML) package is an entry point to a suite of `Julia` packages for numerically solving differential equations in `Julia` and other languages. A common interface is implemented that flexibly adjusts to the many different problems and algorithms covered by this suite of packages. In this section, we review a very informative [post](https://discourse.julialang.org/t/function-depending-on-the-global-variable-inside-module/64322/10) by discourse user `@genkuroki` which very nicely demonstrates the usefulness of the problem-algorithm-solve approach used with `DifferentialEquations.jl`. We slightly modify the presentation below for our needs, but suggest a perusal of the original post.
##### Example: FreeFall
The motion of an object under a uniform gravitational field is of interest.
The parameters that govern the equation of motions are the gravitational constant, `g`; the initial height, `y0`; and the initial velocity, `v0`. The time span for which a solution is sought is `tspan`.
A problem consists of these parameters. Typical `Julia` usage would be to create a structure to hold the parameters, which may be done as follows:
```julia
struct Problem{G, Y0, V0, TS}
g::G
y0::Y0
v0::V0
tspan::TS
end
Problem(;g=9.80665, y0=0.0, v0=30.0, tspan=(0.0,8.0)) = Problem(g, y0, v0, tspan)
```
The above creates a type, `Problem`, *and* a default constructor with default values. (The original uses a more sophisticated setup that allows the two things above to be combined.)
Just calling `Problem()` will create a problem suitable for the earth, passing different values for `g` would be possible for other planets.
To solve differential equations there are many different possible algorithms. Here is the construction of two types to indicate two algorithms:
```julia
struct EulerMethod{T}
dt::T
end
EulerMethod(; dt=0.1) = EulerMethod(dt)
struct ExactFormula{T}
dt::T
end
ExactFormula(; dt=0.1) = ExactFormula(dt)
```
The above just specifies a type for dispatch --- the directions indicating what code to use to solve the problem. As seen, each specifies a size for a time step with default of `0.1`.
A type for solutions is useful for different `show` methods or other methods. One can be created through:
```julia
struct Solution{Y, V, T, P<:Problem, A}
y::Y
v::V
t::T
prob::P
alg::A
end
```
The different algorithms then can be implemented as part of a generic `solve` function. Following the post we have:
```julia
solve(prob::Problem) = solve(prob, default_algorithm(prob))
default_algorithm(prob::Problem) = EulerMethod()
function solve(prob::Problem, alg::ExactFormula)
g, y0, v0, tspan = prob.g, prob.y0, prob.v0, prob.tspan
dt = alg.dt
t0, t1 = tspan
t = range(t0, t1 + dt/2; step = dt)
y(t) = y0 + v0*(t - t0) - g*(t - t0)^2/2
v(t) = v0 - g*(t - t0)
Solution(y.(t), v.(t), t, prob, alg)
end
function solve(prob::Problem, alg::EulerMethod)
g, y0, v0, tspan = prob.g, prob.y0, prob.v0, prob.tspan
dt = alg.dt
t0, t1 = tspan
t = range(t0, t1 + dt/2; step = dt)
n = length(t)
y = Vector{typeof(y0)}(undef, n)
v = Vector{typeof(v0)}(undef, n)
y[1] = y0
v[1] = v0
for i in 1:n-1
v[i+1] = v[i] - g*dt # F*h step of Euler
y[i+1] = y[i] + v[i]*dt # F*h step of Euler
end
Solution(y, v, t, prob, alg)
end
```
The post has a more elegant means to unpack the parameters from the structures, but for each of the above, the parameters are unpacked, and then the corresponding algorithm employed. As of version `v1.7` of `Julia`, the syntax `(;g,y0,v0,tspan) = prob` could also be employed.
The exact formulas, ` y(t) = y0 + v0*(t - t0) - g*(t - t0)^2/2` and `v(t) = v0 - g*(t - t0)`, follow from well-known physics formulas. Each answer is wrapped in a `Solution` type so that the answers found can be easily extracted in a uniform manner.
For example, plots of each can be obtained through:
```julia
earth = Problem()
sol_euler = solve(earth)
sol_exact = solve(earth, ExactFormula())
plot(sol_euler.t, sol_euler.y;
label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)
```
Following the post, since the time step `dt = 0.1` is not small enough, the error of the Euler method is rather large. Next we change the algorithm parameter, `dt`, to be smaller:
```julia
earth₂ = Problem()
sol_euler₂ = solve(earth₂, EulerMethod(dt = 0.01))
sol_exact₂ = solve(earth₂, ExactFormula())
plot(sol_euler₂.t, sol_euler₂.y;
label="Euler's method (dt = $(sol_euler₂.alg.dt))", ls=:auto)
plot!(sol_exact₂.t, sol_exact₂.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)
```
It is worth noting that only the first line is modified, and only the method requires modification.
Were the moon to be considered, the gravitational constant would need adjustment. This parameter is part of the problem, not the solution algorithm.
Such adjustments are made by passing different values to the `Problem`
constructor:
```julia
moon = Problem(g = 1.62, tspan = (0.0, 40.0))
sol_eulerₘ = solve(moon)
sol_exactₘ = solve(moon, ExactFormula(dt = sol_euler.alg.dt))
plot(sol_eulerₘ.t, sol_eulerₘ.y;
label="Euler's method (dt = $(sol_eulerₘ.alg.dt))", ls=:auto)
plot!(sol_exactₘ.t, sol_exactₘ.y; label="exact solution", ls=:auto)
title!("On the Moon"; xlabel="t", legend=:bottomleft)
```
The code above also adjusts the time span in addition to the
graviational constant. The algorithm for exact formula is set to use
the `dt` value used in the `euler` formula, for easier
comparison. Otherwise, outside of the labels, the patterns are the
same. Only those things that need changing are changed, the rest comes
from defaults.
The above shows the benefits of using a common interface. Next, the post illustrates how *other* authors could extend this code, simply by adding a *new* `solve` method. For example,
```julia
struct Symplectic2ndOrder{T}
dt::T
end
Symplectic2ndOrder(;dt=0.1) = Symplectic2ndOrder(dt)
function solve(prob::Problem, alg::Symplectic2ndOrder)
g, y0, v0, tspan = prob.g, prob.y0, prob.v0, prob.tspan
dt = alg.dt
t0, t1 = tspan
t = range(t0, t1 + dt/2; step = dt)
n = length(t)
y = Vector{typeof(y0)}(undef, n)
v = Vector{typeof(v0)}(undef, n)
y[1] = y0
v[1] = v0
for i in 1:n-1
ytmp = y[i] + v[i]*dt/2
v[i+1] = v[i] - g*dt
y[i+1] = ytmp + v[i+1]*dt/2
end
Solution(y, v, t, prob, alg)
end
```
Had the two prior methods been in a package, the other user could still extend the interface, as above, with just a slight standard modification.
The same approach works for this new type:
```julia
earth₃ = Problem()
sol_sympl₃ = solve(earth₃, Symplectic2ndOrder(dt = 2.0))
sol_exact₃ = solve(earth₃, ExactFormula())
plot(sol_sympl₃.t, sol_sympl₃.y; label="2nd order symplectic (dt = $(sol_sympl₃.alg.dt))", ls=:auto)
plot!(sol_exact₃.t, sol_exact₃.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)
```
Finally, the author of the post shows how the interface can compose with other packages in the `Julia` package ecosystem. This example uses the external package `MonteCarloMeasurements` which plots the behavior of the system for perturbations of the initial value:
```julia
earth₄ = Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)
sol_euler₄ = solve(earth₄)
sol_sympl₄ = solve(earth₄, Symplectic2ndOrder(dt = 2.0))
sol_exact₄ = solve(earth₄, ExactFormula())
ylim = (-100, 60)
P = plot(sol_euler₄.t, sol_euler₄.y;
label="Euler's method (dt = $(sol_euler₄.alg.dt))", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
Q = plot(sol_sympl₄.t, sol_sympl₄.y;
label="2nd order symplectic (dt = $(sol_sympl₄.alg.dt))", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
R = plot(sol_exact₄.t, sol_exact₄.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
plot(P, Q, R; size=(720, 600))
```
The only change was in the problem, `Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)`, where a different number type is used which accounts for uncertainty. The rest follows the same pattern.
This example, shows the flexibility of the problem-algorithm-solver pattern while maintaining a consistent pattern for execution.