Numerical Methods
In this course, we have solved by hand certain types of differential equations (linear equations, separable equations, and exact equations) using analytic techniques. Unfortunately, the vast majority of differential equations of the form ^{dy}⁄_{dt} = ƒ(t,y) are not linear, separable, or exact. For this reason, we need to think about other possible ways to analyze firstorder initial value problems. The techniques for finding approximate solutions to differential equation problems using estimation and calculation rather than analytic manipulation are collectively called numerical methods. When we use numerical methods, we don't get a nice formula as a solution to our initial value problem. We do, however, get a way to calculate an approximate value for the solution at any point we choose.
Consider the following initial value problem:
(1) 
^{dy}⁄_{dx} = x^{2}  y^{3} y(0) = 0 
The differential equation in (1) is not linear (because of the y^{3}), nor is it separable or exact. Convince yourself of this with a pencil and paper.
We have no symbolic methods for solving this equation. Thus we will need to use the numerical methods mentioned above. One way we can get an idea of what the solutions to this (or any) differential equation look like is with a direction field. We've already seen in Assignment 2 how to use slopefield
and drawode
to get a picture of our solution.
>> f = @(x,y) x^2  y^3;
slopefield(f, [10,10], [10,10], 30)
hold on
drawode(f, [10,10], 0, 0)
hold off
When you run these commands, you will get a warning: MATLAB says it is "unable to meet integration tolerances without reducing the step size." In this case, the curve becomes so steep on the left side that the computer cannot accurately compute more points to plot in a reasonable time. What does that suggest to us about this solution curve?
Pictures like these are very informative for understanding the behavior of the solutions to even very complicated differential equations. Unfortunately, such graphical techniques are limited to three dimensions at best. Indeed, most technology problems involve highdimensional systems of ODEs, so that it is impossible to visualize their solutions. Even with two or threedimensional systems, visual techniques might fail us when we want to get very precise information. For example, we might want to find out the value of y in our solution above when x = 10. From the direction field, we can guess that it is somewhere around 4.3. By changing the axes in slopefield
, we might be able to zoom in and improve this estimate, but that can be an awkward way to obtain a highly accurate answer.
Thinking about direction fields and their limitations leads us to two questions.
 Q1. What do we do for large systems?
 Q2. How does
