## Section 2.6: Solving Equations on the Computer

Section 2.6: Solving Equations on the Computer

Solving Equations on the Computer

This course is about how to get as much out of your dynamical systems without having to solve any, or at least too many, differential equations.

Euler’s Method

Most of this is going to take the form of these geometrical methods that we’ve been describing - phase portraits, fixed points, stability analysis, and similar. We want to be able to get a computer to do the hard work for you. One of the simplest ways is so called Euler’s method, and it works particularly well for these single, first order differential equations. Let’s take a particular example we’ve seen before. The logistic equation:

P

where we have chosen the growth rate and the carrying capacity both to be 1.

Let’s say we want to know what’s going to happen if we start with a population of 0.2 (i.e. ). (let’s say that the units are in millions of... Hmmm... Rabbits). Let’s plot this point as the starting point of our trajectory:

P(0)=0.2

In[]:=

plfp=ListPlot[{{0,0.2}},PlotStylePointSize[0.02],PlotRange{{-0.1,5},{0,1.5}},AxesLabel{Style["t",18],Style["P(t)",18]}]

Out[]=

What would the rate of change of the population at this point by? Well we have an equation for that!

P

Clearly, if the population doesn’t change much, then the gradient isn’t going to change much. So let’s roll forward time with a rate of change of P of 0.16 for, let’s say 1 timestep. Where would this take us to? (Euler’s Method)

P=0.2+0.161=0.36

That’s now our new population. Let’s put that on our trajectory, and join the points:

In[]:=

plfp=Show[ListPlot[{{0,0.2},{1,0.36}},PlotStylePointSize[0.02],PlotRange{{-0.1,5},{0,1.5}},AxesLabel{Style["t",18],Style["P(t)",18]}],ListLinePlot[{{0,0.2},{1,0.36}}]]

Out[]=

Now we do the same thing again. Now we’re at , so what’s the rate of change now?

P=0.36

P

The gradient has increased a bit. Let’s use this gradient and roll forward again for one time step:

P=0.36+0.23041=0.5904

Which gives us:

In[]:=

plfp=Show[ListPlot[{{0,0.2},{1,0.36},{2,0.5904}},PlotStylePointSize[0.02],PlotRange{{-0.1,5},{0,1.5}},AxesLabel{Style["t",18],Style["P(t)",18]}],ListLinePlot[{{0,0.2},{1,0.36},{2,0.5904}}]]

Out[]=

We can continue with this... Using a computer program we can automate this process. Let’s do this for another few steps:

In[]:=

p={{0,0.2}};deltat=1For[i=1,i≤5,i++,pd=p[[-1,2]](1-p[[-1,2]]);p=Append[p,{i,p[[-1,2]]+pddeltat}]]plfp=Show[ListPlot[p,PlotStylePointSize[0.02],PlotRange{{-0.1,5.2},{0,1.5}},AxesLabel{Style["t",18],Style["P(t)",18]}],ListLinePlot[p]]

Out[]=

1

Out[]=

We see that it’s doing what we expect from a fixed point analysis. We know that is a fixed point, and we can see that we are gradually going towards this value.

P=1

The idea is to use the equation to approximate the gradient over the next short time step as the actual gradient at the start and moving in a straight line with that gradient through the interval and then repeat.

We could have taken smaller time steps which would be slightly more accurate. This is because the gradient of the smooth curve changes less over a smaller interval. Let’s go with time steps of 0.1:

In[]:=

p={{0,0.2}};deltat=0.1For[i=1,i≤50,i++,pd=p[[-1,2]](1-p[[-1,2]]);p=Append[p,{ideltat,p[[-1,2]]+pddeltat}]]plfp=Show[ListPlot[p,PlotStylePointSize[0.015],PlotRange{{-0.1,5.2},{0,1.5}},AxesLabel{Style["t",18],Style["P(t)",18]}],ListLinePlot[p]]

Out[]=

0.1

Out[]=

This is actually already very close to the exact solution and we didn’t have to do any integration! We can do the same thing for a range of different starting values of the population and we get this:

In[]:=

p0s={0.2,0.6,1.4,2};allps={};deltat=0.1;For[j=1,j≤Length[p0s],j++,p={{0,p0s[[j]]}};For[i=1,i≤50,i++,pd=p[[-1,2]](1-p[[-1,2]]);p=Append[p,{ideltat,p[[-1,2]]+pddeltat}]];allps=Append[allps,p]]plfp=Show[ListPlot[allps,PlotStylePointSize[0.015],PlotRange{{-0.1,5.2},{0,2.1}},AxesLabel{Style["t",18],Style["P(t)",18]}],ListLinePlot[allps,PlotRange{{-0.1,5.2},{0,2.1}}]]

I've done this in Mathematica, but now you try and do this using your Python skills.

Euler's method, is one of the simplest ways to explore these differential equations. One can do the same thing with higher order differential equations. In fact there are hundreds of books written on the subject of numerically solving differential equations, but we’re not going to explore that too much here.

Direction Fields

There is another computational method that we can use which is a so-called Direction Field. We’ll look at a non-autonomous differential equation, as it’s pretty simple, and you can see that it’s even simpler if the system is autonomous.

Let’s look at the equation:

So, what we are going to do is draw in a short line segment, let’s say of length 1 at this point, with gradient 1.

The idea here is that if you knew the true trajectory, as it passed through the point (0,1) it would look something like this. We can do the same thing at lots of other points in the plane. For instance, at the point (3,1), we would get a short line of gradient 4 which would look like:

Now let's do this all over the plane:

Let’s do it for more points, and with slightly shorter lines:

We can picture this like looking at a stream, where you have dropped a little bit of coloured dye into many points, and are watching the flow. You can imagine that you start at some point, move along that line, then find the nearest next line, and follow that, etc. This gives you a trajectory.

These lines are known as direction fields and are another way of figuring out the behaviour without actually having to solve the differential equation.

As I said, this is for a non-autonomous system. Had we gone with an autonomous system, let’s say

Our direction field would look like this:

Soon we will move on to section 3 and explore bifurcations, which is really about the way the fixed points can themselves behave. First of all, let’s have a quick recap of phase spaces and phase portraits.