From the series: Solving ODEs in MATLAB
Cleve Moler, MathWorks
ODE1 implements Euler's method. It provides an introduction to numerical methods for ODEs and to the MATLAB suite of ODE solvers. Exponential growth and compound interest are used as examples.
Hello. I'm Cleve Moler, one of the founders and chief mathematician at the MathWorks. This series of videos is about solving ordinary differential equations in MATLAB. We can begin by recalling the definition of derivative.
The derivative of a function at a point is the slope of the tangent line to the graph of the function at that point. Our numerical approximations will rely upon the slope of the secant to the graph. That's a line through two points separated by a distance h.
We'll have a lot to say about the step size h as we go along. What's important to realize is that as h goes to 0, the slope of the secant approaches the slope of tangent. The wiggly equals sign means approximately equal to.
T0 is the point where we are finding the approximation. The value of the derivative at t0 is approximately equal to the slope of the secant. The slope of the secant is the change in the y value over the change in the t value. The change in y value is the difference between the two values of y.
The change in the t value is the step size h. If we rewrite this, we get the value of y at the point t0 plus h is approximately equal to the value of y at t0 plus h times the value of y prime at t0. This is the basis for our first numerical method, Euler's method.
Leonhard Euler was a 18th century Swiss mathematician. Probably the most influential mathematician of his era. He made important contributions to a wide range of fields of mathematics, physics, and astronomy. He invented the notion of a function, for example.
The differential equation is given by this function f of two variables t and y. And the task in general is to find the function y whose derivative is equal to f. Now, there's lots of functions y whose derivative is equal to f. And so there's an initial condition, a point t naught or t0, and a value y0. And the initial condition is that y at t0 should be equal to y0.
Here's some examples. The compound interest problem is just the interest rate times y. Here, the function of t and y doesn't actually depend upon t. And it's linear in y. The initial condition is at time 0, there's a specified value of y, like $100.
That's the compound interest problem. Here's the logistic equation. Nonlinear equation. Here, f of t and y, again, doesn't depend upon t. And it's a constant times y minus another constant times y squared. That's the logistic equation.
And again, the value is specified at 0. Let's say y at 0 is equal to 1. Here's another nonlinear equation. F of t and y is t squared plus y squared. It's not possible to find an analytic solution to this equation.
We'll use these numerical methods to find a solution to this equation. Initial condition y at 0 is equal to 0. That's an example of a function of t and y. Euler's method actually isn't a practical numerical method in general.
We're just using it to get us started thinking about the ideas underlying numerical methods. Euler's method involves a sequence of points t sub n separated by a fixed step size h. And then y sub n is an approximation to the value of the solution at t sub n.
The approximation comes from the slope of the secant. The ratio of the difference of the values of y and to the step size h. The differential equation says that this ratio should be the value of the function at t sub n.
And if we rearrange this equation, we get Euler's method. That yn plus 1 is yn plus h times the function f evaluated at t sub n and y sub n. This is Euler's method. We're now ready for our first MATLAB program, ODE1. It's called ODE1 because it's our first program and because it evaluates the function f that defines the differential equation once per step.
There are five input arguments. The first is f, a function that defines the differential equation. This is something called an anonymous function. I'll talk more about that in a moment. The other four are scalar numerical values.
The first three define the interval of integration. We're going from t0 in steps of h to t final. The fifth input argument is y0, the initial value. The output is a vector. Vector y out is the values of the solution at the points in the interval.
We start by putting y0, the initial value, into y and then putting y into the output vector. The body of the function is a four loop, t goes from T0 not in steps of H up to one step short of t final and then the final passage through the body of the code takes t up to t final.
We evaluate the function f at t and y. That gives us a slope s, s is for slope. Here's the Euler step. Take the current value of y, add h, times the slope. That gives us this new value of y. And then y is appended to y out.
This MATLAB construction with the square brackets takes a vector y, adds another value to it, making it one element longer and puts the resulting y out back in y out. This is the entire code. This is it. This is ODE1 that implements Euler's method.
The first argument to any of the MATLAB ODE solvers is the name of a function that specifies the differential equation. This is known as a function handle. The easiest way to get a function handle is to make use of an anonymous function created with the ampersand or at sign.
All of the differential equations involve anonymous functions of two variables, t and y. And so we have f equals at parenthesis t comma y closed parenthesis. This is followed by any expression involving either t or y.
And many of them don't depend upon t. So here is an anonymous function defining our interest problem. And we can just evaluate this like any ordinary function. When y is equal to 1, f of 1 is 0.06.
Here's an example of a function that depends upon both t and y. The functions can involve constants that have values. So here, we can define two constants. And then we can use those two constants to define the logistic equation f of a times y minus b times y squared.
Again, this is an autonomous equation that doesn't actually depend upon t. Let's see how Euler's method and ODE1 work on this simple example. Y prime equals 2y with the initial condition y of 0 equals 10 on the interval t between 0 and 3.
We define the anonymous function f of t and y is equal to 2y. The initial condition is t0 equals 0. We're going to take a step size of 1. Go to t final equals 3 starting at y0 equals 10 and here's our call to ODE1. We have an animation that shows these steps.
Start at t0 equals 0 and y0 equals 10. Here's our first point. We evaluate the function there. We get a slope of 20. That's 2 times 10. We take an Euler step of length 1 across the first step. That brings us to the second point, 30.
Evaluate the function there. Two times 30 is 60. That's our slope. Take the second step to get to y2. Y2 is 90. Evaluate the function there. Get 2 times 90 is 180. That gives us a slope. Take a step across the interval with that slope would get us to a third point.
The third point is 270 and that's the end of the integration. So that's three Euler steps to get from t0 to t final. Euler's method is actually the same as computing compound interest. So let's do a compound interest problem.
Define the interest rate. Define our anonymous function using that interest rate. Start at time 0. Go in steps of a month. Go for 10 years. Start with $100. And here's our result of using ODE1 to compute compound interest.
That's 121 numbers. MATLAB actually has a format for looking at dollars and cents. And so here they are as dollars and cents starting with $100 and compounding every month. We get up to just over $180.
I'm going to plot that. So I want a time vector months. And I actually want to compare it with simple interest. So here's how you compute simple interest. $0.50 a month. And now let's plot those two.
So the straight line is simple interest getting up to $160. And the blue line is the compound interest. There is a slight upward curvature, getting us up to $180. There's a dot every month here as we show the results of Euler's method, which as I said is the same thing as computing compound interest.
Finally, here's an exercise. Find the differential equation that produces linear growth and rerun this example using ODE1 twice, once to compute the compound interest and once to compute the simple interest.