6 Solutions of ordinary differential equations

He felt a restless, vague ambition
A craving for a change of air
(A most unfortunate condition,
A cross not many choose to bear.)
– Alexander Pushkin, Eugene Onegin

In the last chapter we considered discrete time models, in which time is counted in integers. This worked well to describe processes that happen in periodic cycles, like cell division or heart pumping. Many biological systems do not work this way. Change can happen continuously, that is, at any moment in time. For instance, the concentration of a biological molecule in the cell changes gradually, as does the voltage across the cell membrane in a neuron.

There are at least two good reasons to use differential equations for many applications. First, they are often more realistic than discrete time models, because some events happen very frequently and non-periodically. The second reason is mathematical: it turns out that dynamical systems with continuous time, described by differential equations, are better behaved than difference equations. This has to do with the essential “jumpiness” of difference equations. Even for simple nonlinear equations, the value of the variable after one time step can be far removed from its last value. This can lead to highly complicated solutions, as we saw in some of the numerical solutions in the last chapter.

That kind of erratic behavior is impossible in (reasonable) differential equations. Solutions of ODEs are generally smooth and predictable, as long they are defined by continuous defining function, the flow produces smooth solutions. Sometimes, they can be written down exactly, as a mathematical function, called an analytic solution. For most ODEs, particularly nonlinear ones, that is not possible, but then computational methods can be used to produce an approximate numeric solution that can visualize the behavior of the solutions.

build differential equations based on stated assumptions
find and verify analytic solutions of linear differential equations
understand the concepts of numeric solutions of differential equations

6.1 Building differential equations

6.1.1 from discrete time to continuous

In this chapter we investigate continuous time dynamical systems, for which it does not make sense to break time up into equal intervals. Instead of equations describing the increments in the dependent variable from one time step to the next, we will see equations with the instantaneous rate of the change (derivative) of the variable. Let us see the connection between the discrete and continuous dynamic models by reducing the step size of the bacteria-division population model.

First, suppose that instead of dividing every hour, the population of bacteria divide every half-hour, but only half of the population does. That half is chosen randomly, so we don’t have to keep track of whether each bacterium divided the last time around or not. Therefore, each half-hour exactly half of the population is added to the current population: \[ N(t+0.5) = N(t) + 0.5N(t) = 1.5N(t)\] The solution for this model can be figured out from the linear difference equation solution we derived in section $\ref{sec:math14}$ Every half-hour, the population is multiplied by 1.5, so we can write: \[ N(t) = 1.5^{2t} N(0) = (1.5^2)^t N(0)\] Compare this solution with the one for the every-hour model, $N(t) = 2^t N(0)$ by plugging in a few numbers for $t$. The half-hour model grows faster, because it has the base of 2.25 instead of 2.

Now, suppose that the bacteria can divide four times an hour, but only a quarter of the population reproduces at any given time. The model can be written similarly: \[N(t+0.25) = N(t) + 0.25N(t) = 1.25N(t)\] The solution for this model is once again exponential, with the difference that each half contains 4 division events: \[N(t) = 1.25^{4t}N(0) = (1.25^4)^t N(0)\] This solution has the exponential base is $1.25^4$, which is larger than $1.5^2$. So what happens when we take this further?

Suppose the bacteria divide $m$ times an hour, with time step $1/m$. Then extending our models above, we can write down the model and the solution: \[ N(t + 1/m) = N(t) + 1/m N(t) = (1+1/m) N(t) \] \[ N(t) = (1+1/m)^{mt} N(0) = [(1+1/m)^m]^t N(0)\]

Now we can do what mathematicians enjoy the most: take things to the limit. What if $m$ were 100? A million? A gazillion? Let us re-write the model equation: \[ N(t+1/m) - N(t) = 1/m N(t) \Rightarrow \frac{N(t+1/m) - N(t) }{1/m} = N(t)\]

The expression on the left is known as Newton’s quotient that you encounter in the definition of a derivative. It measures the rate of change of the population $N$ from some time $t$ to the next time step $t+1/m$. If $m$ is increased to make the time step smaller, this makes both the numerator and the denominator smaller, and the quotient approaches the instantaneous rate of change of $N(t)$. So, if bacteria divide at any point in time, with the , the model becomes a differential equation: \[ \frac{dN}{dt} = N(t)\]

We can do a similar procedure to the formula of the solution of the model. The dependence on $m$ is all on the left-hand side, in the expression $(1+1/m)^m$, which is the base of the exponential function. What happens to this number as $m$ becomes larger? Does it increase without bound? You can investigate this numerically by plugging in progressively larger numbers $m$, and see that the number approaches a specific value: 2.71828… This is the special constant $e$, called the base of the natural logarithm. So, if bacteria divide at any point in time, with the average rate of 1 per hour, the solution of the model becomes: \[ N(t) = e^t N(0)\]

6.1.2 Exercises (optional)

Here we will explore the effect of changing the step size on the solution of a discrete time dynamic model. We will use a very simple model of bacterial population growth, in which we assume that bacteria divide once an hour and there are no deaths.

Calculate the solution for this population, assuming that all bacteria divide exactly once an hour - in other words, a birth rate of one per individual. Starting with one bacterium use a for loop to calculate the solution for 10 hours and print out the last value.

# YOUR CODE HERE

Suppose that these bacteria can divide twice an hour, but only half of the population divides each time - in other words, a per capita birth rate of 0.5 per half an hour. Change your model so it calculates a solution vector with the time step of 30 minutes over 10 hours, print out the number of bacteria after 10 hours and compare it with the previous value.

# YOUR CODE HERE

Suppose that these bacteria divide every 15 minutes, but only one quarter of the population divides each time - in other words, a per capita birth rate of 0.25 per quarter hour. Change your model so it calculates a solution vector with the time step of 15 minutes over 10 hours, print out the number of bacteria after 10 hours and compare it with the previous value.

