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 9 Ordinary Differential Equations9.1 Differential equationsA differential equation (DE) is an equation that describes the derivatives of an unknown function. “Solving a DE” means finding a function whose derivatives satisfy the equation. For example, when bacteria grow in particularly bacteriafriendly 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 f to be a function that maps from time, t, to population y, such that y = f(t). We don’t know what it is, but we can write a differential equation that describes it:
where a is a constant that characterizes how quickly the population increases. 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 ordinary differential equation (ODE) because all the derivatives involved are taken with respect to the same variable. If the equation related derivatives with respect to different variables (partial derivatives), it would be a partial differential equation. This equation is first order because it involves only first derivatives. If it involved second derivatives, it would be second order, and so on. This equation is linear because each term involves t, f, or df/dt raised to the first power; if any of the terms involved products or powers of t, f, and df/dt it would be nonlinear. Linear, first order ODEs can be solved analytically; that is, we can express the solution as a function of t. This particular ODE has an infinite number of solutions, but they all have this form:
For any values of a and b, the function (3) satisfies the ODE (1). If you don’t believe me, take its derivative and check. 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 f(0) = 5 billion cells, then we can write
and solve for b, which is 5. That determines what we wanted to know:
The extra bit of information that determines b is called the initial condition (although it isn’t always specified at t=0). 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. 9.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 t and measure the current population, y, and the rate of change, r. What do you think the population will be after some period of time Δ t has elapsed? If you said y + r Δ t, congratulations! You just invented Euler’s method. This estimate is based on the assumption that r is constant, but in general it’s not, so we only expect the estimate to be good if r changes slowly and Δ t is small. But let’s assume (for now) that the ODE we are interested in can be written so that
where g is some function that maps (t, y) onto r; that is, given the time and current population, it computes the rate of change. Then we can advance from one point in time to the next using these equations:
Here {T_{i}} is a sequence of times where we estimate the value of f, and {F_{i}} is the sequence of estimates. For each index i, F_{i} is an estimate of f(T_{i}). The interval Δ t is called the time step. Assuming that we start at t=0 and we have an initial condition f(0) = y_{0} (where y_{0} denotes a particular, known value), we set T_{0} = 0 and F_{0} = y_{0}, and then use equations (7) and (8) to compute values of T_{i} and F_{i} until T_{i} gets to the value of 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 Δ t and then run it again with time step Δ t/2. If the results are the same, they are probably accurate; otherwise, cut the time step again. Euler’s method is first order, which means that each time you cut the time step in half, you expect the estimation error to drop by half. With a secondorder method, you expect the error to drop by a factor of 4; thirdorder drops by 8, etc. The price of higher order methods is that they have to evaluate g more times per time step. 9.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 and their meanings, which are presented according to their types:
So f is a function that computes the population as a function of time, f(t) is the function evaluated at a particular time, and if we assign f(t) to a variable, we usually call that variable y. Similarly, g is a “rate function” that computes the rate of change as a function of time and population. df/dt is the first derivative of f, and it maps from t to a rate. If we assign df/dt(t) to a variable, we can call that variable r. It is easy to get df/dt confused with g, but notice that they are not even the same type. df/dt is more specific: it computes the rate of change for the population f(t) at any time t; g is more general: it computes the rate of change for any population y at any time t. 9.4 ode45A limitation of Euler’s method is that the time step is constant from one iteration to the next. But some parts of the solution are harder to estimate than others; if the time step is small enough to get the hard parts right, it is doing more work than necessary on the easy parts. The ideal solution is to adjust the time step as you go along. Methods that do that are called adaptive, and one of the best adaptive methods is the DormandPrince pair of RungeKutta formulas. You don’t have to know what that means, because the nice people at Mathworks have implemented it in a function called ode45. The ode stands for “ordinary differential equation [solver]” and the 45 indicates that it uses a combination of 4th and 5th order formulas. In order to use ode45, you have to write a MATLAB function that evaluates g as a function of t and y. As an example, suppose that the rate of population growth for rats depends on the current population and the availability of food, which varies over the course of the year. The governing equation might be something like
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 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)
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:
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 xaxis shows time from 0 to 365 days; the yaxis 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. 9.5 Multiple output variablesode45 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(end1) is the secondtolast 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? 9.6 Analytic or numerical?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 f(t). 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 the output with the simple function above: This figure shows the output of a more complicated script that zooms 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}, F_{i}). Notice that ode45 typically evaluates g several times for each estimate. This allows it to improve the estimates, for one thing, but also to detect places where the errors are increasing so it can decrease the time step (or the other way around). 9.7 What could go wrong?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) Not enough input arguments. Error in rats (line 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:
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_8_7(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 Output argument "res" (and maybe others) not assigned during call to "rats_8_7". Error in odearguments (line 87) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0. Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin); 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 odearguments (line 21) When the first argument to ode45 is a function handle, the tspan argument must have at least two elements. Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin); 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). 9.8 StiffnessThere 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:
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. Exercise 1
Write a rate function for this ODE and use
ode45 to solve it with the given initial condition and interval.
Start with δ = 0.1 and decrease it by multiples of 10. If
you get tired of waiting for a computation to complete, you can
press the Stop button in the Figure window or press ControlC in
the Command Window.
Now replace ode45 with ode23s and try again! 9.9 Glossary
9.10 ExercisesExercise 2
Suppose that you are given an 8 ounce cup of coffee at 90 ^{∘}C and
a 1 ounce container of cream at room temperature, which is 20 ^{∘}C.
You have learned from bitter experience that the hottest coffee you
can drink comfortably is 60 ^{∘}C.
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:
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.^{2} 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.

Are you using one of our books in a class?We'd like to know about it. Please consider filling out this short survey. 