This site is devoted to mathematics and its applications. Created and run by Peter Saveliev.

# Differential equations

## Contents

- 1 Discrete models: how to set up differential equations
- 2 Solution sets of ODEs
- 3 Change of variables in ODEs
- 4 Separation of variables in ODEs
- 5 The method of integrating factors
- 6 Euler's method: back to discrete
- 7 Qualitative analysis of ODEs
- 8 How large is the difference between discrete and continuous?
- 9 Linearization of ODEs
- 10 Motion under forces: ODEs of second order

## Discrete models: how to set up differential equations

Below, we produce discrete models and then differential equations from verbal descriptions.

Recall first that an *augmented partition* of an interval $[a,b]$ is a sequence of $n$ *primary nodes*, $t_{i}$, alternating with the *secondary nodes*, $c_i$, of the partition:
$$ a=t_{0}\le c_1\le t_{1}\le c_2\le t_{2}\le ... \le t_{n-1}\le c_n\le t_{n}=b.$$
The *increments* of $t$ are:
$$\Delta t_i = t_i-t_{i-1},\ i=1,2,...,n.$$

Typically, we have:

- the intervals of the partition of equal length, $\Delta t_i=\Delta t=h$, and
- the secondary nodes placed in the middle of each interval, $c_i=\frac{1}{2}(t_{i-1}+t_i)$.

Furthermore, the *difference* of a function $y$ defined at the primary nodes is a function defined at the secondary nodes of the partition:
$$\Delta y\, (c_{k})=y(t_{k})-y(t_{k-1}),\ k=1,2,...,n.$$
The indices are often omitted, $\Delta y\, (c_{k}) =\Delta y\, (c)$.

**Example (uniform motion).** Description: “the velocity is constant”. This means that the increment of the location $\Delta y$ is proportional to the increment of time $\Delta t$:
$$\Delta y=k\cdot \Delta t,$$
for some constant $k$. Suppose we also have an initial condition:
$$y(t_0)=y_0.$$
Then our equation gives us a sequence via this *recursive formula*:
$$y_{n+1}=y_n+k\cdot \Delta t.$$
We can also think of this sequence as a function defined at the nodes of the partition:
$$y(t_{n+1})=y(t_n)+k\cdot \Delta t,$$
where
$$t_{n+1}=t_n+\Delta t.$$
It is then a $0$-form. For example, here are a few solutions of the equation:
$$\Delta y=2 \Delta t,$$

They are just arithmetic progressions (with the same increment $2\Delta t$). Generally, they don't have to be -- when the partition is uneven:

But the points still lie on the same straight lines! $\square$

**Example (location from velocity).** The formula has been used in Parts I, II, and III to model uniform motion and acquire location from this velocity. More generally, we took any function $z=f(t)$ defined at the secondary nodes of the partition and think of $f(c_i)$ as the value of the velocity during the time interval $[t_{i-1},t_i]$. It is then a $1$-form. We set up a recursive equation:
$$y(t_{n+1})=y(t_n)+f(c_n)\cdot \Delta t,$$
where
$$t_{n+1}=t_n+\Delta t.$$
What kind of function is $f$? It is function of one variable $z=f(t)$. However, the graph of a solution $y_n=y(t_n)$ is to be plotted on the $ty$-plane. It is, therefore, beneficial to think of $f$ as a function of two variables $z=f(t,y)$ that just happens to have the same value for each $t$ regardless of $y$:

The visualization that shows the values of $f$ in terms of color -- blue for negative and red for positive -- is the most convenient of these. Indeed, as $y$ is plotted point by point, we know that we go up when the area is red and down when it's blue:

In the meantime, the recursive formula is seen as a *differential* equation. There are two equivalent forms: one in terms of *differences*,
$$\Delta y=f(c)\cdot \Delta t,$$
and another in terms of *difference quotients*,
$$\frac{\Delta y}{\Delta t}=f.$$
Furthermore, if we think of $y$ as a sampled continuously differentiable function (and $f$ as a sampled continuous function), we can take the limit of the latter equation, $\Delta t\to 0$. The result is the third equation, in terms of *derivatives*,
$$y'=f(t).$$
$\square$

**Example (population growth).** Description: “the rate of growth/decay of the population is proportional to the current population”. Then the increment of the population $\Delta y$ is proportional to $y$... and $\Delta t$:
$$\Delta y=k\cdot y\cdot \Delta t,$$
for some constant $k$.
The “right-hand” function is very simple:
$$f(y)=ky.$$
Once again, however, the graph of a solution $y_n=y(t_n)$ is plotted on the $ty$-plane and, therefore, it is better to think of $f$ as a function of two variables $z=f(t,y)$ that just happens to have the same value for each $y$ regardless of $t$:

Suppose we also have an *initial condition*:
$$y(t_0)=y_0.$$
Then the equation gives us a discrete model via this recursive formula:
$$y(t_{n+1})=y(t_n)+ky(t_n)\cdot \Delta t.$$
The formula has been used in Chapter 10 to study exponential models and now we will apply it to models with the rate of change of a quantity depending on the quantity's current value:

We see how we progress in the upward direction which makes the values of $f$ higher and, as a result, the growth of $y$ faster. Our discrete differential equations are: $$\Delta y=ky\cdot \Delta t\ \text{ and }\ \frac{\Delta y}{\Delta t}=k y .$$ In the meantime, if we think of $y$ as a sampled differentiable function, the equation is converted to a new equation: $$\lim_{\Delta t\to 0}\frac{\Delta y}{\Delta t}=ky \ \Longrightarrow\ y'=ky.$$ $\square$

**Example (logistic growth).** Description: “the rate of growth of the population is proportional to the current population and to the remaining population of the total that can be sustained”. The increment of the population $\Delta y$ is proportional to $y$... and $T-y$, where $T$ is the total possible population:
$$\Delta y=k\cdot y\cdot (T-y)\cdot \Delta t.$$
This is our function ($k=5,\ T=1$):

We see low level of potential growth in the areas where the population is low or where it is too close to be exhausting the resources. With an initial condition: $$y(t_0)=y_0,$$ we produce a discrete model via this recursive formula: $$y(t_{n+1})=y(t_n)+ky(t_n)(T-y(t_n))\cdot \Delta t,\ k>0.$$ The model describes a restricted growth ($k=5,\ \Delta t=.05$):

We see the fastest growth when the population is about half of what can be sustained. All solutions exhibit this behavior:

Choosing $k$ large relative to $\Delta t$ may break the model ($k=25,\ \Delta t=.05$):

Then the discrete differential equations are: $$\Delta y=ky(T-y)\cdot \Delta t\ \text{ and }\ \frac{\Delta y}{\Delta t}=y(T-y) .$$ If $y$ is a sampled differentiable function, then the limit of the latter equation produces the following: $$y'=ky(T-y).$$ $\square$