drawode
actually compute solution curves? After all, it is dangerous to run programs without having a solid idea how they work; if you do, eventually you will make a serious mistake.
In this lab we shall try to give you the basic idea behind answers to both Q1 and Q2, although going there from our sketch here is a pretty big leap. (Numerical analysis is, after all, an entire branch of mathematics!) We will explore a couple of numerical methods, beginning with a relatively simple method called Euler's Method. It is the most easilyunderstood example of a numerical approach to solving differential equations.
Euler's Formula: A Numerical Method
Euler's Method is one of the simplest and oldest numerical methods for approximating solutions to differential equations that cannot be solved with a nice formula. Euler's Method is also called the tangent line method, and in essence it is an algorithmic way of plotting an approximate solution to an initial value problem through the direction field. It can also be used to get an estimated value for that solution at a given value for x. Here is how it works:
Say we have a generic initial value problem of the form ^{dy}⁄_{dx} = f(x,y) with y(x_{0}) = y_{0}. We may not be able to find a formula for the solution y = φ(x), but if f(x,y) is a continuous function, we do know that such a solution exists. (See Theorem 2.8.1 in Boyce & DiPrima.) We can therefore talk about the solution y = φ(x) and its properties. Our goal is to say as much as we can about this solution, even though we may not be able to write down a formula for it.
At any pair of values for x and y, the direction field gives us the slope of the tangent, and therefore the direction, of the solution curve y = φ(x) passing through that point. We know that the point (x_{0 }, y_{0}) lies on the solution curve, and at that point, the slope of the tangent is f(x_{0 }, y_{0}). Therefore, the equation of the tangent line to the solution at the point (x_{0 }, y_{0}) is:
y = y_{0} + f(x_{0 }, y_{0})(x  x_{0}).
Since a tangent line is a good approximation to a curve near the point (x_{0 }, y_{0}), we can move a short distance along the tangent line by increasing x by some small amount h; then we arrive at a new point, which we will call (x_{1 }, y_{1}). That is, we choose x_{1} = x_{0} + h, and then we obtain y_{1} by plugging x_{1} into the equation for the tangent line. Since we arrived at this new point by moving from (x_{0 }, y_{0}) a small distance along the tangent line at that point, we can say that the point (x_{1 }, y_{1}) is almost on the solution curve y = φ(x).
We can now find the tangent line to the solution at this new point, y = y_{1} + f(x_{1 }, y_{1})(x  x_{1}), and then increase x by the small amount h again, and come up with a new point (x_{2 }, y_{2}). Continuing with this method, we use the value of y found at each step to calculate the next value of y. This whole process can be stated more succinctly by Euler's Formula:
y_{n+1} = y_{n} + f(x_{n }, y_{n})(x_{n+1}  x_{n}).
Another way of writing this is as follows:
x_{n+1} = x_{n} + h;
y_{n+1} = y_{n} + f(x_{n }, y_{n}) · h.
A more complete description of the formula is given in the textbook. As a side remark, note that a natural variant of this algorithm changes the x equation to
x_{n+1} = x_{n} + g(x_{n }, y_{n}) · h,
which is general enough to handle firstorder twodimensional systems.
To see Euler's Method in use, let's try an example.
Consider the following initial value problem:
(2) 
^{dy}⁄_{dx} = x^{2}  1 y(0) = 1 
Suppose we want to use Euler's Method to graph an estimate for the initial value problem with f(x,y) = x^{2}  1 given above, over the interval 0 ≤ x ≤ 2. From the initial value condition, we know that when x = 0, the value of y is 1. Hence we will start at the initial point (x_{0 }, y_{0}) = (0,1). The tangent line at this point is y = 1  x. If we use a "step size" of h = 1, then our x coordinates will be x_{0} = 0, x_{1} = 1, and x_{2} = 2. Euler's Method will give us estimates for the y values corresponding to each of the x values. For this step size, Euler's Method takes just two steps:
x_{0} = 0, x_{1} = 1, x_{2} = 2, 
y_{0} = 1 y_{1} = 1 + f(0,1) · 1 = 0 y_{2} = 0 + f(1,0) · 1 = 0 
So for h = 1, Euler's Method is estimating that our solution curves goes through the three points (0,1), (1,0), and (2,0).
What if we used a smaller step size, say, h = 0.5? Then we would get y values corresponding to x_{0} = 0, x_{1} = 0.5, x_{2} = 1, x_{3} = 1.5, and x_{4} = 2:
x_{0} = 0, x_{1} = 0.5, x_{2} = 1, x_{3} = 1.5, x_{4} = 2, 
y_{0} = 1 y_{1} = 1 + f(0,1) · (0.5) = 0.5 y_{2} = 0.5 + f(0.5, 0.5) · (0.5) = 0.125 y_{3} = 0.125 + f(1, 0.125) · (0.5) = 0.125 y_{4} = 0.125 + f(1.5, 0.125) · (0.5) = 0.75 
With a step size of h = 0.5, we now get the five points (0,1), (0.5, 0.5), (1, 0.125), (1.5, 0.125), and (2, 0.75).
Keep this example in mind; we will soon verify these calculations using MATLAB.
Implementing Euler's Formula in MATLAB
We are going to execute Euler's Method in MATLAB using MFiles. It is possible to simply type the necessary commands into the MATLAB window, but this is inefficient. As we saw back in Assignment 1, the advantage to using MFiles instead is that each time we want to run Euler's Method, we can do so with a single command, rather than retyping a long list of commands. In the same way that we defined new functions with MFiles before, we will now define the (more complicated) function Euler
. The content of Euler
may seem daunting at first, but we will return to explaining how it works after a quick application.
We will need an MFile that runs Euler's Method, so in MATLAB, use the button New or New Script. This will open a new Mfile. In this MFile, we will create a function Euler
that will take five parameters: h
, x0
, y0
, interval_length
, and func
. The first four parameters should be numerical values, where h
is a number greater than 0. The last parameter func
should be a function of two variables, which we will identify as x and y. Using these inputs, Euler
will apply Euler's Method to the differential equation ^{dy}⁄_{dx} = func
on the interval [x0
, x0
+ interval_length
], using a step size of h
, going through the point (x0
, y0
).
When the MFile editor appears, enter (or copy and paste) the following blue text:
1 2 3 4 5 6 7 8 9 10 