# YOUR CODE HERE

Suppose that these bacteria divide every 1 minute, but only 1/60 of the population divides each time - in other words, a per capita birth rate of 1/60 per minute. Change your model so it calculates a solution vector with the time step of 1 minute over 10 hours, print out the number of bacteria after 10 hours and compare it with the previous value.

# YOUR CODE HERE

Suppose that these bacteria divide every second, but only 1/3600 of the population divides each time - in other words, a per capita birth rate of 1/3600 per second. Change your model so it calculates a solution vector with the time step of 1 second over 10 hours, print out the number of bacteria after 10 hours and compare it with the previous value.

# YOUR CODE HERE

Produce a plot of the five solutions of bacterial population dividing with different time steps. Take the five code chunks from above, and copy them all into the chunk below. For each calculation add a time vector that corresponds to each time step (e.g. one for every hour for the first one, one for every second for the last one) and make a plot of each of the solutions as function of time on the same plot - use plot() for the first one and lines() for all the rest, with different colors and add a legend indicating different time steps.

# YOUR CODE HERE

What behaviors do you see for the solutions with different time steps? What effect does shrinking the time step have on the solution? What do you expect would happen if the time step were a millisecond, or a microsecond?

6.1.3 growth proportional to population size

We will now build some common differential equations models. First, a simple population growth model with a constant growth rate. Suppose that in a population each individual reproduces with the average reproductive rate $r$. This is reflected in the following differential equation:

\[\begin{equation} \frac{d x} {dt} = \dot x = r x \label{eq:linear_ode} \end{equation}\]

This expression states that the rate of change of $x$, which we take to be population size, is proportional to $x$ with multiplicative constant $r$. We will sometimes use the notation $\dot x$ for the time derivative of $x$ (which was invented by Newton) for aesthetic reasons.

First, we apply dimensional analysis to this model. The units of the derivative are population per time, as can be deduced from the Newton’s quotient definition. Thus, the units in the equation have the following relationship: \[ \frac{[population]}{[time]} = [r] [population] = \frac{1}{[time]}[population] \]

This shows that as in the discrete time models, the dimension of the population growth rate $r$ is inverse time, or frequency. The difference with the discrete time population models lies in the time scope of the rate. In the case of the difference equation, $r$ is the rate of change per one time step of the model. In the differential equation, $r$ is the instantaneous rate of population growth. It is less intuitive than the growth rate per single reproductive cycle, just like the slope of a curve is less intuitive than the slope of a line. The population growth happens continuously, so the growth rate of $r$ individuals per year does not mean that if we start with one individual, there will be $r$ after one year. In order to make quantitative predictions, we need to find the solution of the equation, which we will see in the next section.

6.1.4 chemical kinetics

Reactions between molecules in cells occur continuously, driven by molecular collisions and physical forces. In order to model this complex behavior, it is generally assumed that reactions occur with a particular speed, known as the kinetic rate constant. As mentioned in chapter 2, a simple reaction of conversion from one type of molecule ($A$) to another ($B$) can be written as follows: \[ A \xrightarrow{k} B \] In this equation the parameter $k$ is the kinetic rate rate constant, describing the speed of conversion of $A$ into $B$, per concentration of $A$.

Chemists and biochemists use differential equations to describe the change in molecular concentration during a reaction. These equations are known as the laws of mass action. For the reaction above, the concentration of molecule $A$ decreases continuously proportionally to itself, and the concentration of molecule $B$ increases continuously proportionally to the concentration of $A$. This is expressed by the following two differential equations:

\[\begin{eqnarray} \label{eq:lin_chem_kin} \dot A &=& - k A \\ \dot B &=& kA \end{eqnarray}\]

Several conclusions are apparent by inspection of the equations. First, the dynamics depend only on the concentration of $A$, so keeping track of the concentration of $B$ is superfluous. The second observation reinforces the first: the sum of the concentrations of $A$ and $B$ is constant. This is mathematically demonstrated by adding the two equations together to obtain the following:

\[ \dot A + \dot B = -kA + kA = 0 \]

One of the basic properties of the derivative is that the sum of derivatives is the same as the derivative of the sum:

\[ \dot A + \dot B = \frac{d(A+B)}{dt} = 0 \]

This means that the sum of the concentrations of $A$ and $B$ is a constant. This is a mathematical expression of the law of conservation in chemistry: molecules can change from one type to another, but they cannot appear or disappear in other ways. In this case, a single molecule of $A$ becomes a single molecule of $B$, so it follows that the sum of the two has to remain the same. If the reaction were instead two molecules of $A$ converting to a molecule of $B$, then the conserved quantity is $2A + B$. The concept of conserved quantity is very useful for the analysis of differential equations. We will see in later chapters how it can help us find solutions, and explain the behavior of complex dynamical systems.

6.1.5 building nonlinear ODEs

The simple, linear population growth models we have seen in the last two chapters assume that the per capita birth and death rates are constant, that is, they stay the same regardless of population size. The solutions for these models either grow or decay exponentially, but in reality, populations do not grow without bounds. It is generally true that the larger a population grows, the more scarce the resources, and survival becomes more difficult. For larger populations, this could lead to higher death rates, or lower birth rates, or both.

How can we incorporate this effect into a quantitative model? We will assume there are separate birth and death rates, and that the birth rate declines as the population grows, while the death rate increases. Suppose there are inherent birth rates $b$ and $d$, and the overall birth and death rates $B$ and $D$ depend linearly on population size $P$: $B = b - aP$ and $D = d + cP$.

To model the rate of change of the population, we need to multiply the rates $B$ and $D$ by the population size $P$, since each individual can reproduce or die. Also, since the death rate $D$ decreases the population, we need to put a negative sign on it. The resulting model is:

