0:09

Now we're going to move on to the second

Â part of our lecture on Introduction to Dynamical Systems.

Â And in this lecture, what we're going to focus on.

Â Is a, a classic method for numerically solving o,

Â o, ordinary differential equation systems known as Euler's method.

Â 0:26

And we're going to go through the general concepts that underlie in Euler's method.

Â We're going to run through an example with a

Â single variable and how you can use Euler's method

Â to solve ODEs and then we're going to talk a

Â little bit about how one can expand Euler's method.

Â To multivariable systems.

Â So the example with the single variable

Â is good to illustrate how Euler's method works.

Â Obviously, for biological simulations [UNKNOWN] you never

Â actually have one that's just of a

Â single variable, so we'll discuss how you

Â can easily expand Euler's method for multivariable systems.

Â Now I'm going to reintroduce the system that we looked at in the last

Â lecture, how one can describe biochemical sys signaling with a system of ODEs.

Â Right, last lecture we discussed the following three

Â component like network, where you have, protein A, a

Â protein B and a protein C that can all

Â be phosphorylated or unphosphorylated, and they regulate one another.

Â And what we discussed in the last lecture is how this

Â system can be described through this set of three ordinary differential equations.

Â 1:28

And this is these are relatively complicated nonlinear equations.

Â They cannot be solved analytically.

Â So they have to be solved numerically.

Â So what we want to address now in this lecture is

Â how can you simulate this numerically using the MATLAB Programming Environment.

Â So, so to demonstrate the general principals that underlie

Â the numerical solution of ODE's, we're going to discuss Euler's method.

Â This is named after a Swiss mathematician

Â Leonhard Euler, who lived in the 18th century.

Â 2:03

The way that you can define a differential equation

Â system is like this, you have a differential equation

Â of dx, dt for one variable x and this

Â differential equation is defined as some function of x.

Â And then the second thing you need

Â in a differential equation system is initial condition.

Â So, in order to solve this you need to you need to

Â know what describes the the kind derivative of x at all times.

Â And that's defined by this function f of x.

Â And then you also need to know what's the value of x at

Â some time, and that time is usually defined as time t equals zero.

Â So, at t equals zero.

Â The value of x is sum value of zero, which you know.

Â 2:43

And Euler had a, a very profound insights, which is that the derivative with

Â respect to time, dx dt, is roughly equal to the value x of the

Â at the next time point, T plus delta T, minus the value of x

Â at the current time point Divided by whatever that time difference is, delta t.

Â So the idea here is that if you're, if you, you

Â know, if you know x at time, some time t, say ten.

Â Then if, if you, if you want to, then you say, well what's the value of x going to

Â be at some future time, time ten plus one ten plus zero point one, ten point one.

Â Then the, the derivative.

Â Is approximately the difference in x between these two you know very closely

Â space time points divided by the time difference between those two time points.

Â And if we remember from calculus, this equation right here in the limit that

Â delta t goes to zero that's actually the formal definition of the derivative right.

Â This is how we define the derivative in calculus.

Â As delta t goes to zero this formula here is equal to the derivative.

Â So Euler's Euler's inside was okay, this is what happens when delta t in

Â the limit delta t goes to zero, so if delta t is not exactly

Â zero but its, its close to zero then this becomes it goes to approximation

Â so we can say the derivative is

Â approximately equal to this formula right here.

Â 4:07

And in our diff, ordinary differential equation system we know that we

Â can computer the derivative at any for any value of x, right?

Â That's our, our function f of x.

Â So what Euler did is instead of saying that

Â the derivative is equal to our function f of

Â x, is he took this approximation of the derivative

Â and said that was equal to f of x.

Â But why did he do that.

Â The reason he did that is, then he could take this equation here and very

Â simply he could rearrange that, so that on the left hand side you had something

Â that was unknown, something in the future and on the right hand side you have

Â three terms here that are all known, that are all at the current moment in time.

Â Right?

Â Alright, of t is something that is at, in the present.

Â So you know that.

Â F of x is something you can calculate if you know x.

Â And delta t is something you know.