Now select the Save button and save the file as Euler.m. (This name will most likely appear in the box for you.) The file should be saved in a directory where MATLAB can access it. If you are working on an ETS machine, you should be able to use your My Documents folder. If you are not working on an ETS machine, make sure you are saving the file into MATLAB's current directory on your machine. (Remember that if you choose to work on nonETS computers, you do so at your own risk. The TAs cannot provide technical support anywhere outside the designated MATLAB lab.)
We can now use the command Euler
from our MATLAB command window. In the previous section above, we went through Euler's Formula for the initial value problem (2) on the interval 0 ≤ x ≤ 2, with h = 0.5. Let's run this again using our new MATLAB command. In the MATLAB command window, type:
f = @(x,y) x^2  1
This defines our function ^{dy}⁄_{dx} = f(x,y) = x^{2}  1 as a function of two variables, x and y. Now let's run an iteration of Euler's Method:
>> h = 0.5;
[x,y] = Euler(h, 0, 1, 2, f);
[x,y]
The results from running Euler's Method are contained in two arrays, x and y. When we enter the last command [x,y]
(note the absence of a semicolon), MATLAB outputs the x and y coordinates of the points computed by Euler's Method. Note that for this example, the output matches what we got at the end of Example 3.2. This is a good sign; Euler
can compute the points predicted by Euler's Method without the need to do computations by hand. We can also plot these points graphically by typing:
>> plot(x, y, 'r')
hold on
The 'r'
tells MATLAB to draw the graph in red, and as we've seen before, the hold on
command will keep the existing graph while you plot more graphs on top of it in the following exercise.
 Run Euler's Method again, this time for h = 0.2. Plot the results in blue this time. (Use
'b'
instead of'r'
.) Look at the arrays x and y sidebyside by typing in
>>Copy the long table of data that you just produced into your Word document. What estimates did the Euler Method come up with for y at x = 1 and x = 2 using this smaller value of h? Now run Euler's Method for h = 0.1 and again for h = 0.01, and plot the results in green ([x,y]
'g'
) and cyan ('c'
), respectively. What are the estimates for y at x = 1 and x = 2 for these values of h? Write these down in your Word document, along with your graph with all four colored curves. (You need not include the data tables for these last two values of h.)  The differential equation given in (2) is separable. Calculate the solution to the initial value problem by hand, and use MATLAB (or a calculator) to compute the actual values for y at x = 1 and x = 2. Make sure to include your manual work in your Word document. (You can simply type it.)
 Does Euler's estimate appear to give better or worse estimates for the solution as h decreases? Give a brief reasoning for why this might be the case. Does Euler's estimate appear to give better or worse estimates for the actual solution as you move farther from the initial point?
Now that we have seen Euler
in action, let's return to examining the content of the MFile Euler.m. We have already explained the first line, where we defined the parameters our function takes. To see the meaning behind the third and fourth lines, type:
>> x = zeros(10,1);
y = zeros(10,1);
[x,y]
Thus we can see that the third and fourth lines of our MFile zero out the contents of our arrays x and y before we begin. Now, you'll continue examining our code.
In your Word document, briefly explain what is happening in each remaining line of the MFile Euler.m.
The algorithm we used to implement Euler's Method works when we have a single differential equation. It is worth mentioning, however, that Euler's Method readily generalizes to large systems of differential equations.
Experimenting with Euler's Formula
Let's return to our first example:
(1) 
^{dy}⁄_{dx} = x^{2}  y^{3} y(0) = 0 
We noticed before that there is no way for us to solve this system by any analytic method, but
since x^{2}  y^{3} is a continuous function, we know that a solution exists on some interval around the initial value x = 0. We saw using slopefield
and drawode
roughly what the solution curve looks like. Now let's use our new numerical method to estimate the solution.
 Run Euler's Method in MATLAB for the inital value problem (1) on the interval 0 ≤ x ≤ 4. Use different colors to plot the estimates for the following values of h: 1, 0.5, 0.25, and 0.1. Remember to use