\[ \dot P = BP - DP = [(b-d)-(a+c)P]P \]

The parameters of the model, the constants $a,b,c,d$, have different meanings. Performing dimensional analysis, we find that $b$ and $d$ have the dimensions of $1/[t]$, the same as the rate $r$ in the exponential growth model. However, the dimensions of $a$ (and $c$) must obey the relation: $[P]/[t] = [a][P]^2$, and thus,

\[[a]=[c] = \frac{1}{[t][P]}\]

This shows that the constants $a$ and $c$ have to be treated differently than $b$ and $d$. Let us define the inherent growth rate of the population, to be $r_0=b-d$ (if the death rate is greater than the birth rate, the population will inherently decline). Then let us introduce another constant $K$, such that $(a+c)=r_0/K$. It should be clear from the dimensional analysis that $K$ has units of $P$, population size. Now we can write down the logistic equation in the canonical form:

\[\begin{equation} \dot P = r\left(1-\frac{P}{K}\right)P \label{eq:log_cont_model} \end{equation}\]

This model can be re-written as $\dot P = aP -bP^2$, so it is clear that there is a linear term ($aP$) and a nonlinear term ($-bP^2$). When $P$ is sufficiently small (and positive) the linear term is greater, and the population grows. When $P$ is large enough, the nonlinear term wins and the population declines.

It should be apparent that there are two fixed points, at $P=0$ and at $P=K$. The first one corresponds to a population with no individuals. On the other hand, $K$ signifies the population at which the negative effect of population size balances out the inherent population growth rate, and is called the carrying capacity of a population in its environment . We will analyze the qualitative behavior of the solution, without writing it down, in the next section of this chapter.

6.2 Solutions of ordinary differential equations

In this section we will investigate how to write down analytic solutions for ordinary differential equations (ODEs). Let us first define the mathematical terms that will be used in this discussion.

Definition

An ordinary differential equation is an equation that contains derivatives of the dependent variable (e.g. $x$) with respect to an independent variable (e.g. $t$). For example: \[ \frac{dx^2}{dt^2}+ 0.2 \frac{dx}{dt} - 25 = 0 \]

For the time being, we will restrict ourselves to ODEs with the highest derivative being of first order, called first-order ODEs. Second-order ODEs are common in models derived from physics, but can actually be converted to first-order ODEs, though it requires an additional dependent variable. Third and higher order ODEs are very uncommon. To be precise:

Definition

A first-order ODE is one where the derivative $dx/dt$ is equal to a defining function $f(x,t)$, like this: \[\frac{dx} {dt} = \dot x = f(x,t)\]

The defining function may potentially depend on both the dependent variable $x$ and the independent variable $t$. If it only depends on $x$, it is called an autonomous ODE, for example:

\[ \frac{dx}{dt} = \dot x = rx \] or

\[\frac{dx}{dt} = \dot x = 5x -4\]

On the other hand, if the defining function depends only on the independent variable $t$, it may be called a pure-time ODE, for example:

\[ \frac{dx}{dt} = \dot x = 5t \] or

\[ \frac{dx}{dt} = \dot x = 20 - 0.3 \sin(4 \pi t) \]

An ODE is if every term involves either the dependent variable $x$ or its derivative. \end{mosdef} For instance, $\dot x = x^2 + \sin(x)$ is homogeneous, while $\dot x = -x + 5t$ is not.

Most simple biological models that we will encounter in the next two chapters are autonomous, homogeneous ODEs. However, nonhomogeneous equations are important in many applications, and we will encounter them at the end of the present section.

6.2.1 separate and integrate method (optional)

Definition

The analytic (or exact) solution of an ordinary differential equation is a function of the independent variable that satisfies the equation. If no initial value is given, then the general solution function will contain an unknown integration constant. If an initial value is specified, the integration constant can be found to obtain a specific solution.

This means that the solution function obeys the relationship between the derivative and the defining function that is specified by the ODE. To verify that a function is a solution of a given ODE, take its derivative and check whether it matches the other side of the equation.

Example. The function $x(t) = 3t^2 +C$ is a general solution of the ODE $\dot x = 6t$, which can be verified by taking the derivative: $\dot x (t) = 6t$. Since this matches the right-hand side of the ODE, the solution is valid.

Example. The function $x(t) = Ce^{5t}$ is a general solution of the ODE $\dot x = 5x$. This can be verified by the taking the derivative: $\dot x = 5C e^{5t}$ and comparing it with the right-hand side of the ODE: $5x = 5 Ce^{5t}$. Since the two sides of the equation agree, the solution is valid.

In contrast with algebraic equations, we cannot simply isolate $x$ on one side of the equal sign and find the solutions as one, or a few numbers. Instead, solving ordinary differential equations is very tricky, and no general strategy for solving an arbitrary ODE exists. Moreover, a solution for an ODE is not guaranteed to exist at all, or not for all values of $t$. We will discuss some of the difficulties later, but let us start with equations that we can solve.

The most obvious strategy for solving an ODE is integration. Since a differential equation contains derivatives, integrating it can remove the derivative. In the case of the general first order equation, we can integrate both sides to obtain the following: \[ \int \frac{dx}{dt} dt = \int f(x,t) dt \Rightarrow x(t) + C = \int f(x,t) dt\] The constant of integration $C$ appears as in the standard antiderivative definition. It can be specified by an initial condition for the solution $x(t)$. Unless the function $f(x,t)$ depends only on $t$, it is not possible to evaluate the integral above. Instead, various tricks are used to find the analytic solution. The simplest method of analytic solution of a first-order ODEs, which I call separate-and-integrate consists of the following steps:

