u0, μ, κ = 1, 1., .05
u(t) = ( μ/( κ*u0*(exp(μ*t) - 1) + μ) ) * u0 * exp( μ * t )
plot(u, 0, 10,
title=L"u'=\mu u(t) - \kappa u(t)^2", legend=false,
xlabel=L"t", ylabel=L"u(t)")2 Theory of IVPs
Lecture 2: Theory of Initial Value Problems
- What are initial value problems?;
- Examples: what can happen? What can go wrong?;
- Existence and Uniqueness: when does it go right?;
- Well-posedness;
I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.
These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.
This lecture is mostly based on Burden, Faires, and Burden (2015) Section 5.1
Differential equations appear everywhere in science and engineering. Variables that change in time or space such as particles, planets, populations, …. can be modelled with differential equations. When there is one independent variable, these equations are known as Ordinary Differential Equations (ODEs). In this chapter, we study scalar initial value problems:
2.1 Examples of Initial Value Problems
Scalar initial value problems are of the form: Seek u \colon [0,T] \to \mathbb R such that
\begin{align} u'(t) &= f\big( t, u(t) \big) \nonumber\\ u(0) &= u_0 \tag{IVP} \end{align}
Remarks:
- t is the independent variable (often thought of as time),
- u_0 is the initial condition,
- u(t) can be thought of as the solution at time t,
- If f(t,u) = \alpha(t) + u \beta(t) is linear in u, we say the ODE is linear, otherwise it is nonlinear.
There is lots of different notation for derivatives:
\begin{align} u'(t) = \frac{\mathrm{d}u(t)}{\mathrm{d}t} = \frac{\mathrm{d}}{\mathrm{d}t}u(t) = \dot{u}(t) = u_t = ... \end{align}
Example 2.1 (Exponential decay/growth) Consider the following linear IVP
\begin{align} &u(0) = u_0, \\ &u'(t) = \mu u(t) \quad \text{on } (0,T). \end{align}
The rate of change of some quantity u is a constant multiple of the quanity u (e.g. growth of bacteria for small time scales, money in a fixed rate savings account, very simple model for the size of a population). One can see that there exists a unique solution given by u(t) = u_0 e^{\mu t}. This equation therefore models expoential growth or decay.
Example 2.2 (Logistic equation) Population do not grow expentially. A more realistic model of population growth includes effects limiting the size of populations: for example
\begin{align} &u(0) = u_0, \\ &u'(t) = \mu u(t) - \kappa u(t)^2 \quad \text{on } (0,T). \end{align}
One can show that there exists a unique solution to this equation given by
\begin{align} u(t) = \frac{\mu}{\kappa u_0 (e^{\mu t}-1) + \mu } u_0 e^{\mu t} \end{align}
Exercise 2.1 (Logistic equation) Let u be the solution from Example 3.2. Explain why u(t) \to \frac{\mu}{\kappa} as t\to \infty.
In general, there is no analytical formula for solutions to (\text{IVP}) and so, in practice, one needs to implement numerical methods to find solutions. For linear IVPs, one can find numerical solutions by approximating integrals:
Exercise 2.2 (Integrating Factor) For linear IVPs with f(t,u) = \alpha(t) + u \beta(t), one can rewrite the solution as an integral involving the so-called integrating factor I(t) := \exp( - \int_0^t \beta(s) \mathrm{d}s ). Show that
\begin{align} I(t) u(t) = u_0 + \int_0^t I(s) \alpha(s) \mathrm{d}s. \end{align}
In general, there is no analytical formula for this integral; one may approximate it using quadrature rules (see Math 5485).
In the following, we will use the package OrdinaryDiffEq to numerically solve some problems and see what solutions look like in practice:
using OrdinaryDiffEqExample 2.3 Consider the following linear IVP
\begin{align} &u(0) = 0, \\ &u'(t) = t^2 - u \sin(t) \quad \text{on } (0,10) \end{align}
We numerically solve this equation:
f(u,p,t) = t^2 - u * sin(t)
# p is for a set of parameters that we don't need at the moment
u0, tspan = 0., (0., 10.)
ivp = ODEProblem(f, u0, tspan)
sol = solve(ivp, Tsit5()); #Tsit5() is a particular method
plot(sol,
title=L"u'=f(t,u)", legend=false,
xlabel=L"t", ylabel=L"u(t)")
scatter!(sol.t, sol.u)The points (t_j, u(t_j)) in the plot are calculated using some approximate scheme and the graph above is an interpolation of these points.
In this case, there exists a unique solution for all time.
Exercise 2.3 Explain what happens to the solution to Example 4.1 for large t. Verify your hypothesis numerically.
tspan = (0., 1000.)
ivp = ODEProblem(f, u0, tspan)
sol = solve(ivp, Tsit5()); #Tsit5() is a particular method
plot(sol,
title=L"u'=f(t,u)", legend=false,
xlabel=L"t", ylabel=L"u(t)")Example 2.4 Consider the following (nonlinear) IVP
\begin{align} &u(0) = 1, \\ &u'(t) = u(t)^3 \quad \text{on } (0,1) \end{align}
We try to solve this problem numerically:
f(u,p,t) = u^3
u0, tspan = 1., (0., 1.)
ivp = ODEProblem(f, u0, tspan)
sol = solve(ivp, Tsit5())
plot(sol, yaxis=:log,
title=L"u'=u^3", legend=false,
xlabel=L"t", ylabel=L"u(t)")┌ Warning: At t=0.4999988132731494, dt was forced below floating point epsilon 5.551115123125783e-17, and step error estimate = 6311.956553629581. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64). └ @ SciMLBase C:\Users\math5\.julia\packages\SciMLBase\TZ9Rx\src\integrator_interface.jl:671
Notice that the solution appears to blow-up in finite time! This is not an artifact of the numerical scheme and, in fact, the solution that exists for small t cannot be extended past T = 0.5.
Exercise 2.4 Show that the solution to Example 4.2 is
\begin{align} u(t) = \frac{1}{\sqrt{1-2t}} \end{align}
which does indeed have a singularity at T = \frac12.
Example 2.5 Consider the following IVP
\begin{align} &u(0) = 0, \\ &u'(t) = \sqrt{u(t)} \quad \text{on } (0,1) \end{align}
We try to solve this problem numerically:
f(u,p,t) = sqrt(u)
u0, tspan = 0., (0., 1.)
ivp = ODEProblem(f, u0, tspan)
sol = solve(ivp, Tsit5())
plot(sol,
title=L"u'=√u", legend=false,
xlabel=L"t", ylabel=L"u(t)", ylims=(-1,1))Notice that u(t) = 0 is indeed a solution to this differential equation. However, this solution is not unique.
Exercise 2.5 Show that \tilde{u}(t) := \frac14 t^{2} satisfies the IVP in Example 2.5. (Notice that the numerical scheme implemented above did not find this solution).
In these examples, we have seen IVPs that (i) has a unique solution for all time, (ii) has a unique solution for some time interval [0,T) that cannot be extended to [0,T] (finite time blow-up), and (iii) has non-unique solutions. It would be nice if we had some analytical results that tell us if/when there exists a unique solution, and over what kind of time interval; this is the content of the next section.
2.2 Existence and Uniqueness Theory
Pages 259–262 of Burden, Faires, and Burden (2015).
Theorem 2.1 (Picard–Lindelöf) Let D := \{ (t,u) : t \in [0,T], u \in \mathbb R \}. Suppose that f is continuous on D and Lipschitz in its second argument (with constant independent of t\in [0,T]). Then, (\text{IVP}) has a unique solution u : [0,T] \to \mathbb R.
Ideas in the Proof (non-examinable). Define
\begin{align} \mathcal Pu(t) := u_0(t) + \int_0^t f\big(s,u(s)\big) \mathrm{d}s. \end{align}
Notice that if u solves (\text{IVP}), then u = \mathcal Pu. The hypotheses on f imply that \mathcal P is a contraction and so we can apply contraction mapping theorem (a more general version of the result we saw in Math 5485) to conclude that there exists a unique fixed point u = \mathcal P u and the iteration u_{n+1} := \mathcal P u_n converges to u. This is known as Picard’s iteration and (u_n) is the sequence of Picard iterates.
Exercise 2.6 What do you notice about the Picard iterates corresponding to the problem
\begin{align} &u(0) = 1 \\ &u'(t) = u(t). \end{align}
Exercise 2.7 Recall that in Example 4.1, we saw that the IVP
\begin{align} &u(0) = 0, \\ &u'(t) = t^2 - u \sin(t) \quad \text{on } (0,10) \end{align}
has a (numerical) solution. Show that this solution is unique and exists for all T \geq 0.
Exercise 2.8 Recall that in Example 4.2, we saw that the IVP
\begin{align} &u(0) = 1, \\ &u'(t) = u^3 \end{align}
has finite-time blow up at T=\frac{1}{2}. Explain why this does not contradict the existence and uniqueness theorem we have just seen.
Exercise 2.9 Recall that in Example 2.5, we saw that the IVP
\begin{align} &u(0) = 0, \\ &u'(t) = \sqrt{u} \quad \text{on } (0,1) \end{align}
has non-unique solutions. Explain why this does not contradict the existence and uniqueness theorem.
2.3 Well-posedness
We will use the notation \| g \|_{L^\infty([0,T])} := \max_{t \in [0,T]} |g(t)|.
In practice, there is some error in the data (u_0 and f). Moreover, when we implement numerical schemes, we introduce round-off errors (coming from the fact we are working in floating point arithmetic, see Math 5485). We therefore wish to know that the true solution is continous in the data: i.e. small changes in the data leads only to small changes in the output.
Definition 2.1 We say (\text{IVP}) is well-posed if (i) there exists a unique solution u and (ii) the solution is continuous with respect to small changes in the data: there exists \epsilon > 0 and C > 0 such that if \delta_0 > 0 and \delta:[0,T]\to \mathbb R is continuous with
- | \delta_0 | \leq \epsilon,
- \| \delta \|_{L^\infty([0,T])} \leq \epsilon,
then, the IVP
\begin{align} &v(0) = u_0 + \delta_0 \nonumber\\ &v'(t) = f\big( t, v(t) \big) + \delta(t) \end{align}
has a unique solution v : [0,T] \to \mathbb R satisfying
\begin{align} \| u - v \|_{L^\infty([0,T])} \leq C \epsilon. \end{align}
This is called a stability estimate.
Theorem 2.2 Under the assumptions of Theorem 2.1, (\text{IVP}) is well-posed.
Suppose that \text{(IVP)} satisfies the conditions of Theorem 2.1 and we only perturb the initial condition u_0:
\begin{align} &v(0) = u_0 + \delta \nonumber\\ &v'(t) = f\big(t,v(t)\big) \end{align}
Then, \| u - v \|_{L^\infty([0,T])} \leq |\delta| e^{LT} for all sufficiently small |\delta| > 0.
This bound is often very pessimistic. However:
Example 2.6 Consider again the example
\begin{align} &u(0) = u_0 \nonumber\\ &u'(t) = \mu u(t). \end{align}
Notice, in this case, f(t,u) = \mu u.
(i) Show that f is Lipschitz with constant L,
(ii) Write down the stability estimate from Note 2.1,
(iii) Show that you cannot improve on this estimate.
We look again at Example 4.1, perturb the initial condition, and plot the solutions to \text{(IVP)} with u(0) = 0 and v(0)= \delta:
δ = 10
f(u,p,t) = t^2 - u * sin(t)
L = 1.0
u0, tspan = 0., (0., 10.)
e(t) = δ * exp( L*t )
ivp = ODEProblem(f, u0, tspan)
sol = solve(ivp, Tsit5());
ivp2 = ODEProblem(f, u0+δ, tspan)
sol2 = solve(ivp2, Tsit5());
plot(sol,
title=L"u'=f(t,u)", legend=false,
xlabel=L"t", ylabel=L"u(t)")
plot!(sol2)
# plot!( t->sol(t)+e(t), linestyle=:dash )
# plot!( t->sol(t)-e(t), linestyle=:dash )You can see that the bound \| u - v \|_{L^\infty([0,T])} \leq |\delta| e^{LT} is pessimistic in this case.
- (If you haven’t already) Read pages 259–266 of Burden, Faires, and Burden (2015) - please come to the lecture on Wednesday with any questions you have,
- A1 due on Wednesday,
- (Sign-up to office hours if you would like)
Next: Euler’s method