Â So he's rearranged this so that something you don't

Â know in the future, is on the left hand side.

Â And something you do know, is three things that

Â you do know is on the right hand side.

Â 5:41

Then you can compute x at time three delta t, et cetera, etc etera.

Â And this is the basis of, of Euler's method.

Â It's simply this approximation right here.

Â The derivative is approximately equal to this formula.

Â 5:57

So now that we've discussed the principles underlining, underlying Euler's method,

Â we can show how you can implement Euler's method using MATLAB.

Â Let's assume that we have this ordinary differential equation.

Â dx dt equals a minus b times x.

Â And we know that at time t equals zero, x is equal to some value c.

Â And these are arbitrary numbers here for the values of a, b, and c.

Â If a, b, and c were 20, two, and five.

Â How could we write MATLAB code

Â to solve this ordinary differential equation numerically?

Â 6:33

A.

Â First we would define our constants, a, b, and c.

Â We would have to define a time step, dt equals zero point zero five.

Â We would have to define, okay, how long do you want to

Â numerically simulate this for, and that's

Â what this that the variable tlast represents.

Â 6:49

Now I've computed the number of iterations.

Â The number of iterations is, okay, how many time steps do we need to take.

Â So if you did it for 100 seconds, and the

Â time step was one second, then you would have 100 iterations.

Â If the time step was half of a second then you would have 200 iterations.

Â So that's all I'm computing here.

Â [INAUDIBLE] divided by delta t.

Â Then I set up this vector called x all which is going to be

Â set up to hold all the values of x for the different values of t.

Â And then this little for loop here is really the core of, of the program.

Â This is where we use Eulars' method to, to perform the numerical solution.

Â You set the value of x equal to c because.

Â The value, because the value c is equal to the initial condition of x.

Â And then we're going to perform iterations number of, of runs through the for loop.

Â So we say for i equals one, the iterations, then we perform

Â these calculations a certain number of times defined by the variable iterations.

Â And each time in the fore loop we just, we just do three things that are very simple.

Â We say, the current value of x all is equal to the current value of

Â x, or the ith value of x all is equal to the value of x.

Â So the first time through this is going to tell us the first value of x all.

Â The second time through it'll tells us, it'll

Â define a second value of x all, et cetera.

Â Then you compute the derivative, dxdt, as a minus b times x.

Â Very simple formula.

Â And then you compute the new value of x, is equal to the current value of x, plus

Â the derivative, times the time stack, which is what

Â we got from the last slide on, on Euler's method.

Â And this is all we do we just do this 100 times

Â or 2,000 times or however many times are defined by the variable iterations.

Â 8:27

Now we define a vector called time which from the MatLab commands

Â we discussed in the first couple lectures you can understand how this works.

Â This creates evenly spaced integers from zero to iterations minus one.

Â And then multiplies them by the time step to, account for the fact that,

Â your, you know, your delta t could be, could be much smaller than one.

Â And then you create a figure.

Â And then you plot time on the x axis, and x all on the y axis.

Â And this is the way to numerically solve differential equations using MatLab.

Â 9:10

Remember this is our, our system we're solving.

Â Dx dt equals a minus bx with a.

Â B, and the initial condition for, for x, to find as follows.

Â So, x starts at five.

Â And here, we're going to show how x evolves as

Â a function of time, as we perform these calculations.

Â 9:40

Then we compute the derivative of x with respect to time.

Â We have a formula for that, right?

Â Dx/ dt is equal to 20 minus two times five.

Â So it's equal to ten.

Â So what that tells us is that where we start here, x equals five.

Â 9:53

T6HESlope of this line, remember that the, you

Â know, that's the, what a derivative is, instantaneous slope.

Â Is equal to ten.

Â So then we can compute the new value of x, is equal

Â to the old value of x plus the derivative, times the time step.

Â So the new value of x is equal to five point five.

Â So now, where we are is up here.

Â 10:19

The next value of the derivative.

Â So now it's 20 minus two times five point five.

Â So now the the slope or the derivative is equal to nine.