use algebra to place the dependent and independent variables on different sides of the equations, including the differentials (e.g. $dx$ and $dt$)
integrate both sides with respect to the different variables, don’t forget the integration constant
solve for the dependent variable (e.g. $x$) to find the general solution
plug in $t=0$ and use the initial value $x(0)$ to solve for the integration constant and find the the specific solution

Example. Consider a very simple differential equation: $\dot x = a$, where $\dot x$ stands for the time derivative of the dependent variable $x$, and $a$ is a constant. It can be solved by integration: \[ \int \frac{dx}{dt} dt = \int a dt \Rightarrow x(t) + C = at \]

This solution contains an undetermined integration constant; if an initial condition is specified, we can determine the complete solution. Generally speaking, if the initial condition is $x(0) = x_0$, we need to solve an algebraic equation to determine $C$: $x_0 = a*0 - C$, which results in $C = -x_0$. The complete solution is then $x(t) = at + x_0$. To make the example more specific, if $a = 5$ and the initial condition is $x(0) = -3$, the solution is

\[x(t) = 5t -3\]

Example. Let us solve the linear population growth model in equation $\ref{eq:linear_ode}$: $\dot x = rx$. The equation can be solved by first dividing both sides by $x$ and then integrating: \[ \int \frac{1}{x} \frac{d x}{dt} dt = \int \frac{dx}{x} = \int r dt \Longrightarrow \log |x| = rt + C \Longrightarrow x = e^{rt+C} = Ae^{rt}\]

We used basic algebra to solve for $x$, exponentiating both sides to get rid of the logarithm on the left side. As a result, the additive constant $C$ gave rise to the multiplicative constant $A=e^C$. Once again, the solution contains a constant which can be determined by specifying an initial condition $x(0) = x_0$. In this case, the relationship is quite straightforward: $x(0) = A e^0 = A$. Thus, the complete solution for equation $\ref{eq:linear_ode}$ is: \[ x(t) = x_0e^{rt}\]

6.2.2 behavior of solutions of linear ODEs

As in the case of the discrete-time models, population growth with a constant birth rate has exponential form. Once again, please pause and consider this fact, because the exponential solution of linear equations is one of the most basic and powerful tools in applied mathematics. Immediately, it allows us to classify the behavior of linear ODE into three categories:

$r > 0$: $x(t)$ grows without bound
$r < 0$: $x(t)$ decays to 0
$r = 0$: $x(t)$ remains constant at the initial value

The rate $r$ being positive reflects the dominance of birth rate over death rate in the population, leading to unlimited population growth. If the death rate is greater, the population will decline and die out. If the two are exactly matched, the population size will remain unchanged.

Example. The solution for the biochemical kinetic model in equation $\ref{eq:lin_chem_kin}$ is identical except for the sign: $A(t) = A_0 e^{-kt}$. When the reaction rate $k$ is positive, as it is in chemistry, the concentration of $A$ decays to 0 over time. This should be obvious from our model, since there is no back reaction, and the only chemical process is conversion of $A$ into $B$. The concentration of $B$ can be found by using the fact that the total concentration of molecules in the model is conserved. Let us call it $C$. Then $B(t) = C - A(t) = C- A_0e^{-kt}$. The concentration of $B$ increases to the asymptotic limit of $C$, meaning that all molecules of $A$ have been converted to $B$.

6.2.3 solutions of nonhomogeneous ODEs

ODEs that contain at least one term without the dependent variable are a bit more complicated. If the defining function is $f(x,t)$ is linear in the dependent variable $x$, they can be solved on paper using the same separate-and-integrate method, modified slightly to handle the constant term. Here are the steps to solve the generic linear ODE with a constant term $\dot x = ax +b$:

separate the dependent and independent variables on different sides of the equations, by dividing both sides by the right hand side $ax+b$, and multiplying both sides by the differential $dt$
integrate both sides with respect to the different variables, don’t forget the integration constant!
solve for the dependent variable (e.g. $x$)
plug in $t=0$ and use the initial value $x(0)$ to solve for the integration constant

Example. Let us solve the following ODE model using separate and integrate with the given initial value: \[\frac{dx}{dt} = 4x -100; \; x(0) = 30\]

separate the dependent and independent variables: \[ \frac{dx}{4x - 100} = dt\]
integrate both sides: \[ \int \frac{dx}{4x -100} = \int dt \Rightarrow \frac{1}{4} \int \frac{du}{u} = \frac{1}{4} \ln | 4x- 100 | = t + C\] The integration used the substitution of the new variable $u=4x -100$, with the concurrent substitution of $dx = du/4$.
solve for the dependent variable: \[ \ln | 4x- 100 | = 4t + C \Rightarrow 4x-100 = e^{4t} B \Rightarrow x = 25 + Be^{4t}\] Here the first step was to multiply both sides by 4, and the second to use both sides as the exponents of $e$, removing the natural log from the left hand side, and finally simple algebra to solve for $x$ as a function of $t$.
solve for the integration constant: \[ x(0) = 25 + B = 30 \Rightarrow B = 5\] Here the exponential “disappeared” because $e^0=1$.

Therefore, the complete solution of the ODE with the given initial value is \[x(t) = 25 + 5e^{4t}\]

At this point, you might have noticed something about solutions of linear ODEs $dx/dt = ax+b$: they always involve an exponential term, with time in the exponent. This can be verified by plugging the solution of the form $x(t) = Ce^{at} + D$ into the ODE to see if it satisfies the equation. First, take the derivative of the solution to get the left-hand side of the ODE: $\frac{dx}{dt} = Ca e^{at}$; the plug in $x(t) = $ into the right hand side of the ODE: $aCe^{at} + aD +b$. Setting the two sides equal, we get: $Ca e^{at} = aCe^{at} + aD +b$, which is satisfied if $aD + b = 0$, which means $D= -b/a$. This is consistent with the example above, the additive constant in the solution was 25, which is $-b/a= -(-100)/4 = 25$.

