This HTML version of the book is provided as a convenience, but some math equations are not translated correctly. The PDF version is more reliable. ## Chapter 8 Ordinary Differential Equations## 8.1 Differential equationsA For example, when bacteria grow in particularly bacteria-friendly
conditions, the rate of growth at any point in time is proportional to
the current population. What we might like to know is the population
as a function of time. Toward that end, let’s define
where Notice that both sides of the equation are functions. To say that two functions are equal is to say that their values are equal at all times. In other words:
This is an This equation is This equation is Linear, first order ODEs can be solved analytically; that is, we
can express the solution as a function of
For any value of If we know the population of bacteria at a particular point in time, we can use that additional information to determine which of the infinite solutions is the (unique) one that describes a particular population over time. For example, if we know that
and solve for
The extra bit of information that determines Unfortunately, most interesting physical systems are described by nonlinear DEs, most of which can’t be solved analytically. The alternative is to solve them numerically. ## 8.2 Euler’s methodThe simplest numerical method for ODEs is Euler’s method. Here’s a
test to see if you are as smart as Euler. Let’s say that you
arrive at time If you said This estimate is based on the assumption that But let’s assume (for now) that the ODE we are interested in can be written so that
where
Here { f, and {F} is the sequence of estimates. For each
index _{i}i, F is an estimate of _{i}f(T).
The interval Δ _{i}t is called the time step.Assuming that we start at F until _{i}T
gets to the value of _{i}t we are interested in.If the rate doesn’t change too fast and the time step isn’t
too big, Euler’s method is accurate enough for most purposes. One
way to check is to run it once with time step Δ Euler’s method is ## 8.3 Another note on notationThere’s a lot of math notation in this chapter so I want to pause to review what we have so far. Here are the variables, their meanings, and their types:
So Similarly,
It is easy to get ## 8.4 |

| (t) = a f(t) | ⎡ ⎣ | 1 + sin(ω t) | ⎤ ⎦ |

where *t* is time in days and *f*(*t*) is the population at time *t*.

*a* and ω are **parameters**. A parameter is a value that
quantifies a physical aspect of the scenario being modeled. For
example, in Exercise 6.2 we used parameters `rho` and
`r` to quantify the density and radius of a duck. Parameters are
often constants, but in some models they vary in time.

In this example, *a* characterizes the reproductive rate, and
ω is the frequency of a periodic function that describes
the effect of varying food supply on reproduction.

This equation specifies a relationship between a
function and its derivative. In order to estimate values of
*f* numerically, we have to transform it into a rate function.

The first step is to introduce a variable,
*y*, as another name for *f*(*t*)

| (t) = a y | ⎡ ⎣ | 1 + sin(ω t) | ⎤ ⎦ |

This equation means that if we are given *t* and *y*, we can
compute *df*/*dt*(*t*), which is the rate of change of *f*.
The next step is to express that computation as a function called
*g*:

g(t, y) = a y | ⎡ ⎣ | 1 + sin(ω t) | ⎤ ⎦ |

Writing the function this way is useful because we can use it
with Euler’s method
or `ode45` to estimate values of *f*. All we have to
do is write a MATLAB function that evaluates *g*. Here’s what
that looks like using the values *a* = 0.01
and ω = 2 π/365 (one cycle per year):

function res = rats(t, y) a = 0.01; omega = 2 * pi / 365; res = a * y * (1 + sin(omega * t)); end

You can test this function from the Command Window by calling it with
different values of `t` and `y`; the result is the rate of
change (in units of rats per day):

>> r = rats(0, 2) r = 0.0200

So if there are two rats on January 1, we expect them to reproduce at a rate that would produce 2 more rats per hundred days. But if we come back in April, the rate has almost doubled:

>> r = rats(120, 2) r = 0.0376

Since the rate is constantly changing, it is not easy to predict
the future rat population, but that is exactly what `ode45` does.
Here’s how you would use it:

>> ode45(@rats, [0, 365], 2)

The first argument is a handle for the function that
computes *g*. The second argument is the interval we are interested
in, one year. The third argument is the initial population, *f*(0) = 2.

When you call `ode45` without assigning the result to a variable,
MATLAB displays the result in a figure:

The x-axis shows time from 0 to 365 days; the y-axis shows the rat population, which starts at 2 and grows to almost 80. The rate of growth is slow in the winter and summer, and faster in the spring and fall, but it also accelerates as the population grows.

`ode45` is one of many MATLAB functions that return more
than one output variable. The syntax for calling it and saving
the results is

>> [T, Y] = ode45(@rats, [0, 365], 2);

The first return value is assigned to `T`; the second is assigned
to `Y`. Each element of `T` is a time,
*t*, where `ode45` estimated the population; each element of `Y` is an estimate of *f*(*t*).

If you assign the output values to variables,
`ode45` doesn’t draw the figure;
you have to do it yourself:

>> plot(T, Y, 'bo-')

If you plot the elements of `T`, you’ll see that the
space between the points is not quite even. They are closer
together at the beginning of the interval and farther apart at the end.

To see the population at the end of the year, you can display the last element from each vector:

>> [T(end), Y(end)] ans = 365.0000 76.9530

`end` is a special word in MATLAB; when it appears as an index,
it means “the index of the last element.” You can use it in an
expression, so `Y(end-1)` is the second-to-last element of
`Y`.

How much does the final population change if you double the initial
population? How much does it change if you double the interval
to two years? How much does it change if you double the value
of *a*?

When you solve an ODE analytically, the result is a function, *f*,
that allows you to compute the population, *f*(*t*), for any value of
*t*. When you solve an ODE numerically, you get two vectors. You can
think of these vectors as a discrete approximation of the continuous
function *f*: “discrete” because it is only defined for certain
values of *t*, and “approximate” because each value *F _{i}*
is only an estimate of the true value

So those are the limitations of numerical solutions. The primary advantage is that you can compute numerical solutions to ODEs that don’t have analytic solutions, which is the vast majority of nonlinear ODEs.

If you are curious to know more about how `ode45` works, you
can modify `rats` to display the points, (*t*, *y*), where
`ode45` evaluates *g*. Here is a simple version:

function res = rats(t, y) plot(t, y, 'bo') a = 0.01; omega = 2 * pi / 365; res = a * y * (1 + sin(omega * t)); end

Each time `rats` is called, it plots one data point; in order
to see all of the data points, you have to use `hold on`.

>> clf; hold on >> [T, Y] = ode45(@rats, [0, 10], 2);

This figure shows part of the output, zoomed in on the range from Day 100 to 170:

The circles show the points where `ode45` called `rats`.
The lines through the circles show the slope (rate of change) calculated
at each point. The rectangles show the locations of the estimates
(*T _{i}*,

Don’t forget the `@` on the function handle.
If you leave it out, MATLAB treats the first argument as a function
call, and calls `rats` without providing arguments.

>> ode45(rats, [0,365], 2) ??? Input argument "y" is undefined. Error in ==> rats at 4 res = a * y * (1 + sin(omega * t));

Again, the error message is confusing, because it looks like the problem
is in `rats`. You’ve been warned!

Also, remember that the function you write will be called by
`ode45`, which means it has to have the signature `ode45`
expects: it should take two input variables, `t` and `y`,
in that order, and return one output variable, `r`.

If you are working with a rate function like this:

g(t, y) = a y |

You might be tempted to write this:

function res = rate_func(y) % WRONG a = 0.1 res = a * y end

But that would be wrong. So very wrong. Why? Because
when `ode45` calls `rate_func`, it provides two arguments.
If you only take one input variable, you’ll get an error. So
you have to write a function that takes `t` as an input
variable, even if you don’t use it.

function res = rate_func(t, y) % RIGHT a = 0.1 res = a * y end

Another common error is to write a function that doesn’t make an assignment to the output variable. If you write something like this:

function res = rats(t, y) a = 0.01; omega = 2 * pi / 365; r = a * y * (1 + sin(omega * t)) % WRONG end

And then call it from `ode45`, you get

>> ode45(@rats, [0,365], 2) ??? Output argument "res" (and maybe others) not assigned during call to "/home/downey/rats.m (rats)". Error in ==> rats at 2 a = 0.01; Error in ==> funfun/private/odearguments at 110 f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0. Error in ==> ode45 at 173 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

This might be a scary message, but if you read the first line and ignore the rest, you’ll get the idea.

Yet another mistake that people make with `ode45` is leaving
out the brackets on the second argument. In that case, MATLAB
thinks there are four arguments, and you get

>> ode45(@rats, 0, 365, 2) ??? Error using ==> funfun/private/odearguments When the first argument to ode45 is a function handle, the tspan argument must have at least two elements. Error in ==> ode45 at 173 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Again, if you read the first line, you should be able to figure
out the problem (`tspan` stands for “time span”, which we
have been calling the interval).

There is yet another problem you might encounter, but if it makes you
feel better, it might not be your fault: the problem you are trying to
solve might be **stiff**^{1}.

I won’t give a technical explanation of stiffness here, except to say that for some problems (over some intervals with some initial conditions) the time step needed to control the error is very small, which means that the computation takes a long time. Here’s one example:

| = f^{2} − f^{3} |

If you solve this ODE with the initial condition *f*(0) = δ over
the interval from 0 to 2/δ, with δ = 0.01, you should
see something like this:

After the transition from 0 to 1, the time step is very small and the computation goes slowly. For smaller values of δ, the situation is even worse.

In this case, the problem is easy to fix: instead of `ode45` you can
use `ode23s`, an ODE solver that tends to perform well on stiff
problems (that’s what the “s” stands for).

In general, if you find that `ode45` is taking a long time, you
might want to try one of the stiff solvers. It won’t always solve
the problem, but if the problem is stiffness, the improvement can
be striking.

*Now replace ode45 with ode23s and try again!
*

**differential equation (DE):**- An equation that relates the derivatives of an unknown function.
**ordinary DE:**- A DE in which all derivatives are taken with respect to the same variable.
**partial DE:**- A DE that includes derivatives with respect to more than one variable
**first order (ODE):**- A DE that includes only first derivatives.
**linear:**- A DE that includes no products or powers of the function and its derivatives.
**time step:**- The interval in time between successive estimates in the numerical solution of a DE.
**first order (numerical method):**- A method whose error is expected to halve when the time step is halved.
**adaptive:**- A method that adjusts the time step to control error.
**stiffness:**- A characteristic of some ODEs that makes some ODE
solvers run slowly (or generate bad estimates). Some ODE solvers,
like
`ode23s`, are designed to work on stiff problems. **parameter:**- A value that appears in a model to quantify some physical aspect of the scenario being modeled.

*Assuming that you take cream in your coffee, and that you would like
to start drinking as soon as possible, are you better off adding
the cream immediately or waiting? And if you should wait, then how
long?*

*To answer this question, you have to model the cooling process
of a hot liquid in air. Hot coffee transfers heat to the environment
by conduction, radiation, and evaporative cooling. Quantifying these
effects individually would be challenging and unnecessary to answer
the question as posed.*

*As a simplification, we can use Newton’s Law of
Cooling ^{2}:*

| = −r (f − e) |

*
where f is the temperature of the coffee as a function of time and
df/dt is its time derivative; e is the temperature of the
environment, which is a constant in this case, and r is a parameter
(also constant) that characterizes the rate of heat transfer.*

*It would be easy to estimate r for a given coffee cup by making
a few measurements over time. Let’s assume that that has been
done and r has been found to be 0.001 in units of inverse
seconds, 1/s.*

*Using mathematical notation, write the rate function,**g*, as a function of*y*, where*y*is the temperature of the coffee at a particular point in time.*Create an M-file named*`coffee`and write a function called`coffee`that takes no input variables and returns no output value. Put a simple statement like`x=5`in the body of the function and invoke`coffee()`from the Command Window.*Add a function called*`rate_func`that takes`t`and`y`and computes*g*(*t*,*y*). Notice that in this case*g*does not actually depend on*t*; nevertheless, your function has to take*t*as the first input argument in order to work with`ode45`.*Test your function by adding a line like*`rate_func(0,90)`to`coffee`, the call`coffee`from the Command Window.*Once you get*`rate_func(0,90)`working, modify`coffee`to use`ode45`to compute the temperature of the coffee (ignoring the cream) for 60 minutes. Confirm that the coffee cools quickly at first, then more slowly, and reaches room temperature (approximately) after about an hour.*Write a function called*`mix_func`that computes the final temperature of a mixture of two liquids. It should take the volumes and temperatures of the liquids as parameters.*In general, the final temperature of a mixture depends on the specific heat of the two substances*^{3}. But if we make the simplifying assumption that coffee and cream have the same density and specific heat, then the final temperature is (*v*_{1}*y*_{1}+*v*_{2}*y*_{2}) / (*v*_{1}+*v*_{2}), where*v*_{1}and*v*_{2}are the volumes of the liquids, and*y*_{1}and*y*_{2}are their temperatures.*Add code to*`coffee`to test`mix_func`.*Use*`mix_func`and`ode45`to compute the time until the coffee is drinkable if you add the cream immediately.*Modify*`coffee`so it takes an input variable*t*that determines how many seconds the coffee is allowed to cool before adding the cream, and returns the temperature of the coffee after mixing.*Use*`fzero`to find the time*t*that causes the temperature of the coffee after mixing to be 60^{∘}C.*What do these results tell you about the answer to the original question? Is the answer what you expected? What simplifying assumptions does this answer depend on? Which of them do you think has the biggest effect? Do you think it is big enough to affect the outcome? Overall, how confident are you that this model can give a definitive answer to this question? What might you do to improve it?*