First-order ODEs

ODEs of first order contain only the first derivative \(y'\), so they take the form

(1)\[ F(x, y, y') = 0 \]

or often

(2)\[ y' = \frac{dy}{dx} = f(x, y) \]

where \(y\) is the dependent variable and \(x\) is the independent variable and \(f\) is any function of them. The ODE in eq. (1) is written in implicit form, while it is written explicitly in eq. (2).

Direction field

import numpy as np
import matplotlib.pyplpot as plt

Direct integration

The basic technique in solving differential equations is integration. We can directly integrate an ODE of the form

\[ \frac{du}{dx} = f(x) \]

Separation of variables

If we can reduce an ODE through algebraic manipulations to the form

\[ \frac{du}{dx} = g(x)h(u) \]

we call it a separable equation. If \(h(u) = 0\) then the solution is \(u =\) const. If \(h(y) \neq 0\) we can divide through by it and integrate with respect to \(x\):

\[ \int \frac{du}{h(u)} = \int g(x) dx + c \]

to find the general solution.

Example 2. Let us find the general solution of

\[ u' = e^{2x-1} u^2 \]

This is a separable equation … integrate both sides w.r.t. \(x\):

\[ \int \frac{du}{u^2} = \int e^{2x-1} dx \]

Example: \(y' + ky = 0, \quad k \in \mathbb{R}\)

Consider a simple, but very important, ODE \(y' = -ky\), where \(k\) is a constant. The question we need to be asking ourselves is “what function when differentiated is equal to itself times some \(-k\)?” The answer is the exponential function

\[ y = Ce^{-kx} \quad \longrightarrow \quad y' = -kC e^{-kx} \]

With practice this type of questions and the answers to them will come more naturally. If an answer is not obvious (and most often it will not be), you could try solving it using some method. Let us do that:

\[\begin{split} \frac{dy}{dx} = -k y \\ \frac{dy}{y} dx = -k \\ \int \frac{dy}{y} = \int -k dx \\ \ln y = -kx + \tilde{C} \\ y = e^{-kx + \tilde{C}} = Ce^{-kx}\end{split}\]

Reduction to separable form

Exact ODEs

If a function \(u(x, y)\) has continuous partial derivatives, its differential (also called its total differential) is

\[ du = \frac{\partial u}{\partial x} dx + \frac{\partial u}{\partial y} dy \]

From here it follows that if \(u(x, y) =\) const, then \(du = 0\).

A first-order ODE \(M(x, y) + N(x, y)y' = 0\), written as

(3)\[ M(x, y) dx + N(x, y) dy = 0 \]

is called an exact differential equation if the differential form \(M(x, y) dx + N(x, y) dy\) is exact, that is, this form is the differential

(4)\[ du = \frac{\partial u}{\partial x} dx + \frac{\partial u}{\partial y} dy \]

of some function \(u(x, y)\). Then (1) can be written

\[ du = 0 \]

By integration we immediately obtain the general solution of (1) in the form

\[ u(x, y) = c \]

This is called an implicit solution

Comparing (1) and (2), we see that (1) is an exact differential equation if there is some function \(u(x, y)\) such that

(5)\[ \frac{\partial u}{\partial x} = M, \qquad \frac{\partial u}{\partial y} = N \]

Necessary and sufficient condition for (1) to be an exact differential equation

\[ \frac{\partial M}{\partial y} = \frac{\partial^2 u}{\partial y \partial x} = \frac{\partial^2 u}{\partial x \partial y} = \frac{\partial N}{\partial x} \]

If (1) is exact, the function \(u(x, y)\) can be found by inspection or in the following systematic way. From 4a we have by integration w.r.t. x

(6)\[ u = \int M dx + k(y) \]

where \(k(y)\) is the “constant” of integration. To determine \(k(y)\), we derive \(\partial u / \partial y\) from (6), use (4b) to get \(dk/dy\), and integrate \(dk/dy\) to get \(k\).

Could have done the same with 4b to get

\[ u = \int N dy + l(x) \]

Physical interpretation: conservative force

Recall that a force \(\mathbf{F}\) is conservative if the work done while moving between two points is independent of the path taken. Then it can be expressed as a gradient of potential \(\mathbf{F} = \nabla \Phi = (\Phi_x, \Phi_y)\).

Consider a force field \(\mathbf{F} = (P, Q)\) along an infinitesimal path \(\mathbf{dr} = (dx, dy)\)

\[ \mathbf{F \ dr} = Pdx + Qdy = \frac{\partial \Phi}{\partial x} dx + \frac{\partial \Phi}{\partial y} dy \]

The exactness condition is

\[ \frac{\partial P}{\partial y} = \frac{\partial^2 \Phi}{\partial x \partial y} = \frac{\partial^2 \Phi}{\partial y \partial x} = \frac{\partial Q}{\partial x} \]

which is satisfied. Therefore, the exactness of a differential equation can be identified with physical systems in which there is a potential.

Reduction to exact form. Integrating factors

Nonexact equation

(7)\[ P(x, y) dx + Q(x, y) dy = 0 \]

multiply by a function F that will in general be a function of both x and y

\[ FP dx + FQ dy = 0 \]

which is exact. F(x, y) is then called an integrating factor of that equation.

HOW TO FIND INT FACTORS

In simpler cases it may be found by inspection. In general:

The exactness condition \(\partial M / \partial y = \partial N / \partial x\), so for \(FP dx + FQ dy = 0\) the exactness condition is

\[ \frac{\partial}{\partial y} (FP) = \frac{\partial}{\partial x} (FQ) \]

\[ F_y P + FP_y = F_x Q + FQ_x \]

In general this would be difficult to solve, so we try a simpler case. In many practical cases there are factors that depend only on one variable. We let \(F = F(x)\) so the above becomes

\[ FP_y = F'Q + FQ_x \]

(8)\[ \frac{1}{F} \frac{d F}{dx} = R, \quad where \quad R = \frac{1}{Q} \left( \frac{\partial P}{\partial y} - \frac{\partial Q}{\partial x} \right) \]

This proves the following theorem: If (12) is such that the right side R of 16 depends only on x, then (12) has an integrating factor F = F(x) which is obtained by integrating 16 and taking exponents on both sides

\[ F(x) = \exp \int R(x) dx \]

Similarly, if F* = F*(y), then instead of 16

\[ \frac{1}{F*} \frac{d F*}{dy} = R*, \quad where \quad R* = \frac{1}{P} \left( \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \right) \]

and

\[ F*(y) = \exp \int R*(y) dy \]

First order linear differential equation

A general first-order linear ODE is of the form

\[ u' + p(x) u = r(x) \]

In engineering r(x) is often called the input and y(x) the output of the response to the input (and, if given, to the initial condition)

We find the solution to the homogeneous equation \(u_h(x)\) and then one particular solution \(u_p(x)\) to obtain the complete solution \(u\). The corresponding homogeneous equation

\[ u'(x) + p(x) u(x) = 0 \]

is separable, so the solution can be found by separation of parameters

\[ \frac{u'}{u} = -p \]

and integrating from \(x_0\) to \(x\)

\[\begin{split} \ln u(x) - \ln u(x_0) = - \int_{x_0}^x p(s) ds \\ u_h(x) = C e^{- \int_{x_0}^x p(s) ds}\end{split}\]

where \(s\) is a dummy variable that we use so to not confuse the integration variable with the upper integration limit \(x\). The reason why we are integrating from \(x_0\) to \(x\) will become clear very soon.

We find the particular solution by variation of parameters where \(C\) is no longer a constant but a function of \(x\):

\[ u_p(x) = C(x) e^{- \int_{x_0}^x p(s) ds} \]

We plug this back into the original equation and solve for \(C(x)\)

\[\begin{split} u'_p + pu_p = (C' - Cp) e^{- \int p} + pCe^{- \int p} = C' e^{- \int p} = r \\ C' = r e^{\int p} \\ C(x) = \int_{x_0}^x r(s) e^{\int_{x_0}^s p(\sigma) d\sigma} ds \end{split}\]

The solution to the inhomogeneous problem is

\[\begin{split} \begin{aligned} u(x) & = u_p + u_h \\ & = e^{- \int_{x_0}^x p(s) \ ds} \left ( \int_{x_0}^x r(s) \ e^{\int_{x_0}^s p(\sigma) d\sigma} \ ds + C \right ) \\ & = e^{- \int_{x_0}^x p(s) \ ds} \left ( \int_{x_0}^x r(s) \ e^{\int_{x_0}^s p(\sigma) d\sigma} \ ds + u(x_0) \right ) \end{aligned} \end{split}\]

No integration limit

Kreyszig omits the integration limits when writing this. Doing this only changes the constant of integration:

\[ u(x) = e^{-h} \left ( \int e^h r \ dx + c \right), \quad h = \int p(x) \ dx \]

where \(C = u(x_0)\) because for \(x=x_0\) all integrals are \(0\); this is why we integrated from \(x_0\) but, of course, we did not have to. If we did not, only the constant would be different. Rewriting this in a more convenient way

\[ u(x) = e^{-h} \left ( \int_{x_0}^x e^h r \ ds + u(x_0) \right), \quad h = \int_{x_0}^x \ p(s) \ ds \]

Note. We can rewrite the above as

\[ u(x) e^{\int_{x_0}^x \ p(s) \ ds} = \int_{x_0}^x \ r(s) \ e^{\int_{x_0}^s p(\sigma) d\sigma} \ ds + u(x_0) \]

After differentiating both sides we get

\[ u' e^{\int_{x_0}^x \ p(s) \ ds} + p(x) u e^{\int_{x_0}^x \ p(s) \ ds} = r(x) \ e^{\int_{x_0}^x \ p(s) \ ds}\]

We see that \(e^{\int_{x_0}^x p(s) ds}\) is the integrating factor so we could have found the solution by multiplying both sides of the equation by it.

Example: \(xy' + 4y = 8x^4, \quad y(1) = 2\)

(Kreyszig PS 1.5, Q7)

We divide both sides by \(x\) to get standard form:

\[ y' + \frac{4}{x}y = 8x^3, \quad y(1) = 2 \]

Now we see that \(p(x) = \frac{4}{x}, r(x) = 8x^3\) and the condition \(y(x_0) = y(1) = 2\). We first calculate \(h\)

\[ h = \int_{x_0}^x p(s) ds = \int_1^x \frac{4}{s} ds = 4 \ln s \bigg\rvert_{s = 1}^x = 4 \ln x \]

Plug this into our solution

\[\begin{split} \begin{aligned} y(x) & = e^{-h} \left ( \int_{x_0}^x e^h r ds + y(x_0) \right) \\ & = e^{-4 \ln x} \left( \int_1^s e^{4 \ln s} 8s^3 ds + 2 \right) \\ & = \frac{1}{x^4} \left( 8 \int_1^x s^7 ds + 2 \right) \\ & = \frac{1}{x^4} \left( s^8 \bigg\rvert_{s=1}^x + 2 \right) \\ & = x^4 + \frac{1}{x^4} \end{aligned} \end{split}\]