Basics of IVPs
clear all
format short
set(0, 'defaultaxesfontsize', 12)
set(0, 'defaultlinelinewidth', 1.5)
set(0, 'defaultFunctionLinelinewidth', 1.5)
set(0, 'defaultscattermarkerfacecolor', 'flat')
gcf;
set(gcf, 'Position', [0 0 600 350], 'Theme', 'light')
addpath ../FNC_matlab6.1. Basics of IVP s ¶ A scalar first-order initial-value problem IVP ) is
u ′ ( t ) = f ( t , u ( t ) ) , a ≤ t ≤ b , u ( a ) = u 0 . \begin{split}
u'(t) &= f(t,u(t)), \qquad a \le t \le b, \\
u(a) &=u_0.
\end{split} u ′ ( t ) u ( a ) = f ( t , u ( t )) , a ≤ t ≤ b , = u 0 . We call t t t independent variable and u u u dependent variable . If u ′ = f ( t , u ) = g ( t ) + u h ( t ) u'=f(t,u)=g(t)+u h(t) u ′ = f ( t , u ) = g ( t ) + u h ( t ) linear ; otherwise, it is nonlinear .
A solution of an initial-value problem is a function u ( t ) u(t) u ( t ) u ′ ( t ) = f ( t , u ( t ) ) u'(t)=f\bigl(t,u(t)\bigr) u ′ ( t ) = f ( t , u ( t ) ) u ( a ) = u 0 u(a)=u_0 u ( a ) = u 0
When t t t u ˙ \dot{u} u ˙ u ′ u' u ′
Suppose u ( t ) u(t) u ( t ) t t t u u u
d u d t = k u , u ( 0 ) = u 0 \frac{d u}{d t} = k u, \qquad u(0)=u_0 d t d u = k u , u ( 0 ) = u 0 for some k > 0 k>0 k > 0 u ( t ) = e k t u 0 u(t)=e^{kt}u_0 u ( t ) = e k t u 0
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
d u d t = k u − r u 2 , u ( 0 ) = u 0 . \frac{d u}{d t} = ku - ru^2, \qquad u(0)=u_0. d t d u = k u − r u 2 , 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 ) = k / r 1 + ( k r u 0 − 1 ) e − k t . u(t) = \frac{k/r}{ 1 + \left( \frac{k}{r u_0} - 1 \right) e^{-k t} }. u ( t ) = 1 + ( r u 0 k − 1 ) e − k t k / r . For k , r , u 0 > 0 k,r,u_0>0 k , r , u 0 > 0 u 0 u_0 u 0 k / r k/r k / r
Linear problems can be solved in terms of integrals. Defining the integrating factor ρ ( t ) = exp [ ∫ − h ( t ) d t ] \rho(t) = \exp\bigl[\int -h(t)\, dt \bigr] ρ ( t ) = exp [ ∫ − h ( t ) d t ]
ρ ( t ) u ( t ) = u 0 + ∫ a t ρ ( s ) g ( s ) d s . \rho(t) u(t) = u_0 + \int_a^t \rho(s) g(s) \, ds. ρ ( t ) u ( t ) = u 0 + ∫ a t ρ ( s ) g ( s ) d s . In many cases, however, the necessary integrals cannot be done in closed form. Some nonlinear ODE s, 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 ( t , u , u ′ ) u''(t)=f\bigl(t,u,u'\bigr) u ′′ ( t ) = f ( t , u , u ′ ) IVP requires two conditions at the initial time in order to specify a solution completely. As we will see in IVP systemsIVP s in a first-order form, so we will deal with first-order problems exclusively.
6.1.1 Numerical solutions ¶ Let’s use a built-in method to define and solve the initial-value problem
u ′ = sin [ ( u + t ) 2 ] , t ∈ [ 0 , 4 ] , u ( 0 ) = − 1. u'=\sin[(u+t)^2], \quad t \in [0,4], \quad u(0)=-1. u ′ = sin [( u + t ) 2 ] , t ∈ [ 0 , 4 ] , u ( 0 ) = − 1. We must create a function that computes u ′ u' u ′ u u u
f = @(t, u) sin((t + u)^2);
u0 = -1;
tspan = [0, 4];
These ingredients are supplied to a solver function. A standard first choice of solver is called ode45.
[t, u] = ode45(f, [0, 4], u0);
fprintf("length(t) = %d, length(u) = %d\n", length(t), length(u))The solution is represented as a pair of vectors, t and u, where t contains the times at which the solution was computed and u contains the corresponding values of u ( t ) u(t) u ( t )
clf
plot(t, u, '-o')
xlabel("t")
ylabel("u(t)")
title(("Solution of an IVP"));You might want to know the solution at particular times other than the ones selected by the solver. That requires an interpolation. To do this automatically requires calling the solver with just one output argument, which is then supplied to deval to evaluate the solution at any time.
sol = ode45(f, [0, 4], u0);
u = @(t) deval(sol, t)'; % return a column vector
u(0:0.5:2)6.1.2 Existence and uniqueness ¶ There are simple IVP s that do not have solutions at all possible times.
The equation u ′ = ( u + t ) 2 u'=(u+t)^2 u ′ = ( u + t ) 2
f = @(t, u) (t + u)^2;
u0 = 1;
tspan = [0, 1];
[t, u] = ode45(f, tspan, u0);Warning: Failure at t=7.853789e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t. 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 simply may not exist past the indicated time. This is a possibility in nonlinear ODE s.
clf
semilogy(t, u)
xlabel("t")
ylabel("u(t)")
title(("Finite-time blowup"));We can also produce an IVP that 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.
If the derivative ∂ f ∂ u \frac{\partial f}{\partial u} ∂ u ∂ f ∣ ∂ f ∂ u ∣ \left|\frac{\partial f}{\partial u}\right| ∣ ∣ ∂ u ∂ f ∣ ∣ L L L a ≤ t ≤ b a\le t \le b a ≤ t ≤ b u u u (6.1.1) t ∈ [ a , b ] t\in [a,b] t ∈ [ a , b ]
6.1.3 Conditioning of first-order IVP s ¶ In a numerical context we have to be concerned about the conditioning of the IVP . There are two key items in (6.1.1) ODE problem: the function f ( t , u ) f(t,u) f ( t , u ) u 0 u_0 u 0 u 0 u_0 u 0 Theorem 6.1.1
If the derivative ∂ f ∂ u \frac{\partial f}{\partial u} ∂ u ∂ f ∣ ∂ f ∂ u ∣ \left|\frac{\partial f}{\partial u}\right| ∣ ∣ ∂ u ∂ f ∣ ∣ L L L a ≤ t ≤ b a\le t \le b a ≤ t ≤ b u u u u ( t ; u 0 + δ ) u(t;u_0+\delta) u ( t ; u 0 + δ ) u ′ = f ( t , u ) u'=f(t,u) u ′ = f ( t , u ) u ( a ) = u 0 + δ u(a)=u_0+\delta u ( a ) = u 0 + δ
∥ u ( t ; u 0 + δ ) − u ( t ; u 0 ) ∥ ∞ ≤ ∣ δ ∣ e L ( b − a ) \left\|u(t;u_0+\delta)-u(t;u_0)\right\|_\infty \le |\delta| e^{L(b-a)} ∥ u ( t ; u 0 + δ ) − u ( t ; u 0 ) ∥ ∞ ≤ ∣ δ ∣ e L ( b − a ) for all sufficiently small ∣ δ ∣ |\delta| ∣ δ ∣
Numerical solutions of IVP s have errors, and those errors can be seen as perturbations to the solution. Theorem 6.1.2 e L ( b − a ) e^{L(b-a)} e L ( b − a )
Consider the ODE s u ′ = u u'=u u ′ = u u ′ = − u u'=-u u ′ = − u ∂ f / ∂ u = ± 1 \partial f/\partial u = \pm 1 ∂ f / ∂ u = ± 1 Theorem 6.1.2 e b − a e^{b-a} e b − a u ′ = u u'=u u ′ = u
clf
for u0 = [0.7, 1, 1.3] % initial values
fplot(@(t) exp(t) * u0, [0, 3]), hold on
end
xlabel('t')
ylabel('u(t)')
title(('Exponential divergence of solutions'));But with u ′ = − u u'=-u u ′ = − u
clf
for u0 = [0.7, 1, 1.3] % initial values
fplot(@(t) exp(-t) * u0, [0, 3]), hold on
end
xlabel('t')
ylabel('u(t)')
title(('Exponential convergence of solutions'));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 e^{b-a} e b − a
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.
6.1.4 Exercises ¶ ✍ For each IVP , determine whether the problem satisfies the conditions of Theorem 6.1.2 L L L
(a) f ( t , u ) = 3 u , 0 ≤ t ≤ 1 f(t,u) = 3 u,\; 0 \le t \le 1 f ( t , u ) = 3 u , 0 ≤ t ≤ 1
(b) f ( t , u ) = − t sin ( u ) , 0 ≤ t ≤ 5 f(t,u) = -t \sin(u),\; 0 \le t \le 5 f ( t , u ) = − t sin ( u ) , 0 ≤ t ≤ 5
(c) f ( t , u ) = − ( 1 + t 2 ) u 2 , 1 ≤ t ≤ 3 f(t,u) = -(1+t^2) u^2,\; 1 \le t \le 3 f ( t , u ) = − ( 1 + t 2 ) u 2 , 1 ≤ t ≤ 3
(d) f ( t , u ) = u , 0 ≤ t ≤ 1 f(t,u) = \sqrt{u},\; 0 \le t \le 1 f ( t , u ) = u , 0 ≤ t ≤ 1
✍ Use an integrating factor to find the solution of each problem in analytic form.
(a) u ′ = − t u , 0 ≤ t ≤ 5 , u ( 0 ) = 2 u' = -t u,\ 0 \le t \le 5,\ u(0) = 2 u ′ = − t u , 0 ≤ t ≤ 5 , u ( 0 ) = 2
(b) u ′ − 3 u = e − 2 t , 0 ≤ t ≤ 1 , u ( 0 ) = 5 u' - 3 u = e^{-2t},\ 0 \le t \le 1,\ u(0) = 5 u ′ − 3 u = e − 2 t , 0 ≤ t ≤ 1 , u ( 0 ) = 5
✍ Consider the IVP u ′ = u 2 u'=u^2 u ′ = u 2 u ( 0 ) = α u(0)=\alpha u ( 0 ) = α
(a) Does Theorem 6.1.1
(b) Show that u ( t ) = α 1 − α t u(t) = \dfrac{\alpha}{1-\alpha t} u ( t ) = 1 − α t α IVP .
(c) Does this solution necessarily exist for all t ∈ [ 0 , 1 ] t\in[0,1] t ∈ [ 0 , 1 ]
⌨ Using the method in Example 6.1.2 x ( t ) x(t) x ( t ) logistic equation with harvesting,
x ′ = k ( S − x ) ( x − M ) , 0 ≤ t ≤ 10 , x' = k (S-x)(x-M), \qquad 0\le t \le 10, x ′ = k ( S − x ) ( x − M ) , 0 ≤ t ≤ 10 , with k = S = 1 k=S=1 k = S = 1 M = 0.25 M=0.25 M = 0.25 x ( 0 ) = 0.9 M x(0)=0.9M x ( 0 ) = 0.9 M 1.1 M 1.1M 1.1 M 1.5 M 1.5M 1.5 M 0.9 S 0.9S 0.9 S 1.1 S 1.1S 1.1 S 3 S 3S 3 S 0 ≤ x ≤ 3 0\le x \le 3 0 ≤ x ≤ 3 t = 10 t=10 t = 10
⌨ (a) Using the method in Example 6.1.2 IVP u ′ = u cos ( u ) + cos ( 4 t ) u'=u\cos(u) + \cos(4t) u ′ = u cos ( u ) + cos ( 4 t ) 0 ≤ t ≤ 10 0\le t \le 10 0 ≤ t ≤ 10 u ( 0 ) = u 0 u(0) = u_0 u ( 0 ) = u 0 u 0 = − 2 , − 1.5 , − 1 , … , 1.5 , 2 u_0 = -2,-1.5,-1,\ldots,1.5,2 u 0 = − 2 , − 1.5 , − 1 , … , 1.5 , 2
(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 u_0 u 0 ( − 1 , 1 ) (-1,1) ( − 1 , 1 ) u 0 u_0 u 0
⌨ Experimental evidence (see Newton et al. (1981) μ g \mu{\rm g} μ g
x ′ ( t ) = − k x + C ( t ) , x ( 0 ) = 0 , x'(t) = -kx + C(t),\quad x(0)=0, x ′ ( t ) = − k x + C ( t ) , x ( 0 ) = 0 , where k = log ( 2 ) / 6 k=\log(2)/6 k = log ( 2 ) /6
C ( t ) = { 16 , 0 ≤ t ≤ 0.5 , 0 , t > 0.5. C(t) =
\begin{cases}
16, & 0\le t \le 0.5, \\
0, & t > 0.5.
\end{cases} C ( t ) = { 16 , 0 , 0 ≤ t ≤ 0.5 , t > 0.5. Solve the IVP and make a plot of the caffeine concentration for 12 hours. Then change k = log ( 2 ) / 8 k=\log(2)/8 k = log ( 2 ) /8
⌨ A reasonable model of the velocity v ( t ) v(t) v ( t )
d v d t = − g + k ( t ) m v 2 , v ( 0 ) = 0 , \frac{dv}{dt} = -g + \frac{k(t)}{m}v^2, \quad v(0)=0, d t d v = − g + m k ( t ) v 2 , v ( 0 ) = 0 , where g = 9.8 m/sec 2 g=9.8 \text{ m/sec}^2 g = 9.8 m/sec 2 m m m k k k k = 0.4875 k=0.4875 k = 0.4875 t < 13 t<13 t < 13 k = 29.16 k=29.16 k = 29.16 t ≥ 13 t\ge 13 t ≥ 13 Meade & Struthers (1999)
(a) Solve the IVP for v v v t ∈ [ 0 , 200 ] t\in [0,200] t ∈ [ 0 , 200 ]
(b) The total distance fallen up to time t t t ∫ 0 t v ( s ) d s \displaystyle\int_0^t v(s)\, ds ∫ 0 t v ( s ) d s Function 5.7.1
(c) In part (b), you should have found that the altitude becomes negative. Use Function 4.4.1
Newton, R., Broughton, L. J., Lind, M. J., Morrison, P. J., Rogers, H. J., & Bradbrook, I. D. (1981). Plasma and Salivary Pharmacokinetics of Caffeine in Man. European Journal of Clinical Pharmacology , 21 (1), 45–52. 10.1007/BF00609587 Meade, D. B., & Struthers, A. A. (1999). Differential Equations in the New Millennium: The Parachute Problem. International Journal of Engineering Education , 15 (6), 417–424.