93 Analytical solution of differential equations

1. Nonlinear differential equation

The only kind of nonlinear differential equations that we solve analytically is the so-called separable differential equations (including autonomous equations). See AMBS Ch 38–39. Solution method: separate the variables and integrate.

93.1. Find analytical solution formulas for the following initial value problems. In each case sketch the graphs of the solutions and determine the half-life. See: P. Atkins and L. Jones, Chemical Principles. The Quest for Insight. Freeman, New York, second edition, 2002, pp. 698–706.

(a) First order decay rate:

\begin{matrix}u′ = −ku, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(b) Second order decay rate:

\begin{matrix}u′ = −ku^{2}, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(c) Third order decay rate:

\begin{matrix}u′ = −ku^{3}, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

93.2. Find analytical solution formulas for the following initial value problems. Sketch the graphs. (a) First order increase rate:

\begin{matrix}u′ = ku, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(b) Second order increase rate:

\begin{matrix}u′ = ku^{2}, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(c) Third order increase rate:

\begin{matrix}u′ = ku^{3}, t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(d) The logistic equation:

\begin{matrix}u′ = ku(1 − u), t > 0, (k > 0, u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

(e)

\begin{matrix}u′ = −tu^{2}, t > 0, (u_{0} > 0) \\ u(0) = u_{0}. \end{matrix}

2. Linear differential equation

2.1 Linear differential equation—first order

(93.1)  u′ + a(t)u = f(t).

Here u = u(t) is an unknown function of an independent variable t. The equation is called homogeneous if f(t) ≡ 0 and nonhomogeneous otherwise. The differential operator Lu = u′+a(t)u has constant coefficient if a(t) = a is constant and it has variable coefficient otherwise. The equation is said to be a linear equation, because the operator L is a linear operator:

L(αu + βv) = αLu + βLv, (α, β ∈ R, u = u(t), v = v(t))

i.e., it preserves linear combinations of functions. Check this!

Solution method: multiply by the integrating factor e^{A(t)} with A(t) = ∫_{0}^{t} a(s) ds, and integrate. See AMBS Ch 35.1–2.

93.3. (constant coefficient, homogeneous) Solve the following. Sketch the graph of the solution.

(a)

\begin{matrix}u′ + 2u = 0, t > 0, \\ u(0) = u_{0}. \end{matrix}

(b)

\begin{matrix}u′ − 2u = 0, t > 0, \\ u(0) = u_{0}. \end{matrix}

93.4. (constant coefficient, nonhomogeneous) Solve the following.

(a)

\begin{matrix}u′ + 2u = f(t), t > 0, \\ u(0) = u_{0}. \end{matrix}

(b)

\begin{matrix}u′ − 2u = f(t), t > 0, \\ u(0) = u_{0}. \end{matrix}

93.5. (constant coefficient, nonhomogeneous) Solve the following.

\begin{matrix}u′ + au = f(t), t > 0, \\ u(0) = u_{0} \end{matrix}

93.6. (variable coefficient, nonhomogeneous) Solve the following.

\begin{matrix}u′ + 2tu = f(t), t > 0, \\ u(0) = u_{0}. \end{matrix}

2.2 Linear differential equation—second order—constant coefficients

(93.2)  u′′ + a_{1}u′ + a_{0}u = f(t).

The equation is called homogeneous if f(t) ≡ 0 and nonhomogeneous otherwise. We assume that the differential operator Lu = u′′ + a_{1}u′ + a_{0}u has constant coefficients a_{1} and a_{0}. Check that the operator L is linear!

Variable coefficients: Linear differential equations of second order with variable coefficients u′′ + a_{1}(t)u′ + a_{0}(t)u = f(t), cannot be solved analytically, except in some special cases. One such case can be found in AMBS Ch 35.6. We do not discuss this here.

Homogeneous equation

See AMBS Ch 35.3–35.4. The homogeneous equation (93.2) may be written

(93.3)  D^{2}u + a_{1}Du + a_{0}u = 0,

or

P (D)u = 0,

where

P (r) = r^{2} + a_{1}r + a_{0}

is the characteristic polynomial of the equation. The characteristic equation P (r) = 0 has two roots r_{1} and r_{2}. Hence P (r) = (r − r_{1})(r − r_{2}). All solutions of equation (93.2) are obtained as linear combinations

(93.4)  \begin{matrix}u(t) = c_{1}e^{r_{1}t} + c_{2}e^{r_{2}t}, & if r_{1} ≠ r_{2}, \\ u(t) = (c_{1} + c_{2}t)e^{r_{1}t}, & if r_{1} = r_{2}, \end{matrix}

where c_{2}, c_{2} are arbitrary coefficients. The coefficients may be determined from an initial condition of the form

u(0) = u_{0}, u′(0) = u_{1}.

The formula (93.4) is called the general solution of homogeneous linear equation (93.3). Example 93.1. We solve

u′′ + u′ − 12u = 0; u(0) = u_{0}, u′(0) = u_{1}.

The equation is written (D^{2} − D − 12)u = 0 and the characteristic equation is r^{2} + r − 12 = 0 with roots r_{1} = 3, r_{2} = −4. The general solution is

u(t) = c_{1}e^{3t} + c_{2}e^{−4t}

with the derivative

u′(t) = 3c_{1}e^{3t} − 4c_{2}e^{−4t}.

The initial condition gives

u_{0} = u(0) = c_{1} + c_{2}
u_{1} = u′(0) = 3c_{1} − 4c_{2}

which implies c_{1} = (4u_{0} + u_{1})/7, c_{2} = (3u_{0} − u_{1})/7. The solution is

u(t) = \frac{4u_{0} + u_{1}}{7} e^{3t} + \frac{3u_{0} − u_{1}}{7} e^{−4t}.

93.7. Prove the solution formula (93.4) by writing the equation as

P (D)u = (D − r_{1})(D − r_{2})u = 0

and by solving two first order equations (D − r_{1})v = 0 and (D − r_{2})u = v as in Problems 93.3 and 93.4.

93.8. Write the following equations as P (D)u = 0 and solve the initial value problem. Choose numerical values for the constants and sketch the graph of the solution. Solve the problem with your Matlab program my ode.m.

(a) u′′ − u′ − 2u = 0; u(0) = u_{0}, u′(0) = u_{1}.

(b) u′′ − k^{2}u = 0; u(0) = u_{0}, u′(0) = u_{1}.

(c) u′′ + 4u′ + 4u = 0; u(0) = u_{0}, u′(0) = u_{1}.

93.9. Solve the boundary value problem

\begin{matrix}u′′(x) − k^{2}u(x) = 0, 0 < x < L, \\ u(0) = 0, u(L) = u_{L}. \end{matrix}

Complex roots

If the characteristic polynomial P (r) has real coefficients, then its roots are either real numbers or a pair of conjugate complex numbers, see AMBS Ch 22.11. In the latter case we have r_{1} = α + iω and r_{2} = α − iω and the solution (93.4) becomes (see AMBS Ch 33.2 for the definition of exp(z) with a complex variable z)

u(t) = c_{1}e^{(α+iω)t} + c_{2}e^{(α−iω)t}
= e^{αt}(c_{1}e^{iωt} + c_{2}e^{−iωt})
= e^{αt}(c_{1}( cos(ωt) + i sin(ωt)) + c_{2}( cos(ωt) − i sin(ωt)))
= e^{αt}(d_{1} cos(ωt) + d_{2} sin(ωt)),

with d_{1} = c_{1} + c_{2}, d_{2} = i(c_{1} − c_{2}).

93.10. Write the equation as P (D)u = 0 and solve the initial value problem. Choose numerical values for the constants and sketch the graph of the solution. Solve the problem with your Matlab program my ode.m.

(a) u′′ + 4u′ + 13u = 0; u(0) = u_{0}, u′(0) = u_{1}.

(b) u′′ + ω^{2}u = 0; u(0) = u_{0}, u′(0) = u_{1}.

Nonhomogeneous equation

See AMBS Ch 35.5. The solution of the nonhomogeneous equation P (D)u = f(t) is given by

(93.5)  u(t) = u_{h}(t) + u_{p}(t),

where u_{h} is the general solution (93.4) of the corresponding homogeneous equation, i.e., P (D)u_{h} = 0, and u_{p} is a particular solution of the nonhomogeneous equation, i.e., P (D)u_{p} = f(t).

Proof: If u is given by (93.5), then Lu = L(u_{h} + u_{p}) = Lu_{h} + Lu_{p} = 0 + f = f, so that u solves the nonhomogeneous equation. On the other hand: if u_{p} is a particular solution and u is any other solution of the nonhomogeneous equation, then L(u − u_{p}) = Lu − Lu_{p} = f − f = 0, i.e., u − u_{p} solves the homogeneous equation. Thus u − u_{p} = u_{h}, which is (93.5).

A particular solution can sometimes be found by guess-work: we make an Ansatz for u_{p} of the same form as f.

Example 93.2. u′′ − u′ − 2u = t. Here f(t) = t is a polynomial of degree 1 and we make the Ansatz u_{p}(t) = At + B, i.e., a polynomial of degree 1. Substitution into the equation gives −A − 2(At + B) = t. Identification of coefficients gives A = −\frac{1}{2}, B = \frac{1}{4}, so that u_{p}(t) = \frac{1}{4} − \frac{1}{2}t. The general solution of the homogeneous equation is u_{h}(t) = c_{1}e^{−t} + c_{2}e^{2t}, see Problem 93.8 (a). Hence we get

u(t) = u_{h}(t) + u_{p}(t) = c_{1}e^{−t} + c_{2}e^{2t} + \frac{1}{4} − \frac{1}{2}t.

93.11. Solve the following.

(a) u′′ − u′ − 2u = e^{t} Ansatz: u_{p}(t) = Ae^{t}

(b) u′′ − u′ − 2u = cos(t) Ansatz: u_{p}(t) = A cos(t) + B sin(t)

(c) u′′ − u′ − 2u = t^{3} Ansatz: u_{p}(t) = At^{3} + Bt^{2} + Ct + D

(d) u′′ − u′ − 2u = e^{−t} Ansatz: u_{p}(t) = Ate^{−t} Re-writing as a system of first order equations By setting w_{1} = u, w_{2} = u′, w = [\begin{matrix}w_{1} \\ w_{2} \end{matrix}], we can re-write (93.2) as a system of first order equations

w′(t) = Aw(t) + F (t); w(0) = w_{0},

where

w_{0} = [\begin{matrix}u_{0} \\ u_{1} \end{matrix}] , A = [\begin{matrix}0 & 1 \\ −a_{0} & −a_{1} \end{matrix}] , F (t) = [\begin{matrix}0 \\ f(t) \end{matrix}] .

To see this we compute

w′ = [\begin{matrix}u′ \\ u′′ \end{matrix}] = [\begin{matrix}u′ \\ −a_{0}u − a_{1}u′ + f(t) \end{matrix}] = [\begin{matrix}w_{2} \\ −a_{0}w_{1} − a_{1}w_{2} + f(t) \end{matrix}] = [\begin{matrix}0 & 1 \\ −a_{0} & −a_{1} \end{matrix}] [\begin{matrix}w_{1} \\ w_{2} \end{matrix}] + [\begin{matrix}0 \\ f(t) \end{matrix}] .

It is necessary to do this re-writing before we can use our Matlab programs to solve (93.2). 93.12. Write the equation in Problem 93.11(a) as a system of first order.

2.3 System of linear differential equations of first order

Constant coefficients—homogeneous equations We finally mention

(93.6)  \begin{matrix}u′ + Au = 0, t > 0, \\ u(0) = u_{0}, \end{matrix}

where u(t), u_{0} ∈ R^{d}, and A ∈ R^{d×d} is a constant matrix of coefficients. This kind of system will studied by means of eigenvalues and eigenvectors in the following course ALA-C.

Answers and solutions

93.1. Reaction of order 1 (decay rate of order 1):

\{ \begin{matrix}u′ = −ku \\ u(0) = u_{0} \end{matrix}
u(t) = u_{0}e^{−kt}

The half-life T_{1/2} is given by

u(T_{1/2}) = u_{0}e^{−kT_{1/2}} = \frac{1}{2}u_{0},

which leads to

T_{1/2} = \frac{log(2)}{k}.

Reaction of order n > 1 (decay rate of order n > 1):

\{ \begin{matrix}u′ = −ku^{n} \\ u(0) = u_{0} \end{matrix}

\frac{du}{u^{n}} = −k dt

∫_{u_{0}}^{u(T)} u^{−n} du = − ∫_{0}^{T} k dt

[\frac{u^{−n+1}}{−n + 1}]_{u_{0}}^{u(T)} = −kT

u(T)^{−n+1} − u_{0}^{−n+1} = (n − 1)kT

\frac{1}{u(T)^{n−1}} = \frac{1}{u_{0}^{n−1}} + (n − 1)kT = \frac{1 + (n − 1)u_{0}^{n−1}kT}{u_{0}^{n−1}}

u(T) = \frac{u_{0}}{(1 + (n − 1)u_{0}^{n−1}kT )^{1/(n−1)}}

The half-life T_{1/2} is given by

u(T_{1/2}) = \frac{u_{0}}{(1 + (n − 1)u_{0}^{n−1}kT_{1/2})^{1/(n−1)}} = \frac{1}{2}u_{0}

which implies

T_{1/2} = \frac{2^{n−1} − 1}{(n − 1)u_{0}^{n−1}k}

93.2.

(a) u(t) = u_{0}e^{kt}

(b–c) order n > 1: u(t) = \frac{u_{0}}{(1 − (n − 1)u_{0}^{n−1}kt)^{1/(n−1)}} , 0 ≤ t < \frac{1}{(n − 1)u_{0}^{n−1}k}

(d) u(t) = \frac{u_{0}}{u_{0} + (1 − u_{0})e^{−kt}}
(e) u(t) = \frac{u_{0}}{1 + t^{2}u_{0}/2}

93.3.

(a) u(t) = e^{−2t}u_{0}

(b) u(t) = e^{2t}u_{0} 93.4.

(a) u(t) = e^{−2t}u_{0} + ∫_{0}^{t} e^{−2(t−s)}f(s) ds

(b) u(t) = e^{2t}u_{0} + ∫_{0}^{t} e^{2(t−s)}f(s) ds 93.5. u(t) = e^{−at}u_{0} + ∫_{0}^{t} e^{−a(t−s)}f(s) ds 93.6. Integrating factor: e^{t^{2}}. Solution u(t) = e^{−t^{2}}u_{0} + ∫_{0}^{t} e^{−(t^{2}−s^{2})}f(s) ds.

93.7. The equation for v is (D − r_{1})v = 0, or v′ − r_{1}v = 0, v(0) = v_{0}, with unique solution v(t) = v_{0}e^{r_{1}t}. The equation for u is (D − r_{2})u = v, or u′ − r_{2}u = v_{0}e^{r_{1}t}, u(0) = u_{0}, with unique solution, see Problem 93.5,

u(t) = u_{0}e^{r_{2}t} + e^{r_{2}t} ∫_{0}^{t} e^{−r_{2}s}v_{0}e^{r_{1}s} ds

= u_{0}e^{r_{2}t} + v_{0}e^{r_{2}t} ∫_{0}^{t} e^{(r_{1}−r_{2})s} ds

= u_{0}e^{r_{2}t} + v_{0}e^{r_{2}t}[\frac{e^{(r_{1}−r_{2})s}}{r_{1} − r_{2}} ]_{s=0}^{t}

= u_{0}e^{r_{2}t} + \frac{v_{0}}{r_{1} − r_{2}} (e^{r_{1}t} − e^{r_{2}t})

= \frac{v_{0}}{r_{1} − r_{2}} e^{r_{1}t} + (u_{0} − \frac{v_{0}}{r_{1} − r_{2}} )e^{r_{2}t}

= c_{1}e^{r_{1}t} + c_{2}e^{r_{2}t}, if r_{1} ≠ r_{2}.

If r_{1} = r_{2}, then we have instead

u(t) = u_{0}e^{r_{2}t} + e^{r_{2}t} ∫_{0}^{t} e^{−r_{2}s}v_{0}e^{r_{1}s} ds = u_{0}e^{r_{1}t} + v_{0}e^{r_{1}t} ∫_{0}^{t} ds = u_{0}e^{r_{1}t} + v_{0}te^{r_{1}t} = (c_{1} + c_{2}t)e^{r_{1}t}.

93.8.

(a) u(t) = \frac{1}{3}(2u_{0} − u_{1})e^{−t} + \frac{1}{3}(u_{0} + u_{1})e^{2t}.

(b) u(t) = c_{1}e^{kt} + c_{2}e^{−kt} = d_{1} cosh(kt) + d_{2} sinh(kt), d_{1} = c_{1} + c_{2}, d_{2} = c_{1} − c_{2}. The initial condition gives u(t) = \frac{1}{2}(u_{0}+u_{1}/k)e^{kt}+ \frac{1}{2}(u_{0}−u_{1}/k)e^{−kt} or alternatively u(t) = u_{0} cosh(kt)+ (u_{1}/k) sinh(kt).

(c) u(t) = (u_{0} + (2u_{0} + u_{1})t)e^{−2t}. 93.9. u(x) = u_{L} sinh(kx)/ sinh(kL). 93.10.

(a) u(t) = e^{−2t}(u_{0} cos(3t) + \frac{1}{3}(2u_{0} + u_{1}) sin(3t)).

(b) u(t) = u_{0} cos(ωt) + (u_{1}/ω) sin(ωt). Compare to Problem 93.8 (b). 93.11.

(a) u(t) = c_{1}e^{−t} + c_{2}e^{2t} − \frac{1}{2}e^{t}.

(b) u(t) = c_{1}e^{−t} + c_{2}e^{2t} − \frac{3}{10} cos(t) − \frac{1}{10} sin(t).

(c) u(t) = c_{1}e^{−t} + c_{2}e^{2t} − \frac{1}{2}t^{3} + \frac{3}{4}t^{2} − \frac{9}{4}t + \frac{15}{8}.

(d) u(t) = c_{1}e^{−t} + c_{2}e^{2t} − \frac{1}{3}te^{−t}. Note: the Ansatz u_{p}(t) = Ae^{−t} does not work, because e^{−t} is a solution of the homogeneous equation, P (D)e^{−t} = 0, i.e., e^{−t} is contained in u_{h}.

93.12. [\begin{matrix}w′_{1} \\ w′_{2} \end{matrix}] = [\begin{matrix}0 & 1 \\ 2 & 1 \end{matrix}] [\begin{matrix}w_{1} \\ w_{2} \end{matrix}] + [\begin{matrix}0 \\ e^{t} \end{matrix}].

2004-11-21 /stig