Â So a got our slope here that's almost as steep as it was before, but not quite.

Â And then that gives us the new value of x, which is five point nine five.

Â And then we put that point there.

Â And then we can move on to the next derivative, the next value of x.

Â The next derivative, the next value of x, et cetera, et cetera.

Â 10:58

Next, we want to ask how well Euler's method actually did in this case.

Â So on the previous slide, we obtained a numerical solution to this ODE.

Â But is it the right solution?

Â That's what we want to ask now.

Â 11:19

And in this case, because we picked a, a simple ODE.

Â We can obtain an analytical solution for this ODE, which is defined as follows.

Â So we know what x of t is going to b, as a function

Â of, of time, and as a function of our parameters, a, b, and c.

Â For most ODEs of biological interest, we cannot obtain analytical solutions.

Â But because we started with something simple up here.

Â Just a minus bx.

Â We know what the analy, analytical solution to this is in this case.

Â 11:47

So, what we're going to plot here is x versus time.

Â And this blue line here, is the solution defined by the equation listed up here.

Â And then what we can do, is we can plot, the, the Euler's method solution.

Â For a small delta t.

Â And we can see that these red dots follow along the blue line almost perfectly.

Â 12:08

If we increase delta t, we have like a

Â medium delta t, which are the green dots here.

Â We see that it doesn't really match so well.

Â There's, there are errors here, but it's decent.

Â It's, it's pretty close, and then if we make delta t

Â a lot larger, which are the, the pink or magenta dots.

Â We see.

Â This particular, an this particular numerical solution is way off.

Â It really overshoots the value here and then eventually comes back

Â but in the early time points it's not accurate at all.

Â 12:52

One of the things that makes this appealing is

Â that extending it to systems of ODEs is very straightforward.

Â So if you define three variables, d little x d t This is a function of x, y and z.

Â d little y dt is of is some other function g of x, y

Â and z and a dz dt is equal to h of x, y, z.

Â Its really simple to use Euler's method to solve this system of ODEs.

Â You just have to define a vector v.

Â Which consist of our, your three variables, x, y, and z.

Â 13:25

And now you, you can just extend Euler's method to, to vector form

Â where you say the new vector times delta t is equal to the vector

Â times zero plus delta t times the vector of your three derivatives, where the

Â first element will be f of x, y, and z, which defines the xdt.

Â The second will be G.

Â The third one will be H.

Â 13:56

We saw, an example of an ina,ina, interactive solution on a previous slide.

Â But in other cases, solutions can become very, very unstable.

Â So if you're doing this in MatLab, one of two things, usually occurs.

Â One is that you can get variables to just be, be

Â totally absurd values, like 4.3 times ten raised to the 78th power.

Â So, if you see, if you, if you do an Euler's method implementation of a

Â mathmatical model And you see these absurd values

Â like, you know, ten to the 78th power.

Â 14:28

Sometimes that means your, your time step is too large.

Â Other times that means you, you made some other error.

Â For instance if you were supposed to add two variables together, and instead

Â you multiplied them together, that can

Â also make your numerical solution very unstable.

Â And it can give you these, you know, these very absurd values for your variables.

Â And then the other thing that can occur is that you know

Â in a lot of biological models, there's certain values that can't be negative.

Â For instance, if you're simulating

Â concentrations of species, well the concentration

Â can go down to some value very, very close to zero.

Â But it can't become negative and so if you're looking at your variables, they

Â have to be positive or it they have to be non-negative, and they become negative.

Â A lot of times they means that you are, your time staff delta t is too large.

Â So that's usually one of the first things you want to try, in this case.

Â And can happen, pretty commonly with Euler's method.

Â 15:23

Well, it's, important to note that, you know,

Â as we discussed, Euler lived in the 18th century.

Â So it's been a long time since, since he developed this method.

Â And over those intervening, you know, two plus centuries, applied mathematicians

Â have devised many more stable and more accurate methods for solving [INAUDIBLE].

Â ODEs.

Â And these are the algorithms that are implemented in MatLabs built-in solvers.

Â So, we're going to discuss this later, but MATLAB has these functions