Important fact

Any linear ODE of the form: \[ \frac {dx}{dt} = ax +b \] has an analytic solution of the form: \[ x(t) = Ce^{at} - \frac{b}{a} \]

The unknown constant $C$ can be determined from a given initial value. So the upshot is that all linear ODEs have solutions which are exponential in time with exponential constant coming from the slope constant $a$ in the ODE. The dynamics of the solution are determined by the sign of the constant $a$: if $a>0$, the solution grows (or declines) without bound; and if $a<0$, the solution approaches an asymptote at $-b/a$ (from above or below, depending on the initial value). Go back and read section $\ref{sec:math2}$ for a review of exponential functions if this is not clear.

6.2.4 Exercises

Solve the following linear ODEs and use the specified initial values to determine the integration constant. Describe how the solution behaves over a long time (e.g. grows without bound, goes to zero, etc.). Plug the solution back into the ODE to check that it satisfies the equation.

\[ \frac{dx}{dt} = 0.1; \; x(0)= 100 \]
\[ \frac{dx}{dt} = 2\sin(4t) -0.4t; \; x(0)= 5 \]
\[ \frac{dx}{dt} = 3x; \; x(0) = 0.4 \]
\[ \frac{dx}{dt} = -5x; \; x(0) = -300 \]
\[ \frac{dx}{dt} = -0.5x + 100 ; \; x(0) = 20 \]
\[ \frac{dx}{dt} = 1 + x; \; x(0) = 4 \]
\[ \frac{dx}{dt} = -10 - 0.2x; \; x(0) = 10 \]
\[ \frac{dx}{dt} = -4 + 0.5x; \; x(0) = 6 \]

For the following ODE models, verify whether a particular function is a solution:

For the ODE model of concentration $C(t)$:

\[ \frac{dC}{dt} = -rC + A \]

Check whether the function $C(t) = (-r)^t + B$ satisfies the ODE, and if yes, how the parameter B is determined by model parameters r and A.
For the ODE model of membrane potential $V(t)$ with:

\[ \frac{dV}{dt} = -a(V-V_1) - b(V-V_2) \]

Check whether the function $V(t) = Ce^{-(a+b)t} + D$ satisfies the ODE, and if yes, how the parameters C and D are determines by model parameters a, b, $V_1$ and $V_2$.
For the logistic ODE:

\[ \frac{dP}{dt} = r(1-P/K)P \] Check whether the function $P(t) = A/(1+Be^{-rt})$ is a solution of the logistic ODE and relate the parameters A and B in the function to the parameters $r$ and $K$ in the ODE.
For the logistic ODE with harvesting:

\[ \frac{dP}{dt} = r(1-P/K)P - \mu \] Check whether the function $P(t) = A/(1+Be^{-rt}) - \mu t$ is a solution of the logistic ODE and relate the parameters A and B in the function to the parameters $r$ and $K$ in the ODE.

6.3 Numeric solutions and the Forward Euler method (optional)

Analytic solutions are very useful for a modeler because they allow prediction of the variable of interest at any time in the future. However, for many differential equations they are not easy to find, and for many others they simply cannot be written down in a symbolic form. Instead, one can use a numerical approach, which does not require an exact formula for the solution. The idea is to start at a given initial value (e.g. $x(0)$) and use the derivative from the ODE (e.g. $dx/dt$) as the rate of change of the solution (e.g. $x(t)$) to calculate the change or increment for the solution over a time step. Essentially, this means replacing the continuous change of the derivative with a discrete time step, thus converting the differential equation into a difference equation and then solving it. The solution of the difference equation is not the same as the solution of the ODE, so numeric solutions of ODEs are always approximate. I will use the letter $y(t)$ to denote the numerical solution to distinguish it from the exact solution $x(t)$. The fundamental difference between them is that $y(t)$ is not a formula that can be evaluated at any point in time, but instead is a sequence of numbers calculated every time step, which hopefully are close to the exact solution $x(t)$.

Let us introduce all the players: first, we need to pick the time step $\Delta t$, which is the length of time between successive values of $y$. In the difference equation notation one can use $y_i$ to mean $y(i\Delta t)$, the value of the numerical solution after $i$ time steps. Then we need to calculate the derivative, or the rate of change at a particular point in time. For any first-order ODE of the form \[ \frac{d x} {dt} = \dot x = f(x,t)\]

the rate of change depends (potentially) on the values of $x$ and $t$. This rate of change based on the numerical solution after $i$ time steps is $f(y(i\Delta t), i\Delta t) = f(y_i, t_i)$. Finally, to calculate the change of the dependent variable we need to multiply the rate of change by the time step. This should make sense in a practical context: if you drive for two hours (time step) at 60 miles per hour (rate of change), the total distance (increment) is $2*60=120$ miles. By the same token, we can write down how to calculate the next value of the numerical solution $y_{i+1}$ based on the previous one: \[\begin{equation} y_{i+1} = y_i + \Delta t f(y_i, t_i) \label{eq:ch15_FE} \end{equation}\]

This method of computing a numerical solution of an ODE is called the Forward Euler method, after the famous mathematician who first came up with it. It is called a forward method because it uses the value of the dependent variable and its derivative at time step $i$ to predict the value at the next time step $i+1$. The method is iterative, so it needs to be repeated in order to calculate a set of values of the approximate solution $y(t)$. Here are a couple of simple examples of computing numerical solution using FE:

Let us numerically solve the ODE $\dot x = -0.1$ using the Forward Euler method. This means the defining function in the formulation of FE above is $f(x,t)=-0.1$. We can calculate the numeric solution for a couple of steps and compare the values with the exact solution, since we now know that it is $x(t) = x_0 -0.1t$. Let us pick the time step $\Delta t = 0.2$ and begin with the initial value $x(0)=1$. Here are the first three steps using the FE method: \[ y(0.2) = y(0) + \Delta t f(y(0)) = 1 + 0.2*(-0.1) = 0.98\] \[ y(0.4) = y(0.2) + \Delta t f(y(0.2)) = 0.98+ 0.2*(-0.1) = 0.96\] \[ y(0.6) = y(0.4) + \Delta t f(y(0.4)) = 0.96+ 0.2*(-0.1) = 0.94\] Since the rate of change in this ODE is constant, the solution declines by the same amount every time step. In this case, the numerical solution is actually exact, and perfectly matches the analytic solution. Table $\ref{tab:ch15_FE}$ (right) shows the numerical solution for 3 time steps along with the exact solution.

Let us numerically solve the ODE $\dot x = -0.1x$ using the Forward Euler method. This means the defining function in the formulation of FE above is $f(x,t)=-0.1x$. We can calculate the numeric solution for a couple of steps and compare the values with the exact solution, since we now know that it is $x(t) = x_0 e^{-0.1t}$. Let us pick the time step $\Delta t = 0.2$ and begin with the initial value $x(0)=100$. Here are the first three steps using the FE method: \[ y(0.2) = y(0) + \Delta t f(y(0)) = 100 + 0.2*(-0.1*100) = 98\] \[ y(0.4) = y(0.2) + \Delta t f(y(0.2)) = 98+ 0.2*(-0.1*98) \] \[ = 96.04\] \[ y(0.6) = y(0.4) + \Delta t f(y(0.4)) = 96.04+ 0.2*(-0.1*96.04) \approx \] \[ \approx 94.12\] In this case, the derivative is not constant and the numerical solution is not exact, which is demonstrated in table $\ref{tab:ch15_FE}$ (left). The error in the numerical solution grows with time, which may be problematic. We will further investigate how to implement the computation of numerical solutions using R in the next section.

Numerical solutions of the ODE $\dot x = -0.1$ using Forward Euler $y$ calculated for 3 steps of size $\Delta t = 0.2$ as well as the exact solution $x$, both rounded to two digits after the decimal, and the error of the numerical solution.
t	x	y
0	1	1
0.2	0.98	0.98
0.4	0.96	0.96
0.6	0.94	0.94

Numerical solution of the ODEs $\dot x = -0.1x$ (right) using Forward Euler $y$ calculated for 3 steps of size $\Delta t = 0.2$ as well as the exact solution $x$, both rounded to two digits after the decimal, and the error of the numerical solution.
t	x	y	error
0	100	100	0
0.2	98.02	98	0.02
0.4	96.08	96.04	0.04
0.6	94.18	94.12	0.06

6.3.1 Exercises

Use the Forward Euler method to solve the following differential equations with time step $\Delta t=0.5$ for 2 steps to compute $y(1)$ (the value of the numerical solution at $t=1$.)

\[ \frac{dx}{dt} = 0.1; \; x(0)= 100 \]
\[ \frac{dx}{dt} = 2\sin(4t) -0.4t; \; x(0)= 5 \]
\[ \frac{dx}{dt} = 3x; \; x(0) = 0.4 \]
\[ \frac{dx}{dt} = -5x; \; x(0) = -300 \]
\[ \frac{dx}{dt} = -0.5x + 100 ; \; x(0) = 20 \]
\[ \frac{dx}{dt} = 1 + x; \; x(0) = 4 \]
\[ \frac{dx}{dt} = -10 - 0.2x; \; x(0) = 10 \]
\[ \frac{dx}{dt} = -4 + 0.5x; \; x(0) = 6 \]

6.4 Forward Euler method in R (optional)

6.4.1 implementation

In practice, the most common approach to finding solutions for differential equations is using a computer to calculate a numerical solution, for example using the Forward Euler method. This means using a computer program to construct a sequence of values of the dependent variable that approximate the true solution. Below is an outline of the algorithm that can be translated into a programming language, like R, to solve ODEs.

assign the time step dt, length of time Tmax, number of time steps numstep
pre-allocate the vector of numeric solution values y of length numstep+1
assign the initial value for the ODE to the first element of the solution
assign the vector of time values t from 0 to Tmax of length numstep+1
for loop starting at 1 to numstep
- assign the next solution value to be the current solution value plus the time step multiplied by the defining function at the current solution value
plot numeric solution y as a function of time t

To implement the algorithm, one need to know the defining function $f(x,t)$, the initial value, the time step, and the total time. The output is the solution vector $y$, which contains a sequence of values that approximate the solution of the ODE, along with the vector of time values spaced by the time step. Below is an example implementation of the Forward Euler method for the following linear ODE with initial value 1000:

\[ \frac{dx}{dt} = -0.5 x \]

x0 <- 1000
dt <- 0.01 # set time step
Tmax <- 10 # set length of time
numstep <- Tmax/dt # assign number of time steps
pop <- rep(x0, numstep+1) # initialize solution with y0
for (i in 1:numstep) { # do the Euler!
    pop[i+1] <- pop[i]+dt*(-0.5*pop[i])
}
time <- seq(0,Tmax, dt)
plot(time, pop, type='l')

Notice that it is very similar to the script for numerical solution of a difference equation we saw in $\ref{sec:comp14}$ with the major difference being the presence of a time step, whereas in difference equations the time step is always 1. There is one more important point for the implementation: usually one needs to solve the ODE for a particular length of time $T$ with a specified time step $\Delta t$ . This dictates that the required number of iterations be $T/\Delta t$; in other words, for a given time period the number of time steps is inversely proportional to the time step.

6.4.2 Exercises

Consider a slightly different linear ODE:

\[ \frac{dx}{dt} = 0.2 x \]