hold on
! Some of the colors you can use with MATLAB'splot
command are'r'
,'g'
,'b'
,'c'
,'m'
,'y'
, and'k'
.
Notice that the yaxis is very large. Rescale the yaxis to range from 10 to 10. (Assignment 1 showed us how to do this.) Given the results of the plot, what can you say about the accuracy of the h = 0.1 case? You do not need to include the actual data in the table[x,y]
in your writeup, only the plot and your written answer.  Enter
hold off
so that MATLAB will overwrite your previous plot, and this time, run Euler's Method with the h values 1, 0.1, 0.01, and 0.001. Plot the results on one graph using colors of your choice, and rescale so that y ranges from 0 to 2.5. Paste this plot into your writeup.
Did the additional shrinking of h change the results significantly? What does this suggest about the accuracy of your solution for each step size? Comment on the consistency of your plot with the direction field shown at the top of the page, and also describe any inconsistencies between different step sizes.
This example serves as an excellent warning about some of the hazards that may arise when using numerical methods to solve differential equations. Since the answers are by definition not exact, one should be wary of using numerical results without some way to detect inaccuracy.
ODE45
Happily for our sanity, we do not have to go through the steps above to use numerical methods in MATLAB, because MATLAB has a number of numerical methods built in. One of these is ode45
, which runs a numerical method of a type collectively known as the RungeKutta Methods. These methods are generally more powerful than Euler's Method.
We don't want to go into much detail on the mathematics behind RungeKutta, but there is an informative analogy to Taylor expansions. In Euler's Method, we used the tangent line (which is the first order Taylor polynomial) to approximate our solutions. In essence, RungeKutta Methods use higher order Taylor polynomial approximations.
Let's try to use the ode45
command now.
Consider the differential equation ^{dy}⁄_{dx} = sin(x). We'll first define a function g(x,y) to represent ^{dy}⁄_{dx}:
>>g = @(x,y) sin(x)
Now type the following command into MATLAB:
>>[x,y] = ode45(g, [0, 10], 1)
This tells MATLAB to set y, as a function of x, to be the approximation obtained using ode45
to the solution of the ODE we're looking at on the interval [0,10]. The last parameter, the 1
, gives the initial condition y(0) = 1. Now that we've done this, we can use the plot
command the same way we did before to see a graph of ode45
's estimate of the solution.
The differential equation we just looked at is quite easy to solve symbolically. Let's do that.
 As can be easily verified, the solution to the initial value problem in Example 3.3 is
h(x) = cos(x) + 2Now with the solution of the above differential equation, we can compare the results obtained by usingode45
to the actual values. Type the following commands into MATLAB:
>>
h = @(x) cos(x) + 2; [x, y, h(x), abs(y  h(x))]
ode45
and the actual function value at each specified value of x. How do the values in the two columns compare? The fourth column gives the absolute value of the difference between the function and our approximation; this is the error of our approximation.  Now let us compare the results that
ode45
gives us with the results we would get from using Euler's method. Enter the following commands into MATLAB:
>>
[x, z] = Euler(0.25, 0, 1, 10, g); [y, z, h(x), abs(y  h(x)), abs(z  h(x))]
ode45
; in the second, the results of our Euler's Method routine; and in the third, the values at x of the real solution to our differential equation. Columns four and five give the error ofode45
and of Euler's Method, respectively. Which method seems to give more accurate answers in this situation?
We emphasize two features of ode45
. First, the algorithm is more powerful than Euler's Method. Second, go back and take a careful look at the parameters we plugged into ode45
; there are relatively few of them. Notice that there is no step size parameter, nor is there any other parameter that we may vary to adjust the precision of our calculations. These are picked automatically by the program (a feature of the ode45
implementation, not fundamental to RungeKutta Methods). In Exercise 3.4, we were able to test the results of ode45
explicitly because we could solve the differential equation by hand, but in practice, one should remember to search for inaccuracy before using numerical results.
We can now answer one of the questions asked earlier. At the beginning of this lab, we asked how drawode
works. As a matter of fact, drawode
simply calls ode45
twice, once running forwards in time and once running backwards. So we can think of drawode
as using highdegree polynomials to approximate solutions to initial value problems.
Furthermore, RungeKutta algorithms generalize from a single ODE to a system of ODEs with very little conceptual difference. In fact, the command ode45
can also be used on systems of differential equations. If you're curious, type help ode45
.
Conclusion
The power of computing is best exploited when calculating numerical solutions. In most cases, if a differential equation has a solution, then that solution can be at least closely approximated numerically. One must be a little careful, however. There are many equations that are considered stiff, meaning that some numerical methods for solving them exhibit unstable behavior unless the step size is taken to be very small (as in Exercise 3.3). This is certainly true with Euler's method. The ode45
solver is an allaround tool that does a decent job for many of the equations that come up frequently. In practice, though, one deals with stiff equations by using higherorder numerical methods. This is actually what we just did. Euler's Method is a firstorder method (since we use linear approximations), while ode45
is a fourthorder method. MATLAB also offers other solvers, such as ode15s
or ode23s
. If you're interested, you can type doc
and then search there for ODE solvers to learn more about the specific properties of each numerical solver.
In the next lab, we will examine systems of differential equations in more detail. It is important to know that many numerical methods, including both Euler's Method and ode45
, can handle systems of differential equations with only minor modification.
Last Modified: August 2021