Linear first-order equations and integrating factors

Linear equations

A first-order equation is linear when $y$ and $y'$ appear only on their own, each multiplied by a function of $t$ alone. Every such equation can be written with $y'$ first, with coefficient 1:

Bring an equation to standard form before reading off $p$ and $q$:

• $y' = 0.8(2 - y)$ from lesson 6.1.1 becomes $y' + 0.8y = 1.6$, so $p = 0.8$ and $q = 1.6$.
• $y' = t + y$ becomes $y' - y = t$. Moving $y$ to the left flips its sign: $p = -1$.
• $2y' + 4y = 6$ becomes $y' + 2y = 3$ after dividing by 2.

Squares and products of $y$ break linearity. $y' = y^2$, the logistic $y' = y\left(1 - \frac{y}{4}\right)$, $y' = \sin y$ and $y\,y' = t$ are not linear.

When $q = 0$, the equation $y' + p(t)\,y = 0$ is homogeneous. It is separable as well (lesson 6.1.1): $\frac{dy}{y} = -p\,dt$ gives $y = Ce^{-\int p\,dt}$. When $q \neq 0$ it usually does not separate, and needs a new idea.

The integrating factor

The product rule (lesson 2.3.2) gives $(\mu y)' = \mu y' + \mu' y$ for any function $\mu(t)$, written with the Greek letter mu. Multiply the standard form by $\mu$:

The left side is $(\mu y)'$ exactly when $\mu' = p\,\mu$. That is a separable equation for $\mu$, solved by $\mu = e^{\int p\,dt}$, since $\frac{d}{dt}e^{g(t)} = g'(t)\,e^{g(t)}$ (lesson 2.3.4). This $\mu$ is the integrating factor.

A constant $p$

Solve $y' + 2y = 6$ with $y(0) = 0$.

1. 1
   $$\mu = e^{2t}$$

   Why?
2. 2
   $$(e^{2t}y)' = 6e^{2t}$$

   Why?
3. 3
   $$e^{2t}y = 3e^{2t} + C$$

   Why?
4. 4
   $$y = 3 - 3e^{-2t}$$

   Why?

The equation is also $y' = 2(3 - y)$, the equilibrium case of lesson 6.1.1 with $k = 2$ and $M = 3$, and separation gives the same answer. Whenever $p$ and $q$ are constants, the solution is $y = \frac{q}{p} + \left(y_0 - \frac{q}{p}\right)e^{-pt}$.

A forcing that changes: $y' + y = t$

With $y(0) = 1$, the integrating factor is $e^{t}$, so $(e^t y)' = t\,e^t$. Integration by parts (lesson 2.6.2) gives $\int t\,e^t\,dt = (t - 1)e^t + C$.

The start $y(0) = 1$ gives $-1 + C = 1$, so $C = 2$ and $y = t - 1 + 2e^{-t}$. This equation does not separate: $y' = t - y$ is a difference, not a product.

A steady part plus a decaying part

Every solution of $y' + y = t$ is $t - 1$ plus $Ce^{-t}$. The first part is the same for every start. The second carries the start, and it dies away: by $t = 6$ it is 0.004958 of the 5.004958 in the table.

With a constant $p > 0$, the start is forgotten at the rate $p$. Every $\frac{1}{p}$ units of time, the gap to the steady part shrinks by a factor of $e$. For $y' + 2y = 6$ the steady part is the equilibrium 3; for $y' + y = t$ it is a moving line.

When the second part grows

A negative $p$ turns decay into growth. $y' = t + y$ is $y' - y = t$, with $p = -1$ and $\mu = e^{-t}$. Parts gives $\int t\,e^{-t}\,dt = -(t + 1)e^{-t}$, so $y = -t - 1 + Ce^{t}$.

From $y(0) = 1$, $C = 2$ and $y = 2e^t - t - 1$: 3.436564 at $t = 1$, 11.778112 at $t = 2$ and 104.1963 at $t = 4$. Only $C = 0$, the start $y(0) = -1$, stays on the line. Like the equilibrium 0 of $y' = 0.5y$, the line is unstable.

Because the equation is linear, the difference of any two solutions solves the homogeneous equation. Module 6.2 builds the general solutions of second-order equations from the same fact.

A $p$ that changes with $t$

