Basics of IVPs

6.1. Basics of IVPs#

Definition 6.1.1 : Initial-value problem (scalar)

A scalar first-order initial-value problem (IVP) is

(6.1.1)#\[\begin{split}\begin{split} u'(t) &= f(t,u(t)), \qquad a \le t \le b, \\ u(a) &=u_0. \end{split}\end{split}\]

We call \(t\) the independent variable and \(u\) the dependent variable. If \(u'=f(t,u)=g(t)+u h(t)\), the differential equation is linear; otherwise, it is nonlinear.

A solution of an initial-value problem is a function \(u(t)\) that makes both \(u'(t)=f\bigl(t,u(t)\bigr)\) and \(u(a)=u_0\) true equations.

When \(t\) is meant to be time, sometimes we write \(\dot{u}\) (read “u-dot”) instead of \(u'\).

Example 6.1.2

Suppose \(u(t)\) is the size of a population at time \(t\). We idealize by allowing \(u\) to take any real (not just integer) value. If we assume a constant per capita birth rate (births per unit population per unit time), then

\[\frac{d u}{d t} = k u, \qquad u(0)=u_0\]

for some \(k>0\). The solution of this linear equation is \(u(t)=e^{kt}u_0\), which is exponential growth.

A more realistic model would cap the growth due to finite resources. Suppose the death rate is proportional to the size of the population, indicating competition. Then

(6.1.2)#\[ \frac{d u}{d t} = ku - ru^2, \qquad u(0)=u_0.\]

This is the logistic equation. Although crude, it is still useful in population models. The solution relevant for population models has the form

\[ u(t) = \frac{k/r}{ 1 + \left( \frac{k}{r u_0} - 1 \right) e^{-k t} }.\]

For \(k,r,u_0>0\), the solution smoothly varies from the initial population \(u_0\) to a finite population, equal to \(k/r\), that has been limited by competition.

Linear problems can be solved in terms of integrals. Defining the integrating factor \(\rho(t) = \exp\bigl[\int -h(t)\, dt \bigr]\), the solution is derived from

\[ \rho(t) u(t) = u_0 + \int_a^t \rho(s) g(s) \, ds.\]

In many cases, however, the necessary integrals cannot be done in closed form. Some nonlinear ODEs, such as separable equations, may also be solvable with a short formula, perhaps with difficult integrations. Most often, though, there is no analytic formula available for the solution.

An ODE may have higher derivatives of the unknown solution present. For example, a second-order ordinary differential equation is often given in the form \(u''(t)=f\bigl(t,u,u'\bigr)\). A second-order IVP requires two conditions at the initial time in order to specify a solution completely. As we will see in Section 6.3, we are always able to reformulate higher-order IVPs in a first-order form, so we will deal with first-order problems exclusively.

Numerical solutions#

Demo 6.1.3

The DifferentialEquations package offers solvers for IVPs. Let’s use it to define and solve an initial-value problem for \(u'=\sin[(u+t)^2]\) over \(t \in [0,4]\), such that \(u(0)=-1\).

Because many practical problems come with parameters that are fixed within an instance but varied from one instance to another, the syntax for IVPs includes a input argument p that stays fixed throughout the solution. Here we don’t want to use that argument, but it must be in the definition for the solver to work.

To create an initial-value problem for \(u(t)\), you must supply a function that computes \(u'\), an initial value for \(u\), and the endpoints of the interval for \(t\). The \(t\) interval should be defined as (a,b), where at least one of the values is a float.

f = (u,p,t) -> sin((t+u)^2)     # defines du/dt, must include p argument
u₀ = -1.0                       # initial value
tspan = (0.0,4.0)               # t interval 

(0.0, 4.0)

With the data above we define an IVP problem object and then solve it. Here we tell the solver to use the Tsit5 method, which is a good first choice for most problems.

ivp = ODEProblem(f,u₀,tspan)
sol = solve(ivp,Tsit5());

The resulting solution object can be shown using plot.

plot(sol,label="solution",legend=:bottom,
    xlabel="t",ylabel=L"u(t)",title=L"u'=\sin((t+u)^2)")

The solution also acts like any callable function that can be evaluated at different values of \(t\).

@show sol(1.0);

sol(1.0) = -0.7903205813665345

Under the hood, the solution object holds some information about how the values and plot are produced:

[sol.t sol.u]

15×2 Matrix{Float64}:
0        -1.0
0867807  -0.93483
241035   -0.856617
464665   -0.805668
696832   -0.793614
00862    -0.789925
37461    -0.718601
70407    -0.476837
93572    -0.29033
17184    -0.294994
4843     -0.483948
69425    -0.654121
27049    -1.1783
62534    -1.51729
0        -1.88086

The solver initially finds approximate values of the solution (second column above) at some automatically chosen times (first column above). To compute the solution at other times, the object performs an interpolation on those values. This chapter is about how the discrete \(t\) and \(u\) values are computed. For now, just note how we can extract them from the solution object.

scatter!(sol.t,sol.u,label="discrete values")

Existence and uniqueness#

There are simple IVPs that do not have solutions at all possible times.

Demo 6.1.4

The equation \(u'=(u+t)^2\) gives us some trouble.

f = (u,p,t) -> (t+u)^2

ivp = ODEProblem(f,1.0,(0.,1.))
sol = solve(ivp,Tsit5());

┌ Warning: At t=0.7853839417697203, dt was forced below floating point epsilon 1.1102230246251565e-16, and step error estimate = 42.16290621054672. 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 ~/.julia/packages/SciMLBase/CMjVZ/src/integrator_interface.jl:623

The warning message we received can mean that there is a bug in the formulation of the problem. But if everything has been done correctly, it suggests that the solution may not exist past the indicated time. This is a possibility in nonlinear ODEs.

plot(sol, label="",
    xlabel=L"t",yaxis=(:log10,L"u(t)"), title="Finite-time blowup")

We can also produce an IVP that has more than one solution.

Example 6.1.5

The functions \(u(t)=u^2\) and \(u(t)\equiv 0\) both satisfy the differential equation \(u'=2\sqrt{u}\) and the initial condition \(u(0)=0\). Thus the corresponding IVP has more than one solution.

The following standard theorem gives us a condition that is easy to check and guarantees that a unique solution exists. But it is not the most general possible such condition, so there are problems with a unique solution that it cannot detect. We state the theorem without proof.

Theorem 6.1.6 : Existence and uniqueness

If the derivative \(\frac{\partial f}{\partial u}\) exists and \(\left|\frac{\partial f}{\partial u}\right|\) is bounded by a constant \(L\) for all \(a\le t \le b\) and all \(u\), then the initial-value problem (6.1.1) has a unique solution for \(t\in [a,b]\).

Conditioning of first-order IVPs#

In a numerical context we have to be concerned about the conditioning of the IVP. There are two key items in (6.1.1) that we might consider to be the data of the initial-value ODE problem: the function \(f(t,u)\), and the initial value \(u_0\). It’s easier to discuss perturbations to numbers than to functions, so we will focus on the effect of \(u_0\) on the solution, using the following theorem that we give without proof. Happily, its conditions are identical to those in Theorem 6.1.6.

Theorem 6.1.7 : Dependence on initial value

If the derivative \(\frac{\partial f}{\partial u}\) exists and \(\left|\frac{\partial f}{\partial u}\right|\) is bounded by a constant \(L\) for all \(a\le t \le b\) and all \(u\), then the solution \(u(t;u_0+\delta)\) of \(u'=f(t,u)\) with initial condition \(u(0)=u_0+\delta\) satisfies

(6.1.3)#\[\left\|u(t;u_0+\delta)-u(t;u_0)\right\|_\infty \le |\delta| e^{L(b-a)}\]

for all sufficiently small \(|\delta|\).

Numerical solutions of IVPs have errors, and those errors can be seen as perturbations to the solution. Theorem 6.1.7 gives an upper bound of \(e^{L(b-a)}\) on the infinity norm (i.e., pointwise) absolute condition number of the solution with respect to perturbations at an initial time. However, the upper bound may be a terrible overestimate of the actual sensitivity for a particular problem.

Demo 6.1.8

Consider the ODEs \(u'=u\) and \(u'=-u\). In each case we compute \(\partial f/\partial u = \pm 1\), so the condition number bound from Theorem 6.1.7 is \(e^{b-a}\) in both problems. However, they behave quite differently. In the case of exponential growth, \(u'=u\), the bound is the actual condition number.

But with \(u'=-u\), solutions actually get closer together with time.

In this case the actual condition number is one, because the initial difference between solutions is the largest over all time. Hence the exponentially growing bound \(e^{b-a}\) is a gross overestimate.

In general, solutions can diverge from, converge to, or oscillate around the original trajectory in response to perturbations. We won’t fully consider these behaviors and their implications for numerical methods again until a later chapter.

Exercises#

✍ For each IVP, determine whether the problem satisfies the conditions of Theorem 6.1.7). If so, determine the smallest possible value for \(L\).

(a) \(f(t,u) = 3 u,\; 0 \le t \le 1\)

(b) \(f(t,u) = -t \sin(u),\; 0 \le t \le 5\)

(c) \(f(t,u) = -(1+t^2) u^2,\; 1 \le t \le 3\)

(d) \(f(t,u) = \sqrt{u},\; 0 \le t \le 1\)
⌨ For each ODE in the preceding problem, assume that \(u\) is initially equal to \(1\) on the given interval. Solve the resulting IVP with solve and make a plot of the solution.
✍ Use an integrating factor to find the solution of each problem in analytic form.

(a) \(u' = -t u,\ 0 \le t \le 5,\ u(0) = 2\)

(b) \(u' - 3 u = e^{-2t},\ 0 \le t \le 1,\ u(0) = 5\)
✍ Consider the IVP \(u'=u^2\), \(u(0)=\alpha\).

(a) Does Theorem 6.1.6 apply to this problem?

(b) Show that \(u(t) = \alpha/(1-\alpha t)\) is a solution of the IVP.

(c) Does this solution necessarily exist for all \(t\in[0,1]\)?
⌨ Using solve, compute solutions \(x(t)\) to the logistic equation with harvesting,

\[x' = k (S-x)(x-M), \qquad 0\le t \le 10,\]

with \(k=S=1\) and \(M=0.25\), for the initial conditions \(x(0)=0.9M\), \(1.1M\), \(1.5M\), \(0.9S\), \(1.1S\), \(3S\). Show all the solutions together on one plot with \(0\le x \le 3\). (Note: One of the solutions will throw a warning and fail to reach \(t=10\), but you can plot it anyway.)
⌨ (a) Using solve, solve the IVP \(u'=u\cos(u) + \cos(4t)\), \(0\le t \le 10\), \(u(0) = u_0\) for \(u_0 = -2,-1.5,-1,\ldots,1.5,2\). Plot all the solutions on a single graph.

(b) All of the solutions in part (a) eventually settle into one of two periodic oscillations. To two digits of accuracy, find the value of \(u_0\) in \((-1,1)\) at which the selected long-term solution changes. (This will take repeated trials, narrowing down the range for \(u_0\) each time.)
⌨ Experimental evidence (see [NBL+81]) shows that a 300-mg oral dose of caffeine, such as might be found in a large mug of drip-brewed coffee, creates a concentration of about 8 \(\mu{\rm g}\)/mL in blood plasma. This boost is followed by first-order kinetics with a half-life of about 6 hours (although this rate can vary a great deal from person to person). We can model the caffeine concentration due to one drink taken over half an hour via

\[ x'(t) = -kx + C(t),\quad x(0)=0,\]

where \(k=\log(2)/6\) and

\[\begin{split} C(t) = \begin{cases} 16, & 0\le t \le 0.5, \\ 0, & t > 0.5. \end{cases}\end{split}\]

Use solve to make a plot of the caffeine concentration for 12 hours. Then change \(k=\log(2)/8\) (half-life of 8 hours) and plot the solution again.
⌨ A reasonable model of the velocity \(v(t)\) of a skydiver is

\[\frac{dv}{dt} = -g + \frac{k}{m}v^2, \qquad v(0)=0,\]

where \(g=9.8 \text{ m/sec}^2\) is gravitational acceleration, \(m\) is the mass of the skydiver with parachute, and \(k\) quantifies the effect of air resistance. At the US Air Force Academy, a training jump starts at about 1200 m and has \(k=0.4875\) for \(t<13\) and \(k=29.16\) or \(t\ge 13\). (This is an oversimplification; see [MS99].)

(a) Solve the IVP for \(v\) for an 80-kg cadet for \(t\in [0,200]\), and plot the solution.

(b) The total distance fallen up to time \(t\) is \(\displaystyle\int_0^t v(s)\, ds\). Use Function 5.7.2 to calculate and plot the altitude of the cadet as a function of time.

(c) In part (b), you should have found that the altitude becomes negative. Use Function 4.4.4 to determine accurately when the cadet reaches the ground.