Calculate the numeric solution of the ODE for one time step using Forward Euler, for time step dt=0.1, starting with initial value x(0) = 5. Answer: 5.1

# YOUR CODE HERE

Write a script to solve the ODE using the Forward Euler method based on the outline above. Set the time step dt=0.1 and report the solution after 100 time steps. Answer: 36.22323

# YOUR CODE HERE

Change the time step to be dt=0.01 and report the solution after 1000 time steps. Answer: 36.87156.

# YOUR CODE HERE

6.4.3 error analysis

Numeric solutions of ODE are always approximate, because they use discrete time steps to approximate continuous change (derivatives). Thus numeric solutions always have error, which is the difference between the exact or analytic solution and the numeric solution. If we know the exact solution of an ODE, we can calculate the error using vector subtraction in R. For the same linear ODE we solved above: \[ \frac{dx}{dt} = -0.5 x \] The analytic solution is $x(t) = x_0 e^{-0.5t}$, where $x_0$ is the initial value. Here is an example of computing the numeric solution (as we did above) and then calculating the analytic solution and plotting it:

x0 <- 1000
dt <- 0.5 # set time step
Tmax <- 10 # set length of time
numstep <- Tmax/dt # assign number of time steps
pop <- rep(x0, numstep+1) # initialize solution with y0
for (i in 1:numstep) { # do the Euler!
    pop[i+1] <- pop[i]+dt*(-0.5*pop[i])
}
time <- seq(0,Tmax, dt)
plot(time, pop, type='l', main = "Numeric and analytic solutions of an ODE") # plot the numeric solution
exact <- x0*exp(-0.5*time) # calculate the exact solution
lines(time, exact, col = 'red') # plot the exact solution
legend('topright', col=c('black', 'red'), lty=1, legend = c('numeric', 'analytic'))

Now we can calculate the error of the numeric solution and plot it:

error <- abs(exact - pop)
plot(time, error, type = 'l', main='Error of the numeric solution')

What is the sources of this error? There are at least two distinct sources of error in numerical solutions: a) and b) . Roundoff error is caused by computers representing real numbers by a finite string of bits on a computer using what is known as a representation. In many programming languages variables storing real numbers can be single or double precision, which typically support 24 and 53 significant binary digits, respectively. Any arithmetic operation involving floating point numbers is only approximate, with an error that depends on the way the numbers are stored in the memory. Truncation error is caused by approximations inherent in numerical algorithms, such as Forward Euler, which represent instantaneous rate of change in an ODE with discrete steps and thus are always a bit off from the true analytic solution.

In practice, roundoff error is not a big concern in contemporary computation for most modelers. Truncation error, on the other hand, can cause big problems, but luckily it is within your control. One can decrease the error in the case of finite difference methods by choosing smaller time steps, or by choosing an algorithm with a higher order of accuracy, which we’ll leave for a more advanced discussion.

Returning specifically to the Forward Euler method, it is called a first-order method because the total error of the solution (after some number of time steps) depends linearly on the time step $\Delta t$. One can show this by using the Taylor expansion of the solution $y(t)$ to derive the forward Euler method, with $\tau(\Delta t)$ representing the truncation error after one time step: \[y(t+\Delta t) = y(t) + \Delta t \frac {dy(t) }{dt} + \tau(\Delta t)\] As you might have learned in calculus, the error remaining after the linear term in the Taylor series is proportional to the the square of the small deviation $\Delta t$. This only describes the error after 1 time step, but since the errors accumulate every time step, the total error after $N$ time steps accumulates $N \tau(\Delta t)$. As we saw in the implementation above, for a given length of time, $N$ is inversely proportional to $\Delta t$. Therefore, the total error is proportional to the $\Delta t$ and so FE is a first-order method.

The exercise above shows that new errors in FE method accumulate in proportion with the time step. The next question is, what happens to these errors over time? Do they grow or dissipate with more iterations? This is known as the stability of a numerical method, and unlike the above question about the order of accuracy, the answer depends on the particular ODE that one needs to solve. Below I show an example of error analysis for a linear ODE:

Example. To numerically solve the equation $\dot x = ax$, we substitute the function $ax$ for the function $f(x,t)$, and obtain the FE approximation for this particular ODE: \[y_{i+1} = y_i + \Delta t a y_i = (1+a\Delta t) y_i\] The big question is what happens to the truncation error: does it grow or decay? To investigate this question, let us denote the error at time $t_i$ , that is the difference between the true solution $x(t_i)$ and the approximate solution $y(t_i)$, by $\epsilon_i$. It follows that $y_i = x_i + \epsilon_i$. Then we can wrote the following difference equations involving the error:

\[y_{i+1} = x_{i+1} + \epsilon_{i+1} = (x_i + \epsilon_i) (1+a\Delta t) = x_i (1+a\Delta t) + \epsilon_i(1+a\Delta t)\] Let us set aside the terns in the equation that involve $x$ (since it is just the equation for forward Euler). The remaining difference equation for $\epsilon$ describes the change in the error: \[\epsilon_{i+1} = \epsilon_i(1+a\Delta t) \]

This states that the error in this numerical solution is repeatedly multiplied by the constant $(1+a\Delta t)$. As we saw in section $\ref{sec:math14}$, this linear difference equation has an exponential solution $\epsilon_n = (1+a\Delta t)^n \epsilon_0$, which decays to 0 if $|1+a\Delta t| < 1$ or grows without bound if $|1+a\Delta t| > 1$. The first inequality is called the stability condition for the FE scheme, since it guarantees that the old errors decay over time. Since $\Delta t >0$, the only way that the left hand side can be less than 1 is if $a<0$. Therefore, the condition for stability of the FE method for a linear ODE: \[|1 + a\Delta t| < 1 \Rightarrow \Delta t < -2/a\]