Solve $y' + 2ty = 2t$ with $y(0) = 3$. Now $\int 2t\,dt = t^2$, so $\mu = e^{t^2}$ and $(e^{t^2}y)' = 2t\,e^{t^2}$.

The substitution $u = t^2$, $du = 2t\,dt$ (lesson 2.6.1) gives $\int 2t\,e^{t^2}\,dt = e^{t^2} + C$. So $e^{t^2}y = e^{t^2} + C$, then $y = 1 + Ce^{-t^2}$, and the start fixes $C = 2$.

The steady part is 1. The decaying part $2e^{-t^2}$ dies faster than any exponential $e^{-kt}$, because $p = 2t$ keeps growing. This equation is also $y' = 2t(1 - y)$, which separates; both methods give the same answer.

Where $p$ breaks down

Take $y' + \frac{y}{t} = 1$ with $y(1) = 1$. Here $p = \frac{1}{t}$ has no value at $t = 0$, so work on $t > 0$. There $\int \frac{1}{t}\,dt = \ln t$ (lesson 2.5.2), and $\mu = e^{\ln t} = t$.

So $(ty)' = t$, then $ty = \frac{t^2}{2} + C$, and $y(1) = 1$ gives $C = \frac{1}{2}$: $y = \frac{t}{2} + \frac{1}{2t}$. It is 1.25 at $t = 2$ and 5.05 at $t = 0.1$, blowing up as $t$ approaches the break in $p$. The decaying part $\frac{1}{2t}$ fades only like $\frac{1}{t}$, slower than an exponential.

Case study 6.1.4 meets a homogeneous equation with a changing $p$, $y' = -(0.03 + 0.01t)\,y$, as the price of a bond.

Solving numerically with solve_ivp

An integrating factor needs $\int \mu q\,dt$, which often has no formula, and $p$ or $q$ may come only as data. Then the equation is solved numerically. SciPy's solve_ivp(fun, t_span, y0) does it far better than Euler's method of lesson 6.1.1.

• fun(t, y) returns the slope $f(t, y)$, with t first.
• t_span is the interval, such as (0, 6), and y0 the list of starting values. The solver also handles several equations at once, so one value still goes in a list.
• The answer sol.y has one row per unknown, so sol.y[0] holds $y$ at the times sol.t. t_eval asks for chosen times.

It chooses its own steps. Its default method, RK45, estimates each step's error and shrinks the step until that estimate is small enough. Where the solution is smooth, it lengthens the step again.

For $y' = t - y$ with $y(0) = 1$, the exact value at $t = 6$ is 5.004958. Tighter tolerances cost more evaluations of fun, but far fewer than Euler's method needs:

With the tight tolerances, the values at $t$ = 1, 2, 4 and 6 match the table above to all six decimals.

When atol decides the answer

Take $y' = -10y$ with $y(0) = 1$, whose exact value at $t = 3$ is $e^{-30} \approx 9.36 \times 10^{-14}$. With the defaults, solve_ivp returns about $2.12 \times 10^{-7}$, two million times too large.

Nothing failed: the solution fell below $\text{atol} = 10^{-6}$ by $t \approx 1.38$, and from then on any value that small counts as accurate enough. Even $\text{atol} = 10^{-12}$ leaves the answer off by about half, because $10^{-13}$ is still below it. With $\text{atol} = 10^{-25}$, far below the solution, the relative error drops to about $7 \times 10^{-8}$, as $\text{rtol} = 10^{-8}$ roughly asks.

Common mistakes

1. Use the standard form first. $y' = 3y + 1$ is $y' - 3y = 1$, so $p = -3$; $2y' + 4y = 6$ must be divided by 2.
2. The factor is $e^{\int p\,dt}$, not $e^{p}$ or $\int p\,dt$. For $p = 2t$ it is $e^{t^2}$, not $e^{2t}$.
3. Multiply the forcing too. The right side becomes $\mu q$, not $q$.
4. Add $C$ before dividing by $\mu$. $y' + y = t$ gives $t - 1 + Ce^{-t}$, not $t - 1 + C$.
5. The steady part is not the whole answer. The decaying part carries the start, and when $p < 0$ it grows instead.
6. solve_ivp wants fun(t, y) in that order, y0 as a list, and returns $y$ as sol.y[0].
7. The default tolerances give about three digits, and none at all once $|y|$ falls below atol. Set rtol, and atol well below the smallest size of $y$ you care about.