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.