Thus, if $a>0$, the errors will eventually overwhelm the solution. If $a<0$, if the time step is small enough (less than $-2/a$) then FE is stable. Generally speaking, however, Forward Euler is about the worst method to use for practical numerical solutions of ODEs, due to its low accuracy and to its lack of stability under certain conditions.

6.4.4 Exercises

\[ \frac{dx}{dt} = 0.2 x \]

Calculate the error of the numeric solution of this ODE after one time step with dt=0.1 and initial value x(0) = 5 by subtracting it from the exact (analytic) solution $x(t) = e^{0.2t}x(0)$, with the same initial value. Answer: about 0.001.

# YOUR CODE HERE

Compute the error of the two numeric solution over Tmax = 10 by subtracting the numeric solution vector from the analytic solution calculating over the same time vector and report the mean of that error vector. Answers: for dt=0.1 the mean error is 0.722, for dt=0.01 the mean error is 0.0737.

# YOUR CODE HERE

6.5 Applications of linear ODE models

6.5.1 model of pharmacokinetics

Describing and predicting the dynamics of drug concentration in the body is the goal of pharmacokinetics. Any drug that humans take goes through several stages: first it is administered (put into the body), then absorbed, metabolized (transformed), and excreted (removed from the body) . Almost any drug has a dose at which it has a toxic effect, and most can kill a human if the dose is high enough. Drugs which are used for medical purposes have a therapeutic range, which lies between the lowest possible concentration (usually measured in the blood plasma) that achieves the therapeutic effect and the concentration which is toxic. One of the basic questions that medical practitioners need to know is how much and how frequently to administer a drug to maintain drug concentration in the therapeutic range.

The concentration of a drug is a dynamic variable which depends on the rates of several processes, most directly on the rate of administration and the rate of metabolism. Drugs can be administered through various means (e.g. orally or intravenously) which influences their rate of absorption and thus how the concentration increases. Once in the blood plasma, drugs are metabolized primarily by enzymes in the liver, converting drug molecules into compounds that can be excreted through the kidneys or the large intestine. The process of *metabolism proceeds at a rate that depends on both the concentration of the drug and on the enzyme that catalyzes the reaction. For some drugs the metabolic rate may be constant, or independent of the drug concentration, since the enzymes are already working at full capacity and can’t turn over any more reactions, for example alcohol is metabolized at a constant rate of about 1 drink per hours for most humans. Figure $\ref{fig:ch15_drug_conc}$a shows the time plots of the blood alcohol concentration for 4 males who ingested different amounts of alcohol, and the curves are essentially linear with the same slope after the peak. For other drugs, if the plasma concentration is low enough, the enzymes are not occupied all the time and increasing the drug concentration leads to an increase in the rate of metabolism. One can see this behavior in the metabolism of the anti-depressant drug bupropion in figure $\ref{fig:ch15_drug_conc}$b, where the concentration curve shows a faster decay rate for higher concentration of the drug than for lower concentration. In the simplest case, the rate of metabolism is linear, or proportional to the concentration of the drug, with proportionality constant called the first-order metabolic rate.

Blood alcohol content after ingesting different numbers of drinks, from 4 in the top curve to 1 in the bottom (figure from the National Institute on Alcohol Abuse and Alcoholism in public domain)

Blood concentration of bupropion for two different drugs in clinical trials (image by CMBJ based on FDA data under CC-BY 3.0 via Wikimedia Commons)

Let us build an ODE model for a simplified pharmacokinetics situation. Suppose that a drug is administered at a constant rate of $M$ (concentration units per time unit) and that it is metabolized at a rate proportional to its plasma concentration $C$ with metabolic rate constant $k$. Then the ODE model of the concentration of the drug over time $C(t)$ is: \[ \frac{dC}{dt} = M - kC\] The two rate constants $M$ and $k$ have different dimensions, which you should be able to determine yourself. The ODE can be solved using the separate-and-integrate method:

Divide both sides by the right hand side $M-kC$, and multiply both sides by the differential $dt$ \[ \frac{dC}{M-kC} = dt\]
integrate both sides with respect to the different variables, don’t forget the integration constant! \[ \int \frac{dC}{M-kC} = \int dt \Rightarrow\] \[ -\frac{1}{k} \log |M-kC| = t + A\]
solve for the dependent variable $C(t)$ \[ \exp(\log |M-kC| ) = -\exp(kt +A) \Rightarrow \] \[ M - kC = B e^{-kt} \Rightarrow \] \[ C(t) = \frac{M}{k}- Be^{-kt}\] Notice that I changed the values of integration constants $A$ and $B$ during the derivation, which shouldn’t matter because they have not been determined yet.
plug in $t=0$ and use the initial value $x(0)$ to solve for the integration constant If we know the initial value $C(0) = C_0$, then we can plug it in and get the following algebraic expression: \[ C_0 = \frac{M}{k} - B \Rightarrow \] \[ B = C_0 - \frac{M}{k}\] Then the complete solution is: \[ C(t) = \frac{M}{k} - (C_0- \frac{M}{k})e^{-kt}\]

The solution predicts that after a long time the plasma concentration will approach the value $M/k$, since the exponential term decays to zero. Notice that mathematically this is the same type of solution we obtained in equation $\ref{eq:ch15_ode_sol}$ for a generic linear ODE with a constant term.

6.5.2 Discussion questions

The following questions encourage you to think critically about the pharmacokinetic model above.

Describe in words the dependence of the long-term plasma concentration of the drug on the }parameters. Does this prediction make intuitive sense?
Explain in practical terms the assumption that the administration of the drug results in a constant rate of growth of the concentration. Under what circumstances does this match reality?
Explain in practical terms the assumption that the drug metabolism rate is proportional to the plasma concentration. Under what circumstances does this match reality?
Discuss how you could modify the ODE model to describe other circumstances, or to add other effects to it.