## Section 2.6: Flows on the line - Solving equations on the computer

Section 2.6: Flows on the line - Solving equations on the computer

As stated before, this course is really about how to get as much out of your dynamical systems without actually having to solve any, or at least too many, differential equations.

Much of this is going to take the form of these geometrical methods that we’ve been describing - phase portraits, fixed points, stability analysis, and the like. One might also 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. One that 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. ie . (let’s also 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? Well, it would be:

P=0.2+0.16x1=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 P=0.36, so what’s the rate of change now?

P

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

P=0.36+0.2304x1=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 keep on with this...in fact if we know how to write 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 P=1 is a fixed point, and we can see that we are gradually going towards this value.

The idea of all of this is to use the equation to predict the gradient over the next short time step as being the gradient at the start and moving across the interval on a straight line of that gradient, 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 a single 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}}]]

Out[]=

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

This method, called 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 actually 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 like like:

Now let’s do this all over the plane:

Let’s do it for more points still, 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 will give 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.