**Example (Newton's Law of Cooling).** Description: “the rate of cooling is proportional to the difference of the current temperature and the room temperature”. Then the increment of the population $\Delta y$ is proportional to $r-y$, where $r$ is the room temperature, and $\Delta t$ ($k>0$):
$$\Delta y=k\cdot (r-y)\cdot \Delta t.$$
The situation is similar to the last example, but this time
$$f(t,y)=k\cdot (r-y).$$
It is plotted below ($k=10,\ \Delta t=.05$):

We can see that when the temperature is higher than the room temperature, it decreases and when the temperature is lower than the room temperature, it increases. This is why the object's temperature can't “overshoot” that of the room. The equation gives us a discrete model via this recursive formula: $$y(t_{n+1})=y(t_n)+k(r-y(t_n))\cdot \Delta t.$$

The formula was used in Chapter 10 to model Newton's Law of Cooling:

Choosing $k$ large relative to $\Delta t$ may break the model ($k=25,\ \Delta t=.05$):

Then the discrete differential equations are: $$\Delta y=k(r-y)\cdot \Delta t\ \text{ and }\ \frac{\Delta y}{\Delta t}=k(r-y) .$$ In the meantime, the derivative, if any, would satisfy the following: $$y'=k(r-y).$$ $\square$

**Example (flow).** Description: “the velocity of flow along a pipe (or a canal with the water that has the exact same velocity at all locations across it) is given for each location”. Trace a single particle of this stream.

The initial time $t_0$ and the initial location $p_0$ are placed in the first row of the spreadsheet. As we progress in time and space, new numbers are placed in the next row of our spreadsheet: $$t_n,\ v_n,\ y_n,\ n=1,2,3,...$$ The same recursive formulas are used:

- $t_{n+1}=t_n+\Delta t$,
- $y_{n+1}=y_n+v_{n+1}\cdot \Delta t$.

But where does the velocity come from? It may come from a function that gives the velocity at every location: $$v_n=f(y_n).$$ But the velocity of the flow may vary with time! We have then: $$v_n=f(t_n,y_n).$$

Once again, we think of $y_n$ as a function: $y_n=y(t_n)$. We substitute this, as well as $t=t_n$, into the recursive formula for $y_{n+1}$:
$$y(t+\Delta t)=y(t)+f(t,y(t))\cdot \Delta t,$$
or
$$\Delta y=f(t,y,t)\cdot \Delta t.$$
Furthermore,
$$\frac{\Delta y}{\Delta t}=f(t,y).$$
These equations are the most general form of a discrete differential equation. We will use the former for simulation but will not attempt to *solve* the latter. In the rest of this chapter, we will try to solve some of the ODEs for derivatives that emerge from these relations: $y'=f(t,y)$.

This is an example of a possible choice of $f$ and some solutions of the corresponding ODE:

$\square$

**Example (data flow).** The quantities used to compute the next location (or a state) -- current time and velocity -- doesn't have to come from a formula! This is nothing but data and it may come as *strings of numbers* (or it can be fed into the computer continuously). For example, the initial time $t_0$ and the initial location $y_0$ are placed in the first row of the spreadsheet and, as we progress in time and space, new numbers are placed in the next row of our spreadsheet:
$$t_n,\ v_n,\ y_n,\ n=1,2,3,...$$
The same recursive formula is used to find the next location:
$$y_{n+1}=y_n+v_{n+1}\cdot \Delta t.$$
The result is a growing table of values:
$$\begin{array}{c|c|c|c|c|c}
&\text{iteration } n&\text{time }t_n&\text{velocity }v_n&&\text{location }y_n\\
\hline
\text{initial:}&0&3.5&--&&22\\
&1&3.6&33&&25.3\\
&...&...&...&&...\\
&1000&103.5&4&&336\\
&...&...&...&&...\\
\end{array}$$
All but the last column may come as a data file. $\square$

## Solution sets of ODEs

Now, suppose the velocity comes from an explicit formula as a function of the location $z=f(y)$ defined on an interval $J$, is there an explicit formula for the location as a function of time $y=y(t)$ defined on an interval $I$? The last equation is satisfied for *every* $\Delta t>0$ small enough. Taking the limit over $\Delta t\to 0$ gives us the following:
$$y'(t)=f(y(t)),$$
provided $y=y(t)$ is differentiable at $t$ and $z=f(y)$ is continuous at $y(t)$. Such an equation may be *possible to solve*.

Let's review the differential equations we have seen. We've seen many kinds: $$\begin{array}{lllll} &\text{equation}&&&\text{ a solution}\\ \hline 1.&f'(x)&=x^2&\longrightarrow& f(x)=x^3/3&\\ 2.&f'(x)&=f(x)&\longrightarrow& f(x)=e^x&\\ 3.&f' '(x)&=2&\longrightarrow& f(x)=x^2&\\ 4.&f' '(x)&=-f(x)&\longrightarrow& f(x)=\sin x& \end{array}$$ When $x$ represents time and $y=f(x)$ represents location, the equations have simple interpretations:

- 1. the velocity is known at every moment of time: missile;
- 2. the velocity is known at every location: liquid flow;
- 3. the acceleration is known at every moment of time: free fall;
- 4. the acceleration is known at every location: oscillating spring.

This interpretation justifies using $t$ for the name of the independent variable. In the field of differential equations, it is also common *not to name the functions* anymore but use instead the names of the variables under this more compact **notation**:
$$\begin{array}{lllll}
1.&y'&=t^2&\longrightarrow& y=t^3/3 &\text{ is a solution;}\\
2.&y'&=y&\longrightarrow& y=e^t &\text{ is a solution;}\\
3.&y' '&=2&\longrightarrow& y=t^2 &\text{ is a solution;}\\
4.&y' '&=-y&\longrightarrow& y=\sin t &\text{ is a solution.}
\end{array}$$

The last two equations will be studied later in this chapter. The other two are generalized as follows.

First, we are after functions $y=y(t)$ that satisfy the equation:
$$y'(t)=g(t),$$
for some function $g$. It is *location-independent*. Here, for every $t$, we know the slope of the tangent line at $(t,y(t))$. This slope is the same along any given *vertical* line:

When the function $g(x)=k$ happens to be constant, its solution is the family of vertically shifted parabolas:

They are represented by: $$y=-\frac{1}{2}kt^2+C.$$

Second, we are after functions $y=y(x)$ that satisfy the equation:
$$y'(t)=h(y(t)),$$
for some $h$. It is *time-independent*. Here, for every $y$, we know the slope of the tangent line at $(t,y)$. This slope is the same along any given *horizontal* line:

When the function $h(x)=x$ happens to be the identity, its solution is the family of vertically stretched curves:

They are represented by: $$y=Ce^t.$$

**Definition.** Suppose $f$ is some function of two variables. Then a *discrete ordinary differential equation of first order*, or simply an ODE, is defined to be the following:
$$\frac{\Delta y}{\Delta t}=f(t,y);$$
and an *ordinary differential equation of first order*, or simply an ODE, is defined to be the following:
$$\frac{dy}{dt}=f(t,y)\text{ or } y'=f(t,y).$$

**Definition.** Suppose $h$ is some function of three variables. Then a *discrete ordinary differential equation of second order* is defined to be the following:
$$\frac{\Delta^2 y}{\Delta t^2}=h\left(t,y,\frac{\Delta y}{\Delta t}\right),$$
and an *ordinary differential equation of second order* is defined to be the following:
$$\frac{d^2y}{dt^2}=h\left(t,y,\frac{dy}{dt}\right)\text{ or } y' '=h(t,y,y').$$

Here, “ordinary” refers to the fact that there is only one independent variable and only one derivative. The “order” refers to the order of the derivative. We will refer to $f$ and $h$ as the “right-hand side”.

Generally, the solution set of an ODE has neither of the above patterns:

Consider the two examples of ODEs in the last section: $$\begin{array}{lll} y'=-gt+v_0& \Longrightarrow& y(t)=-gt^2/2+v_0t+C;\\ y'=y& \Longrightarrow& y(t)=Ce^t.\\ \end{array}$$ There is one solution for each $C$, as we know. Each of them has domain $(-\infty, \infty)$.

Unlike these two ODEs, some others have solutions that cannot be defined to the whole real line.

**Example.** Some of them simply run away. For example, here is a simple example of such an ODE:
$$y'=1/x.$$
The right-hand side function,
$$f(x,y)=1/x,$$
is undefined at $x=0$. For a function of two variables, its domain consists now of two open half-planes:
$$x<0 \text{ and } x>0.$$
Therefore no solution to our ODE can cross the line $x=0$.

We have, in fact, two families of solutions: $$\begin{array}{lll} y=\ln (-x) +C&\text{ for } x>0;\\ y=\ln (x) +K&\text{ for } x<0. \end{array}$$ They exist separately from each other and should not be combined arbitrarily. $\square$

This is the reason why we take the following point of view on solutions.

**Definition.** A *solution* of ODE $y'=f(t,y)$ is a function $y$ differentiable on an open interval $I$ such that for every $t$ in $I$ we have:
$$y'(t)=f(t,y(t)).$$

Note that we require the domain to be an *open* interval so that the differentiability is properly defined.

The very first observation we should make is that, according to the definition, there is a separate solution for each interval inside $I$, i.e., if $y$ is a solution on interval $I$ then so is its *restriction* $z=y\big|_J$ to any open interval $J$ in $I$.

This new function is defined very simply: $$z(t)=y(t) \text{ for each } t \text{ in }J.$$

For the above example, each segment of a curve seen below is a solution:

Warning: interval $I$ may be infinite: $$I=(-\infty, +\infty),\ (-\infty,b),\ (a,+\infty).$$

Recall that solving an ODE means to face a “field of slopes” and then search for curves with these tangents:

Solving problems about the fall of a ball thrown (Chapter 7) from a particular location at a particular time has a related concept in the theory of ODEs.

**Definition.** For a given ODE $y'=f(t,y)$ and a given pair of values $(t_0,y_0)$, the *initial value problem*, or an IVP, is
$$y'=f(t,y),\quad y(t_0)=y_0;$$
and its *solution* is a solution of the ODE that satisfies the *initial condition*, $y(t_0)=y_0$.

In the example in the last section, the initial value problem is solvable for each $(t_0,y_0)$ in the domain. That is why the curves fill the plane. In the last example, they fill all plane except the line $t=0$.

**Definition.** We say that an ODE satisfies the *existence* property at a point $(t_0,y_0)$ when the IVP:
$$y'=f(t,y),\quad y(t_0)=y_0,$$
has a solution.

It doesn't matter how small is the domain of this solution.

If your model of a real-life process doesn't satisfy this property, it may be inadequate. It is as if the process starts but never continues.

From what we know about antiderivatives, we conclude the following.

**Theorem.** An ODE the right-hand side of which is a function independent of $y$ and integrable with respect to $t$ on an open interval $I$:
$$y'=f(t),$$
satisfies the existence property for every point $(t_0,y_0)$ with $t_0$ in $I$.

**Proof.** Indeed, the solution is just the appropriate antiderivative of $f$. $\blacksquare$

**Definition.** We say that an ODE satisfies the *uniqueness* property at a point $(t_0,y_0)$ if every pair of solutions, $y_1,\ y_2$, of the IVP:
$$y'=f(t,y),\quad y(t_0)=y_0,$$
are equal,
$$y_1(t)=y_2(t),$$
for every $t$ in some open interval that contains $t_0$.

In other words, their restrictions to some interval $J$ are equal: $$y_1\big|_J=y_2\big|_J.$$

If your model of a real-life process doesn't satisfy this property, it may be inadequate. It's as if you have all the data but can't predict even the nearest future.

From what we know about antiderivatives, we conclude the following.

**Theorem.** An ODE the right-hand side of which is a function independent of $y$ and integrable with respect to $t$ on an open interval $I$:
$$y'=f(t),$$
satisfies the uniqueness property for every point $(t_0,y_0)$ with $t_0$ in $I$.

**Proof.** Indeed, the solution is always one of the antiderivatives of $f$. We also know that these antiderivatives always differ by a constant and, therefore, their graphs never intersect. $\blacksquare$

We will see ODEs without the uniqueness property later.

**Definition.** We say that an ODE satisfies the *continuity of IVP* property at $(t_0,y_0)$ if for each solution $y$ of the IVP above there is $t_1 >t_0$ such that for each sequence of solutions, $\{ y_n \}$, of the ODE we have, as $n\to \infty$,
$$y_n(t_0)\to y(t_0) \ \Longrightarrow\ y_n(t_1 )\to y(t_1 ).$$

**Exercise.** What about convergence of *functions* $y_n\to y$? What about uniform convergence?

**Exercise.** Prove that the function $Q$ of the $y$-axis to itself defined by $Q(z)=y(t_1)$, where $y(t_0)=z$, is continuous. Hint:

If your model of a real-life process doesn't satisfy this property, it may be inadequate. It's as if any error in the initial data -- no matter how small -- may cause a significant error in your prediction.

From what we know about antiderivatives, we conclude the following.

**Theorem.** An ODE the right-hand side of which is a function independent of $y$ and integrable with respect to $t$ on an open interval $I$:
$$y'=f(t),$$
satisfies the continuity of IVP property for every point $(t_0,y_0)$ with $t_0$ in $I$.

**Proof.** Indeed, the solution is always one of the antiderivatives of $f$. We also know that these antiderivatives always differ by a constant and, therefore, converge to each other when this constant goes to zero. $\blacksquare$

Next, ODEs produce families of curves as the sets of their solutions... and vice versa: if a family of curves is given by an equation with a single parameter, then differentiating the equation will produce an ODE.

**Example.** The family of vertically shifted graphs,
$$y=g(t)+C,$$
creates an ODE the right-hand side of which is a function independent of $y$:
$$y'=g'(t).$$
The family of stretched exponential graphs,
$$y=Ce^t,$$
creates an ODE the right-hand side of which is a function independent of $t$:
$$y'=Ce^t\ \Longrightarrow\ y'=y.$$
Either of the two equations satisfy the three basic properties. $\square$

**Example (non-existence).** What about these concentric circles?

They are given by $$x^2+y^2=C>0.$$ We differentiate (implicitly) with respect to $x$: $$2x+2yy'=0,$$ or, $$y'=-\frac{x}{y}.$$ The domain of the right-hand side is $y\ne 0$. Therefore, the existence property cannot be satisfied on the $x$-axis. Furthermore, even if this function was defined on the $x$-axis, the existence property would still break down:

As you can see, we cannot extend the solution to the left of this point. Similarly, we cannot extend the solution to the right of the point if it is on the other side of $0$. $\square$

**Example (non-existence).** These hyperbolas are given by these equations:
$$xy=C.$$

(In the case when $C=0$, we have the two axes.) Then, we have an ODE: $$y+xy'=0,$$ or, $$y'=\frac{y}{x}.$$ The domain of the right-hand side is $x\ne 0$. Therefore, the existence property cannot be satisfied on the $y$-axis. Furthermore, even if this function was defined on the $y$-axis, the existence property would still break down:

We cannot have a solution passing through the $y$-axis. $\square$

It appears that the presence or the absence of the existence property depends on the continuity condition of the right-hand side function of the ODE. The proof of the following important theorem lies outside the scope of this book.

**Theorem (Existence).** Suppose $(t_0,y_0)$ is a point on the $ty$-plane and suppose

- $J$ is an open interval that contains $t_0$, and
- $G$ is an open interval that contains $y_0$.

Suppose also that a function $z=f(t,y)$ of two variables is

- continuous with respect to $t$ on $J\times G$, and
- continuous with respect to $y$ on $J\times G$.

Then the ODE $y'=f(t,y)$ satisfies the existence property at $(t_0,y_0)$.

Warning: the domain of the solution doesn't have to be the whole $J$ because it may exit the rectangle through its top or its bottom:

**Exercise.** What if $G=(-\infty,\infty)$?

**Example (non-uniqueness).** Consider this ODE:
$$\frac{dy}{dx}=y^{2/3}.$$
Notice that there is one more solution passing through $(0,0)$ in addition to the trivial solution $y=0$:

It is $$y=\begin{cases} 0&\text{ for } x<0;\\ \frac{1}{27}x^3&\text{ for } x\ge 0. \end{cases}$$ What makes the right-hand side $y^{2/3}$ special in comparison to the previous examples? Its derivative is infinite at $y=0$. $\square$

It appears that the presence or the absence of the uniqueness property depends on the differentiability condition of the right-hand side function of the ODE. The proof of the following important theorem lies outside the scope of this book. The set-up is identical to that of the last theorem except: continuity for $y$ is replaced with differentiability.

**Theorem (Uniqueness).** Suppose $(t_0,y_0)$ is a point on the $ty$-plane and suppose

- $J$ is an open interval that contains $t_0$, and
- $G$ is an open interval that contains $y_0$.

Suppose also that a function $z=f(t,y)$ of two variables is

- continuous with respect to $t$ on $J\times G$, and
- differentiable with respect to $y$ on $J\times G$.

Then the ODE $y'=f(t,y)$ satisfies the uniqueness property at $(t_0,y_0)$.

Warning: this is *not* what non-uniqueness looks like:

It appears that the presence or the absence of the continuity of IVP property depends on the differentiability condition of the right-hand side function of the ODE. The proof of the following important theorem lies outside the scope of this book. The set-up is identical to that of the last theorem.

**Theorem (Continuity of IVP).** Suppose $(t_0,y_0)$ is a point on the $ty$-plane and suppose

- $J$ is an open interval that contains $t_0$, and
- $G$ is an open interval that contains $y_0$.

Suppose also that a function $z=f(t,y)$ of two variables is

- continuous with respect to $t$ on $J\times G$, and
- differentiable with respect to $y$ on $J\times G$.

Then the ODE $y'=f(t,y)$ satisfies the continuity of IVP property at $(t_0,y_0)$.

Next, once we have a solution $y=y(t)$ on interval $I$, any restriction $z=y\big|_J$ of function $y$ to any interval $J$ inside $I$ is also a solution. Conversely, once we have a solution $z$ on interval $J$, there may be an *extension* $y$ of $z$, i.e., $y=y(t)$ is a solution and $z=y\big|_J$. Sometimes it is impossible to further extend our solution.

**Definition.** A solution $y=y(t)$ is called a *maximal solution* if it doesn't have extensions.

As we have seen above, a maximal solution is often defined on $(-\infty,+\infty)$ but also may be defined on other open intervals.

Below, we will refer to maximal solutions as simply *solutions*.

## Change of variables in ODEs

Let's review a simple example of *integration by substitution*:
$$\int xe^{x^2}\, dx=?$$
Recognizing the presence of a composition:
$$y=e^{x^2},$$
we split it. We introduce an intermediate variable, $u$:
$$x\mapsto u\mapsto y,$$
with
$$u=x^2.$$
Then, the differential of $u$ is this differential form:
$$du=2xdx.$$
We substitute these two into the integral:
$$\int xe^{x^2}\, dx=\int xe^{x^2}\Bigg|_{x^2=u}\, dx\Bigg|_{dx=\frac{du}{2x}} =\int xe^u\frac{du}{2x}=\frac{1}{2}\int e^u\, dx.$$
At this point,

*change of variables*is complete, and- the new integral is simpler!

**Exercise.** Carry out the substitution $v=x^3$ in the above integral. Hint: yes, third power.

The idea of finding a new variable that *might* make the problem simpler is pervasive in mathematics. To complete the change, the derivatives are also to be modified, via the Chain Rule.

Let's recast the above construction in terms of differential equations.

The integral may be seen as the answer... to what question? An ODE. The ODE corresponding to the integral is:
$$y'=xe^{x^2}.$$
Let's use Leibniz notation:
$$\frac{dy}{dx}=xe^{x^2}.$$
We have *two variables* related via an unknown function: $y=y(x)$. We introduce a new variable:
$$u=x^2.$$
We have *three variables* related to each other via:

- an unknown function: $y=y(x)$ that is still to be found,
- the change of variables formula: $u=x^2$, and
- another unknown function: $y=y(u)$ to be found first.

The dependence diagram is below:
$$\begin{array}{ccc}
x&\to&u\\
&\searrow&\downarrow\\
&&y
\end{array}$$
The three relations between the three variables have their derivatives:
$$\frac{du}{dx}=2x,\ \frac{dy}{du}=?,\ \frac{dy}{dx}=xe^{x^2},$$
also related. These *related rates* are related via the Chain Rule:
$$\frac{du}{dx}\cdot \frac{dy}{du}=\frac{dy}{dx}.$$
Therefore, our ODE becomes after substitution:
$$\frac{dy}{dx}=xe^{x^2}\ \Longrightarrow\ \frac{du}{dx}\cdot \frac{dy}{du}=xe^u\ \Longrightarrow\ 2x\frac{dy}{du}= xe^u \ \Longrightarrow\ \frac{dy}{du}=\frac{1}{2}e^u.$$
Then,

*change of variables*is complete, and- the new ODE is simpler!

**Exercise.** Carry out the substitution $v=x^3$ in the above ODE.

We solve the new ODE easily:
$$y=\int \frac{1}{2}e^u\, du=\frac{1}{2}e^u+C.$$
Then the *back-substitution*, $u=x^2$, gives us the solution to the original ODE:
$$y=\frac{1}{2}e^{x^2}+C.$$

All of the above substitutions introduce a new *independent* variable. A new *dependent* variable may also simplify things.

The simplest such substitution is the *shift*:
$$y=z+a,$$
where $z$ is the new dependent variable. To analyze the ODE “locally”, we might want to concentrate on a single point $y=a$ (and its vicinity) at a time. The shift allows us to move that point to zero. For example, suppose we have a *linear ODE*:
$$y'=my+b,\ m\ne 0.$$
Let's choose
$$z=y-a, \text{ or } y=z+a, \text{ with } a=-\frac{b}{m}.$$
Then, first,
$$y'=(z+a)'=z',$$
and, second,
$$my+b=m(z+a)+b=mz+ma+b=mz+m\left(-\frac{b}{m}\right)+b=mz.$$
We have a new ODE with respect to $x$:
$$z'=mz.$$
It is simpler!

## Separation of variables in ODEs

Differential forms allow us sometimes to separate $x$, and $dx$, from $y$, and $dy$.

**Example.** We re-write the ODE as an equation of differential forms and then integrate both sides:
$$\frac{dy}{dx}=x^2\ \Longrightarrow\ dy=x^2\, dx\ \Longrightarrow\ \int\, dy=\int x^2\, dx\ \Longrightarrow\ y+C=x^3/3+K.$$
The two indefinite integration produce two constants. Since both are arbitrary, just one is sufficient:
$$y=x^3/3+K.$$
The ODE is solved. $\square$

**Example.** Following the same procedure, we integrate both sides of an equation of differential forms keeping only one constant:
$$\frac{dy}{dx}=y\ \Longrightarrow\ \frac{dy}{y}=dx\ \Longrightarrow\ \int\, \frac{dy}{y}=\int \, dx\ \Longrightarrow\ \ln y=x+C \text{ for }y>0.$$
This is a part of the solution set of the ODE given implicitly. Similarly, we obtain:
$$\ln (-y)=x+K \text{ for }y<0.$$
What's left is the case of $y=0$. That's a whole (constant) solution. We verify this fact by substitution. The picture this result produces is familiar:

$\square$

In general, the method applies whenever our right-hand side function splits into a factor of two functions with one depending only on $x$ and the other only on $y$:
$$\frac{dy}{dx}=p(x)q(y)\ \Longrightarrow\ \frac{dy}{q(y)}=p(x)\, dx.$$
Such ODEs are called *separable*.

**Theorem.** Every solution of a separable ODE $y'=p(x)q(y)$ satisfies the equation:
$$\int\frac{dy}{q(y)}=\int p(x)\, dx.$$

The question then remains whether the two integrals can be evaluated in an elementary fashion.

**Example.** Consider this separable ODE:
$$\frac{dy}{dx}=y\sin x\ \Longrightarrow\ \frac{dy}{y}=\sin x\, dx\ \Longrightarrow\ \int\, \frac{dy}{y}=\int \sin x\, dx\ \Longrightarrow\ \ln y=-\cos x+C \text{ for }y>0.$$

The ODE is solved implicitly. Even though some of the solutions appear to touch the $x$-axis, the uniqueness is satisfied. $\square$

**Exercise.** Prove the last statement.

Let's apply the method to an ODE that we will see later.

**Theorem.** The ODE
$$y'=a(x)y,$$
with an integrable function $a$, satisfies the existence and uniqueness. The solutions of the ODE are given by:
$$y=Ce^{A(x)},$$
where $A$ is any antiderivative of $a$:
$$A(x)=\int a(x)\, dx.$$

Thus the solution set splits into three parts: $$(1)\ y=Ke^{A(x)}\text{ with }K<0,\ (2)\ y=0,\ (3)\ y=Ce^{A(x)}\text{ with }C>0.$$

**Exercise.** Prove the theorem.

## The method of integrating factors

**Example.** Let's take a careful look at the left-hand side of the following (non-separable) equation:
$$y'\cdot \sin x+y\cdot \cos x =xe^{x^2}.$$
We see two pairs of functions present:

- $y$ and its derivative $y'$, and
- $\sin x$ and its derivative $\cos x$.

They are “cross-multiplied”. The expression is familiar; it's the outcome of the Product Rule: $$\left( y \cdot \sin x \right)'=y'\cdot \sin x+y\cdot \cos x .$$ Therefore, our ODE turns into: $$\left( y \cdot \sin x \right)'=xe^{x^2}.$$ We integrate: $$\int \left( y \cdot \sin x \right)' \, dx=\int xe^{x^2}\, dx,$$ and apply the Fundamental Theorem of Calculus: $$y \cdot \sin x=\frac{1}{2}e^{x^2}+C.$$ $\square$

We were lucky.

**Example.** Consider the ODE:
$$y'+y=x.$$
The left-hand side is *not* the outcome of the Product Rule... Can we make it so? After all, the pair $y$ and $y'$ is already there. Unfortunately, what's left doesn't work:
$$y'\cdot 1+y\cdot 1=x.$$
We can't just replace the $1$'s with some functions $F$ and its derivative $F'$, can we? We can if this is the *same* function:
$$\left( e^x \right)'=e^x.$$
We multiply both sides of the equation by this factor $F(x)=e^x$:
$$y'\cdot e^x+y\cdot e^x=xe^x,$$
or
$$y'\cdot e^x+y\cdot \left( e^x \right)'=xe^x.$$
By the Product Rule we have instead:
$$\left( y \cdot e^x \right)'=xe^x.$$
After integration, we have
$$y \cdot e^x =-e^x+xe^x+C,$$
or even simpler:
$$y =-1+x+Ce^{-x}.$$
$\square$

Then, it is promising for this approach to have both $y$ and $y'$, but not $y^2$ or $\sin y'$. In other words, an equation should be *linear* with respect to $y$ and $y'$.

**Example.** We make a small adjustment to the ODE:
$$y'+2y=x.$$
The left-hand side is not the outcome of the Product Rule... even if we multiply by $e^x$:
$$y'\cdot e^x+2y\cdot e^x=xe^x.$$
What choice of a factor $F$ would work:
$$y'\cdot F(x)+2y\cdot F(x)=xF(x)?$$
We'd need this:
$$F'=2F.$$
But that's just another ODE! And a familiar one too:
$$F(x)=e^{2x}.$$
Thus we have:
$$y'\cdot e^{2x}+2y\cdot e^{2x}=xe^{2x},$$
or
$$y'\cdot e^{2x}+y\cdot \left( e^{2x} \right)'=xe^{2x}.$$
Then,
$$\left( y \cdot e^{2x} \right)'=xe^{2x}.$$
After integration, we have
$$y \cdot e^{2x} =-e^{2x}/2+xe^{2x}+C.$$
$\square$

In general, when the linear equation
$$y'+a(x)y=b(x)$$
is multiplied by some factor $F$, we have
$$y'F(x)+a(x)yF(x)=b(x)F(x).$$
Then, the approach via the Product Rule applies when that equation is identical to this equation:
$$y'F(x)+yF'(x)=b(x)F(x),$$
i.e., when $F$ satisfies:
$$F'=a(x)F.$$
This ODE always has a solution according to the last theorem:
$$F=e^{\int a(x)\, dx},$$
as long as function $a$ is integrable. Function $F$ defined this way is called the *integrating factor* of the linear equation.

**Theorem.** Every solution of the ODE
$$y'+a(x)y=b(x)$$
is given by:
$$y=e^{-A(x)}\int b(x)e^{A(x)}\, dx,$$
where $A$ is any antiderivative of $a$:
$$A(x)=\int a(x)\, dx.$$

## Euler's method: back to discrete

Unfortunately, an ODE typically has no algebraic solution!

We would like to approximate solutions of a general IVP: $$y'=f(t,y),\ y(t_0)=y_0.$$ The IVP tells us:

- where we are (the initial condition), and
- the direction we are going (the ODE).

It is then easy to produce a solution if we assume that the change is incremental: as we arrive to a new location, we, once again, know where we are and where we are going.

The method presented here follows the idea from Chapter 10: the unknown solution is replaced with its *best linear approximation*.

**Example.** Let's consider again these concentric circles:

They are the solutions of the ODE: $$y'=f(x,y)=-\frac{x}{y}.$$ We will be solving numerous initial value problems while staying away from the $x$-axis where there are no solutions.

We choose the increment of $x$: $$\Delta x=1.$$

We start with this initial condition: $$x_0=0,\quad y_0=1.$$ We substitute these two numbers into the equation: $$y'=-\frac{0}{1}=0.$$ This is the slope of the line we will follow. How far? As far as the increment of $x$ allows. Therefore, the increment of $y$ is $$\Delta y=0\cdot \Delta x=0\times 1 =0.$$ Our next location on the $xy$-plane is then: $$x_1=x_0+\Delta x=0+1=1,\quad y_1=y_0+\Delta y=1+0=1.$$

This computation gives us a *new initial condition*:
$$x_1=1,\quad y_1=1.$$
We again substitute these two numbers into the equation:
$$y'=-\frac{1}{1}=-1.$$
This is the slope of the line we will follow. Therefore, the increment of $y$ is
$$\Delta y=-1\cdot \Delta x=-1.$$
Our next location on the $xy$-plane is then:
$$x_2=x_1+\Delta x=1+1=2,\quad y_2=y_1+\Delta y=1+(-1)=0.$$
We have ended up on the $x$-axis and stop. These three points form an approximate solution.

We start all over from another initial condition: $$x_0=0,\quad y_0=2.$$ We substitute these two numbers into the equation: $$y'=-\frac{0}{2}=0$$ producing the slope we will follow. The increment of $y$ is $$\Delta y=0\cdot \Delta x=0\cdot 1 =0.$$ Our next location on the $xy$-plane is then: $$x_1=x_0+\Delta x=0+1=1,\ y_1=y_0+\Delta y=2+0=2.$$ A new initial condition appears: $$x_0=1,\quad y_0=2.$$ We again substitute these two numbers into the equation: $$y'=-\frac{1}{2}=-1/2,$$ producing the slope we will follow. The increment of $y$ is $$\Delta y=-1/2\cdot \Delta x=-1/2.$$ Our next location on the $xy$-plane is then: $$x_2=x_1+\Delta x=1+1=2,\quad y_2=y_1+\Delta y=2+(-1/2)=3/2.$$ One more IVP: $$x_2=2,\ y_2=3/2.$$ The increment of $y$ is $$\Delta y=-\frac{x}{y}\cdot \Delta x=-\frac{2}{3/2}\cdot 1 =-4/3.$$ Our next location on the $xy$-plane is then: $$x_3=x_2+\Delta x=2+1=3,\quad y_3=y_2+\Delta y=3/2-4/3=1/6.$$ We would have to pass the $x$-axis with next step and we stop now. These four points, in addition to the previous three, form a very crude approximation of our circular solutions:

$\square$

From the point of view of motion, this is the method's interpretation:

- at our current location and current time, we examine the ODE to find the velocity and then move accordingly to the next location.

**Definition.** The *Euler solution* with increment $\Delta t>0$ of the IVP
$$y'=f(t,y),\ y(t_0)=y_0,$$
is the sequence $\{y_n\}$ of real numbers given by:
$$t_{n+1}=t_n+\Delta t,\quad y_{n+1}=y_n+f(t_n,y_n)\cdot \Delta t.$$

The most important fact to know about Euler's method is that if we derived our ODE from a discrete ODE (via $\Delta t\to 0$), Euler's method will bring us right back to it: $$\newcommand{\la}[1]{\!\!\!\!\!\xleftarrow{\quad#1\quad}\!\!\!\!\!} \begin{array}{cccccc} &&\text{ODE}\\ &\nearrow&&\searrow\\ \text{discrete ODE}&&\la{same!}&&\text{Euler's method} \end{array}$$

**Example.** Let's now carry out this procedure with a spreadsheet. The formula for $y_n$ is:
$$\texttt{=R[-1]C+R[-1]C[1]*R2C1}$$
The results are less crude:

However, the approximations appear to behave erratically when they get close to the $x$-axis. The reason is the division by $y$ which may be very small. In fact the solutions are supposed to stop at the $x$-axis, but, by design of Euler's method, they can't.

Elsewhere, the Euler solutions are close to perfect:

$\square$

Warning: Euler's method might ignore the non-existence property of the ODE.

**Example.** Let's consider again these hyperbolas:
$$xy=C.$$

They are the solutions of the ODE: $$y'=\frac{y}{x}.$$ For larger values of $h$, the Euler solution might cross the $x$-axis demonstrating non-uniqueness:

Several better looking Euler solutions are shown below:

Even then, is this asymptotic convergence toward the $x$-axis or do they merge? $\square$

Warning: Euler's method might introduce non-uniqueness to an ODE.

## Qualitative analysis of ODEs

Unfortunately, Euler's method depends on the value of $h$ and, from simple experimentation, one cannot know whether $h$ is small enough. It's just an approximation! The method can also fail near the boundary of the domain.

With these limitation of the *quantitative* methods, *qualitative* analysis provides fully accurate if broad descriptions of the solutions. The method amounts to gathering information about the solutions without solving the ODE -- either analytically or numerically.

The main tools come from Chapter 9. First, for a differentiable function $y$ defined on an open interval we have: $$\begin{array}{lll} y'>0\ \Longrightarrow\ y\nearrow;\\ y'<0\ \Longrightarrow\ y\searrow. \end{array}$$

**Theorem.** If $y$ is a solution of the ODE $y'=f(t,y)$, then on each open set in the $ty$-plane, we have:
$$\begin{array}{lll}
f(t,y)>0\ \Longrightarrow\ y \text{ is increasing};\\
f(t,y)<0\ \Longrightarrow\ y \text{ is decreasing}.
\end{array}$$

The second tool is also elementary.

**Theorem.** If the right-hand side of an ODE $y'=f(y)$ is independent of $t$ and $f(y_0)=0$, then $y(t)=y_0$ is a (constant) solution.

**Example.** Consider:
$$y'=-\tan (y).$$
First, the right-hand side, the tangent, is undefined at $y=\pi/2+k\pi$ for all integer $k$. Therefore, these horizontal lines on the $ty$-plane will cut the domain into strips so that every solution will stay inside one of them. Furthermore, the sign of the tangent changes at those values of $y$ as well as at $y=k\pi$, in the middle of each strip. For example, this is what we conclude about the strip $[-\pi/2,\pi/2]$:

- for $-\pi/2<y<0$, we have $y'=-\tan y>0$ and, therefore, $y\nearrow$;
- for $y=0$, we have $y'=-\tan y=0$ and, therefore, $y$ is a constant solution;
- for $0<y<\pi/2$, we have $y'=-\tan y<0$ and, therefore, $y\searrow$.

The results of the analysis are summarized below with the monotonicity of the solution matched with the sign of the tangent:

Furthermore, every solution $y$ is decreasing (or increasing) throughout its domain. If we can demonstrate existence, every solution is approaching the middle line of the strip, $y=0$, which is also a solution. We can then conclude that the latter is its horizontal asymptote, $y\to 0$ as $t\to\infty$, if we can demonstrate uniqueness. The conclusions are confirmed with Euler's method:

We notice now that the solutions are running away from the edge of the domain and toward the constant solution. We also notice that the solutions are repeated under a horizontal shift. $\square$

The last observation has important implications.

**Theorem.** If the right-hand side of an ODE $y'=f(y)$ is independent of $t$ and if $y=y(t)$ is a solution of the ODE then so is $y=y(t+s)$ for any real $s$.

**Exercise.** Prove the theorem.

**Exercise.** In the last example, is $y=\pi/2$ also a horizontal asymptote? Confirm your conclusion with Euler's method. What is $y'$ near this line?

Notice that both above and below Euler's method is run for various values of $t_0=0,1,2,3...$ in order to avoid empty spaces created by the asymptotic behavior of the solutions.

**Example.** Consider next:
$$y'=\cos y \cdot \sqrt{|y|}.$$
We start with the “sign analysis” of the right-hand side, $f(y)=\cos y \cdot \sqrt{|y|}$, which is fully determined by $\cos y$ unless $y=0$:
$$\begin{array}{r|c|l|l}
y&f(y)&y=y(t)\\
\hline
3\pi/2&0&y(t)=3\pi/2, \text{ constant}&\to\to\to\to\\
&-&y\searrow&\searrow\searrow\searrow\searrow\\
\pi/2&0&y(t)=\pi/2, \text{ constant}&\to\to\to\to\\
&+&y\nearrow&\nearrow\nearrow\nearrow\nearrow\\
0&0&y(t)=0, \text{ constant}&\to\to\to\to\\
&+&y\nearrow&\nearrow\nearrow\nearrow\nearrow\\
-\pi/2&0&y(t)=-\pi/2, \text{ constant}&\to\to\to\to\\
&-&y\searrow&\searrow\searrow\searrow\searrow\\
-3\pi/2&0&y(t)=-3\pi/2, \text{ constant}&\to\to\to\to\\
\end{array}$$
The results are confirmed with Euler's method:

The plotted Euler solutions cross the path of $y=0$ indicating non-uniqueness. Considering non-differentiability of $f(y)$ at $y=0$, this is a possibility. $\square$

**Example.** The next one is dependent on $t$:
$$y'=\frac{1}{t-y}.$$
The right-hand side is undefined whenever $t-y=0$. Therefore, no solution crosses the diagonal line $y=t$. Furthermore, the solutions above it are decreasing and the ones below it are increasing. Check:

$\square$

**Example.** The next one is defined on the whole plane:

$$y'=\cos (t+y).$$

The sign is changing whenever $t+y=\pi/2+k\pi$ for all integer $k$. Thus the solution change their monotonicity whenever they reach one of the lines $$y=-t+\pi/2+k\pi$$ of slope $-1$. One of these lines however is special and we wouldn't think of it without doing Euler's method first:

Indeed, $y=\pi-t$ is a solution! The rest of them are waves moving diagonally. $\square$

**Example.** The next one is discontinuous:
$$y'=[t+y].$$
Since the function changes abruptly, so does the derivative (the slope) of each solution. For Euler's method we use the $\texttt{FLOOR}$ function:

The existence is visible but to confirm we'd need a stronger version of our existence theorem. $\square$

Suppose the ODE is *time-independent*,
$$y'=f(y).$$
Then it is thought of a liquid flow: in a canal or a pipe.

Especially in the latter case, the right-hand side is recognized as a *one-dimensional vector field*.

With such an equation, the qualitative analysis is much simpler. In fact, one of the ODEs above, $f(y)=\cos y \cdot \sqrt{|y|},$ exhibits all possible patterns of *local* behavior.

We concentrate on what is going on in the vicinity of a given location $y=a$.

The first main possibility is $f(a)\ne 0$. Then, from the continuity of $f$, we conclude that $f(y)>0$ or $f(y)<0$ in some interval $I$ that contains $a$. Then, the solutions located within the band $(-\infty, +\infty)\times I$ are either all increasing or all decreasing respectively:

The behavior is “generic”.

More interesting behaviors are seen around a zero of $f$:

- $f(a)=0\ \Longrightarrow\ y=a $ is a stationary solution (an equilibrium).

Then the pattern in the vicinity of the point, i.e., an open interval $I$, depends on whether this is a maximum of $f$, a minimum, or neither.

First, if $f$ changes its sign from positive to negative, we have a *stable equilibrium*:

Indeed, we have all solutions within the band $(-\infty, +\infty)\times I$, under the uniqueness condition, asymptotically approach $y=a$: $$y\to a \text{ as } x\to +\infty.$$

Second, if $f$ changes its sign from negative to positive, we have an *unstable equilibrium*:

It means that there are no solutions within the band $(-\infty, +\infty)\times I$ that asymptotically approach $y=a$ as $x\to +\infty$. On the other hand, we have $$y\to a \text{ as } x\to -\infty.$$

Third, if $f$ does not change its sign at $y=a$, we have an *semi-stable equilibrium*:

It means that the solutions within one of the two sides of the line $y=a$ within the band $(-\infty, +\infty)\times I$ asymptotically approach this line and the ones on the other side do not.

The last possibility is $f(y)=0$ on the whole $I$. Then all the solutions in the band are stationary. This is a *degenerate equilibrium*.

This *classification of equilibria* is easy to understand in terms of the flow represented by the ODE:

- sink: flow in only (stable equilibrium);
- source: flow out only (unstable equilibrium);
- pass: flow, stop, flow (semi-stable equilibrium).

There are only these three main possibilities in the $1$-dimensional case. We will see in Chapter 24 a wide variety of behaviors around an equilibrium when the space of location is $2$-dimensional, a plane.

## How large is the difference between discrete and continuous?

**Example.** How far apart are the solutions of a pair of matching discrete and continuous ODEs? In other words, how far is the Euler solution from the actual solution of the ODE for the derivatives? For example, consider:
$$\frac{\Delta y}{\Delta t}=-\frac{t}{y} \text{ and } y'=-\frac{t}{y}.$$
The solutions of the latter are circles. So, if one such solution starts at $(0,y_0)$, i.e., $t=0$ and $y=y_0$, then the curve is supposed to end at the same location on the $t$-axis, i.e., $(y_0,0)$. This is not what we see when we lot a solution of the former equation:

Why does this happen? The solutions of the former decrease slower than those of the latter. Of course, the difference diminishes with the value of $h=\Delta t$. $\square$

We consider a solution of a discrete ODE and that of the corresponding continuous ODE with the same initial condition, $y(t_0)=y_0$. Their difference of the two is a function and the absolute value of this function evaluated at some $t>t_0$ is referred to as the *error*.

To estimate the error, we consider the *error bound of the best linear approximation* from Chapter 10. Suppose $y$ is twice differentiable at $t=t_0$ and $L(t)=y(t_0)+y'(t_0)(t-t_0)$ is its best linear approximation at $t_0$. Then, the error satisfies:
$$E(t)=|y(t)-L(t)| \le \tfrac{1}{2}K(t-t_0)^2,$$
where $K$ is a bound of the second derivative on the interval from $t_0$ to $t$:
$$|y' '(c)|\le K \text{ for all } c \text{ in this interval}.$$

Suppose now that $y$ is a solution of the IVP: $$y'=f(t,y),\ y(t_0)=y_0.$$ Provided its derivative is bounded as above, the error of a single step of the corresponding discrete ODE (i.e., Euler method) is bounded by: $$E(t_0+h)=|y(t_0 + h) - y_1| \le \tfrac{1}{2}Kh^2,$$ where $y_1=L_1(t_0+h)$ is the first step of the approximation that comes from $$L_1(t)=y_0+f(t_0,y_0)(t-t_0),$$ the linearization of $y$ at $t_0$.

We conclude that *the error of a single step is decreasing quadratically with $h$*.

Now, let's consider multiple steps. If $y' '$ is bounded by the same constant $K$, we can apply this formula to the second step of the approximation: $$|y(t_1 + h) - y_2| \le \tfrac{1}{2}Kh^2,$$ where $y_2=L_2(t_1+h)$ is the second step of the approximation that comes from $$L_2(t)=y_1+f(t_1,y_1)(t-t_1),$$ the linearization of $y$ at $t_1$.

Notice that the two error bounds apply to two *different* solutions $y$ of two different IVPs:

Therefore, the best bound for the error of the two steps is the sum of the two (in our case, identical) bounds: $$E(t_0+2h)=|y(t_0 + 2h) - y_2| \le 2\cdot \tfrac{1}{2}Kh^2.$$

And so on... The bound for the error of the $n$ steps is $n$ times the bound: $$E(t_0+nh)=|y(t_0 + nh) - y_n| \le n\cdot \tfrac{1}{2}Kh^2.$$

Suppose now that we extend the discrete solution with $n$ steps to $t_0+H$. Then $$n=\frac{H}{h}.$$ Therefore, $$E(t_0+nh)=|y(t_0 + nh) - y_n| \le \frac{H}{h}\cdot \tfrac{1}{2}Kh^2=\tfrac{1}{2}KHh.$$

We conclude that *the error is decreasing linearly with $h$*.

**Theorem (Error bound).** Suppose a differentiable function $z=f(t,y)$ satisfies
$${\partial f\over\partial t}(t, y) + {\partial f\over\partial y}(t, y) \, f(t, y)\le K,$$
for every $(t,y)$ in the rectangle $t_0\le t\le t_0+H,\ A\le y\le B$.
Suppose $y$ is a solution of the IVP
$$y'=f(t,y),\ y(t_0)=y_0,$$
with the domain $[t_0,t_0+H]$ and the range within $[A,B]$. Suppose $y_n$ is the $n$ steps of the solution of the discrete counterpart of this IVP with $h=H/n$. Then their difference satisfies:
$$|y(t_0 + nh) - y_n| \le \tfrac{1}{2}KHh.$$

**Proof.** The proof is above; we just use the *Chain Rule* to derive this:
$$y' '(t) = {\partial f\over\partial t}(t, y(t)) + {\partial f\over\partial y}(t, y(t)) \, f(t, y(t)).$$
$\blacksquare$

For each $h$, the error bound provides a tunnel, of constant width, around the discrete solution that contains the (unknown) solution of the continuous ODE, and vice versa:

**Exercise.** Prove the last statement. Hint: you need more than just $n$ points as in the theorem.

**Exercise.** Apply the theorem to the example of concentric circles in the beginning of the section.

**Example.** In our example of concentric circles, we also saw that the discrete solutions overshoot when they reach the $t$-axis. The reason is *concavity*. Suppose $y$ is a solution above the $t$-axis, then $y>0$ and $y'<0$. Then, we can find the sign of the second derivative of $y$:
$$y' ' =\frac{d}{dt}f(t,y)=\frac{d}{dt}\left(-\frac{t}{y} \right)=-\frac{1\cdot y-t\cdot y'}{y^2}=\frac{t\cdot y'-y}{y^2}<0.$$
It's negative! It's positive below the $t$-axis. Therefore, the solutions are concave down in the former and concave up in the latter case. Then, the linear approximations that we use for Euler's method overestimate the solutions in the former case and underestimate in the latter.

The benefit is that we are cutting the tunnel in half. $\square$

**Corollary (One-sided error bound).** Suppose a differentiable function $z=f(t,y)$ satisfies
$${\partial f\over\partial t}(t, y) + {\partial f\over\partial y}(t, y) \, f(t, y)\le K,$$
for every $(t,y)$ in the rectangle $t_0\le t\le t_0+H,\ A\le y\le B$.
Suppose $y$ is a solution of the IVP
$$y'=f(t,y),\ y(t_0)=y_0,$$
with the domain $[t_0,t_0+H]$ and the range within $[A,B]$. Suppose $h=H/n$ and $\{y_n\}$ is the solution of the corresponding discrete IVP. Then,

- when $f'>0$, we have for each $t$ in $[t_0,t_0+H]$:

$$y_n\le y(t_0 + nh) \le y_n+ \tfrac{1}{2}KHh,$$

- when $f'<0$, we have for each $t$ in $[t_0,t_0+H]$:

$$y_n-\tfrac{1}{2}KHh \le y(t_0 + nh) \le y_n.$$

**Exercise.** Apply the corollary to the example of concentric circles in the beginning of the section.

## Linearization of ODEs

Euler's method approximates solutions of an ODE by linearizing them, one location at a time. We can also approximate the solutions of the ODEs by linearizing, locally, its right-hand side.

The main idea is the following. As we are often unable to solve a differential equation we encounter, let's try to replace it with a simpler one that can be solved. How accurate do we need the approximation to be? At the very least, the new ODE should have the same *qualitative behavior* -- increasing, decreasing, and stationary solutions -- in the vicinity of the chosen point.

We will focus on the time-independent ODE with a continuous right hand-side, $$y'=f(y),$$ in the vicinity of a chosen point $y=a$. We will rely on the results stated previously.

**Theorem.** Suppose $f$ is continuous at $y=a$. Then, there is such an $\varepsilon>0$ that each solution $y$ of the ODE $y'=f(y)$ the graph of which lies within the band $(-\infty,+\infty)\times (a-\varepsilon, a+\varepsilon )$ satisfies:
$$\begin{array}{lll}
f(a)>0\ \Longrightarrow\ y \text{ is increasing};\\
f(a)<0\ \Longrightarrow\ y \text{ is decreasing}.
\end{array}$$

**Proof.** The continuity of $f$ implies that there is $\varepsilon>0$ small enough that for each $y$ within $(a-\varepsilon, a+\varepsilon )$, we have
$$\begin{array}{lll}
f(a)>0\ \Longrightarrow\ f(y)>0;\\
f(a)<0\ \Longrightarrow\ f(y)<0.
\end{array}$$
$\blacksquare$

Before we get to the linear, let's try the *constant approximation*. For a chosen point $y=a$,

- we replace $z=f(y)$ with $z=f(a)$.

In other words, the value $f(a)$ of $f$ at $a$ is used the right-hand side of the new ODE:
$$y'=f(a),$$
for all $y$ in some open interval $I$ that contains $a$. Its solutions are these *linear* functions of $t$:
$$y=f(a)(t+C),$$
for each real $C$.

Does their behavior match the original? There are two main cases:
$$f(a)\ne 0 \text{ and }f(a)= 0.$$
This is what we conclude about *all* solutions of the simplified ODE. In the former case, we have:
$$f(a)>0\ \Longrightarrow\ y\nearrow \text{ and } f(a)<0\ \Longrightarrow\ y\searrow ,$$
for all solutions in the band $(-\infty,+\infty)\times I$.

To compare to the solutions of the original ODE, we apply the last theorem. It's a match!

In the latter case, we have: $$f(a)=0\ \Longrightarrow\ y\text{ is constant}.$$

Then all the solutions in the band are stationary. But the original ODE might have a semi-stable equilibrium, a mismatch!

The answer is the *best linear approximation* of $f$. Suppose $f$ is differentiable at $a$.

The process of linearization starts with moving the point of interest, currently $y=a$, to $0$. Such a substitution is demonstrated earlier in the chapter. This is the new ODE: $$y'=f(y) \text{ with } f(0)=0.$$

The linearized ODE is
$$y'=f'(0)y,$$
for all $y$ in some open interval $I$ that contains $0$. Its solutions are these *exponential* functions of $t$:
$$y=Ce^{f'(0)t},$$
for each real $C$.

Does their behavior match the original? There are two main cases:
$$f'(0)\ne 0 \text{ and }f'(0)= 0.$$
This is what we conclude about *all* solutions of the linearized ODE. First we have:
$$f'(0)<0\ \Longrightarrow\ f\searrow \ \Longrightarrow\ f(y)>0\text{ for }y<0 \text{ and }f(y)<0\text{ for }y>0,$$
for all solutions in the band $(-\infty,+\infty)\times I$. Then, $f$ changes its sign from positive to negative and we have an unstable equilibrium:

Unlike the constant approximation, the linearization has produced a match!

Second, we have within the band: $$f'(0)>0\ \Longrightarrow\ f\nearrow \ \Longrightarrow\ f(y)<0\text{ for }y<0 \text{ and }f(y)>0\text{ for }y>0.$$ Then, $f$ changes its sign from negative to positive and we have a stable equilibrium. Once again, a match!

What about the case $f'(0)=0$? We are in a similar place to the one for the constant approximation. Here, the linearization gives us the ODE $y'=0$ with stationary solutions only. In the meantime, $f$ might have a semi-stable equilibrium:

A mismatch!

The answer is the best *quadratic* approximation, i.e., the second Taylor polynomial $T_2$ of $f$ around $y=0$, provided $f$ is twice differentiable. However, there may still be exceptions that will call for using the *cubic* Taylor polynomial $T_3$ of $f$. And so on. We will need all the Taylor polynomials, i.e., the Taylor *series*. The idea is explained in Chapter 23.

## Motion under forces: ODEs of second order

We have already seen the simplest ODEs for motion: if we know the *velocity* of a moving object, we can predict its dynamics. For example, the ODE:
$$y'=at+v_0,$$
produces uniformly accelerated motion.

Suppose instead what we know is the *forces* affecting the object and, therefore, its acceleration. How can we predict its dynamics?

Let's review the *discrete model*. We start with the following four quantities that come from the setup of the motion:

- the initial time $t_0$,
- the initial velocity $v_0$,
- the initial location $p_0$, and
- the current acceleration $a_1$.

As we progress in time and space, new numbers are placed in the next row of our spreadsheet: $$t_n,\ a_n,\ v_n,\ p_n,\ n=1,2,3,...$$ with the following recursive formulas:

- $t_{n+1}=t_n+\Delta t$;
- $v_{n+1}=v_n+a_n\cdot \Delta t$;
- $p_{n+1}=p_n+v_n\cdot \Delta t$.

The result is a growing table of values: $$\begin{array}{c|c|c|c|c|c} &\text{iteration } n&\text{time }t_n&&\text{acceleration }a_n&\text{velocity }v_n&\text{location }p_n\\ \hline \text{initial:}&0&3.5&&--&33&22\\ &1&3.6&&66&38.5&25.3\\ &...&...&&...&...&...\\ &1000&103.5&&666&4&336\\ &...&...&&...&...&...\\ \end{array}$$

Where does the current acceleration come from? Our main interest is the case when the acceleration depends on the last location, e.g., $a_{n+1}=1/p_n^2$, such as when the gravity depends on the distance to the planet or the force of the spring depends on the distance of its end from the equilibrium.

**Example.** Recall some examples. A rolling ball is unaffected by horizontal forces and the recursive formulas for the horizontal motion simplify as follows:

- the velocity $v_{n+1}=v_n+a_n\cdot \Delta t=v_n=v_0$ is constant;
- the position $p_{n+1}=p_n+v_n\cdot \Delta t=p_n+v_0\cdot \Delta t$ grows at equal increments.

In other words, the position depends linearly on the time. A falling ball is unaffected by horizontal forces and the vertical force is constant: $a_n=a$ for all $n$. The first of the two recursive formulas for the vertical motion simplifies as follows:

- the velocity $v_{n+1}=v_n+a_n\cdot \Delta t=v_n+a\cdot \Delta t$ grows at equal increments;
- the position $p_{n+1}=p_n+v_n\cdot \Delta t$ grows at linearly increasing increments.

$\square$

*Problem:* Given

- the acceleration as a function of the location $z=f(y)$,

represent

- the velocity as a function of time $v=v(t)$ and
- the location as a function of time $y=y(t)$.

Then our two functions have to satisfy:
$$a_n=f(p_n),\ v_n=v(t_n), \text{ and } p_n=y(t_n).$$
We assume that there is a version of our pair recursive relations,
$$\begin{array}{lll}
v_{n+1} &=v_n&+a_n\cdot \Delta t,\\
p_{n+1} &=p_n&+v_n\cdot \Delta t,
\end{array}$$
for *every* $\Delta t>0$ small enough.
We substitute these two, as well as $t=t_n$, into our recursive formulas:
$$\begin{array}{lll}
v(t+\Delta t)&=v(t)&+f(y(t+\Delta t))\cdot \Delta t,\\
y(t+\Delta t)&=y(t)&+v(t+\Delta t)\cdot \Delta t.
\end{array}$$
Then,
$$\begin{array}{lll}
\frac{v(t+\Delta t)-v(t)}{\Delta t}&=f(y(t+\Delta t)),\\
\frac{y(t+\Delta t)-y(t)}{\Delta t}&=v(t+\Delta t).
\end{array}$$
Taking the limit over $\Delta t\to 0$ gives us the following relations between our functions:
$$\begin{array}{lll}
v'(t)&=f(y(t)),\\
y'(t)&=v(t),
\end{array}$$
provided $y=y(t)$ and $v=v(t)$ are differentiable at $t$ and $z=f(y)$ is continuous at $y(t)$.

The outcome is treated either as a *system of differential equations of first order* discussed later:
$$\begin{cases}
y'&=v;\\
v'&=f(y),
\end{cases}$$
which is a vector field, or as a single *differential equation of second order*
$$y' '=f(y).$$

**Example (free fall).** Description: “the acceleration is constant”. Suppose $\Delta t$ is again the increment of time. Then the increment of the velocity $\Delta v$ is proportional to $\Delta t$ as in this equation:
$$\Delta v=k\cdot \Delta t.$$
Suppose we have an initial condition for the velocity:
$$v(t_0)=v_0.$$
Then the equation gives us this recursive formula:
$$v(t_{n+1})=v(t_n)+k\cdot \Delta t.$$
What about the location, $y$? We have this equation for any interval of time:
$$\Delta y=v \cdot \Delta t,$$
where $v$ is known from above. Suppose we also have an initial condition for the location:
$$y(t_0)=y_0.$$
This sets up a recursive formula:
$$y(t_{n+1})=y(t_n)+v(t_n)\cdot \Delta t.$$
The combination of the two recursive formulas gives us a discrete model we have used in Parts I and II to model free fall and generally acquire location from acceleration:

In the meantime, if we think of $y$ and $v$ as sampled differentiable functions, the two equations are converted to ODEs as follows: $$\Delta y=v\cdot \Delta t \ \Longrightarrow\ \frac{\Delta y}{\Delta t}=v \ \Longrightarrow\ \lim_{\Delta t\to 0}\frac{\Delta y}{\Delta t}=v \ \Longrightarrow\ y'=v,$$ and $$\Delta v=k\cdot \Delta t \ \Longrightarrow\ \frac{\Delta v}{\Delta t}=k \ \Longrightarrow\ \lim_{\Delta t\to 0}\frac{\Delta v}{\Delta t}=k \ \Longrightarrow\ v'=k.$$ Therefore, we have a single ODE: $$y' '=k.$$ Its solution is a quadratic function: $$y(t)=kt^2/2+ pt+q,$$ where $p,q$ are two parameters to be determined from the initial conditions. $\square$

**Example (oscillating spring).** Description: “the acceleration of an object on a spring is proportional to the negative of the current location”.

The increment of the velocity $\Delta v$ is proportional to the location $y$ and to $\Delta t$ as in this equation: $$\Delta v=-k\cdot y\cdot \Delta t,\ k>0.$$ Suppose we have an initial condition for the velocity: $$v(t_0)=v_0.$$ Then the equation gives us this recursive formula: $$v(t_{n+1})=v(t_n)-ky(t_n)\cdot \Delta t.$$ For the location, we have this equation for any interval of time: $$\Delta y=v \cdot \Delta t,$$ where $v$ is known from above. Suppose we also have an initial condition for the location: $$y(t_0)=y_0.$$ This sets up a recursive formula: $$y(t_{n+1})=y(t_n)+v(t_n)\cdot \Delta t.$$ The combination of the two recursive formulas gives us a discrete model:

Finally, if we think of $y$ and $v$ as sampled differentiable functions, the two equations are converted to ODEs as follows: $$y'=v, \text{ and } v'=-ky.$$ Therefore, we have a single ODE: $$y' '=-ky,\ k>0.$$ Its solution is a trig function: $$y(t)=A\cos \sqrt{k}t+B\sin \sqrt{k}t,$$ where $A,B$ are two parameters to be determined from the initial conditions. $\square$