Â that are called things like ode23, ode15s, there's a whole series of them.

Â 16:02

In this course we are not going to go into

Â details underlying these more complex numerical methods but I, I

Â thought it would be worthwhile just to illustrate a couple

Â of import points about these numerical methods that are used.

Â 16:14

One classic one is called the Runge-Kutta method.

Â and the idea here is that if you are trying to get from point one,

Â all the way to, to point four if you just approximates the derivative, at point one

Â which is what Euler's method does then you are going to go from point one, you know

Â over here, you are not going to get close to the correct value of point four.

Â But the idea is not to compute, you

Â don't want to only compute the derivative of point 0ne.

Â You also want to compute the derivative at some intervening points.

Â And so if you also computer a derivative at point two and point three, and then you

Â use a weighted average of these derivatives, a

Â derivative at one, derivative at two, derivative at three.

Â Then by computing this weighted average you can get from, from location,

Â from point one all the way to point four in one step.

Â 17:03

So the basis of Runge-Kutta is not to just

Â compute a derivative at current time point where you

Â are but to compute multiple derivatives at multiple locations,

Â take a weighted average, and then undertake your time stuff.

Â And then another Strategy that's used in these

Â algorithms is to take a variable type step.

Â So what we discussed with Euler's method and,

Â and with Runge-Kutta is you're going to have a same,

Â you know, the same delta t each time you move on to the next moment in time.

Â But intuitively, we can look at a function that looks

Â like this, which is changing very rapidly at the beginning.

Â But then it gets very flat at the end.

Â 17:42

And intuitively, it, it makes sense that when the function is changing very

Â rapidly, when the variable is changing rapidly, you can't take a very large

Â time step, because if you try to approximate the derivative here, then you

Â could miss this local maximum if you take a time step that's too big.

Â So here you're going to need a small delta t.

Â But over here where it's flat and it's not changing very rapidly

Â at all, here you might be able to take a big time step.

Â And this is a way you can really speed up

Â numerical,um, solutions of ODE's, is through these variable time step methods

Â which are based on this, this, underlying principle that when things

Â are changing rapidly you have to take a small time step.

Â When things are changing slowly it's okay to take a bigger time step.

Â 18:37

We showed the, the actual MATLAB Code for Euler's method in a few slides ago.

Â But let's review what the structure of this is.

Â First we define the constants.

Â The first things we set in this code were

Â a, b, and c were equal to uh,particular values.

Â 18:54

Then we set the time step and

Â the simulation time, and the number of iterations.

Â We said time step is zero point zero five.

Â Do it for two seconds.

Â So you have to define a time step,

Â you have to define how long you're going to go.

Â 19:07

Then we set the initial conditions our value of x

Â was equal to the, the variable that we defined as

Â c, and then the core of our program was a

Â for loop, and the for loop simulated the evolution of time.

Â And at each time step what you do is.

Â Write the output if, if needed that's we did we

Â saw we saved our current value of x in a

Â variable that we called x all, then you compute dx,

Â dt and in this case I have made it a capital.

Â So what I, what I mean here is that capital x

Â here refers to the vector of all the state variables The sol,

Â the strategy that we use to solve ODE for one variable, can

Â easily be expanded to co compute the numerical solutions in multiple variables.

Â And by, when I say, compute dx, dt, I mean compute all the derivatives.

Â And then, what you want to do is you want to compute x.

Â Again, all your variables, at the next time stop.

Â And then finally, when your for loop is

Â over, you plot the results and you output them.

Â [SOUND]

Â [BLANK_AUDIO]

Â Now to summarize.

Â 20:55

Or, it doesn't necessarily, or sometimes you just have

Â to make your delta t very, very, very small.

Â And then it becomes impractical to solve your, your system

Â with with a very, very small delta t like that.

Â So, because of that more complex numerical algorithms

Â are prefered for most models of, of biological interest.

Â I fi, I think it's helpful to go

Â through the simple algorithm and the classic algorithm but

Â for most of the solutions that we're going to

Â solve later in this class and most of the.

Â