Application of Direction Field to Autonomous Equations

Working by hand and using the results of how derivatives affect the shape of a graph, one may often be able to construct global pictures of solutions to autonomous equations without much difficulty. This was the focus in class. Autonomous differential equations are also separable, so provided one can perform the integral the solution can also be found without much difficulty by hand.

In this Mathematica notebook, we will look at using Mathematica to construct direction fields, solve the DE, and get plots of the solutions to the DE.

Logistic Equation

The logistic equation can be used to model populations that have a carrying capacity (maximum population the resources can support).

Clear[y,y0,field,lines,sol,functionlist]

To begin, let's look at the logistic equation. This equation has the following form:

L[y_]:=  D[y,t]== r (1 - y/K) y

L[y[t]]

where r>0 and K>0.  Let's construct the direction field for this equation, for the specific values of r = 0.5, K= 7.

field=VectorPlot[Evaluate[{1,Last[L[y]]/.{K->7,r->0.5`}}],{t,0,10},{y,-5,15},ImageSize->{400,400},VectorScale->{Tiny,Small,None}]

Now we can solve the initial value problem using Mathematica:

sol = DSolve[{L[y[t]],y[0]==y0},y[t],t]

This assigns the solution to y(t):

y[t_] = y[t]/.First[sol];

Here we create a list of some of the solutions from the family of curves:

functionlist=y[t]/. {y0 -> Range[0.1,15,1] , K->7,r->0.5}

lines = Plot[Evaluate[functionlist],{t,0,10},
PlotRange -> {-3,15}, ImageSize -> {400,Automatic},AxesLabel-> {"t","General Solutions"}]

Show[lines,field,BaseStyle->16]

Notice how the solution (from viewing the direction field) for y(0) < 0 move away from y=0 exponentially. Above y=0 they also move away from y=0, and this leads us to conclude that y(t)=0 is an unstable equilibrium solution. From the graph we can also say that y(t)=K=7 is a stable equilibrium point.

Critical Threshold

The critical threshold equation can be used to model populations that have an critical threshold (minimum population below which the population  becomes extinct).

Clear[y,y0,field,lines,sol,functionlist]

The critical threshold differential equation has the following form:

L[y_]:=  D[y,t]==- r (1 - y/T) y

where r>0 and T>0.  Let's construct the direction field for this equation, for the specific values of r = 0.5, T= 3.

field=VectorPlot[Evaluate[{1,Last[L[y]]/.{T->3,r->0.5`}}],{t,0,10},{y,-5,15},ImageSize->{400,400},VectorScale->{Tiny,Small,None}]

Now we can solve the initial value problem using Mathematica:

sol = DSolve[{L[y[t]],y[0]==y0},y[t],t]

This assigns the solution to y(t):

y[t_] = y[t]/.First[sol];

Here we create a list of some of the solutions from the family of curves:

functionlist=y[t]/. {y0 -> {0.1,1,1.5,2.5,2.9,5} , T->3,r->0.5}

lines = Plot[Evaluate[functionlist],{t,0,10},
PlotRange -> {-3,15}, ImageSize -> {400,Automatic},AxesLabel-> {"t","General Solutions"}]

Notice how the solution for y(0)>T can approach infinity at a finite time. This is unexpected from looking at the original differential equation! This appearance of a singularity that is not evident from the differential equation is a possible feature of nonlinear differential equations.

Show[lines,field,BaseStyle->16]

Notice how the solution (from viewing the direction field) for y(0) < 0 move towards y=0 exponentially. Above y=0 they also move towards y=0, and this leads us to conclude that y(t)=0 is a stable equilibrium solution. From the graph we can also say that y(t)=T is an ustable equilibrium point.

Logistic Equation with a Critical Threshold

This was the situation we looked at in class.

Global picture of the solution

Clear[y,y0,field,lines,sol,functionlist]

Let's look more closely at the direction field for Logistic growth with a critical threshold. The Logistic growth with critical threshold is used to model populations, and is given by

L[y_]:=  D[y,t]== -r (1 - y/K) ( 1 - y/T) y /. {K->7 ,T-> 3 ,r->0.5}

for r>0 and 0<T<K. We have specificed values for the parameters.

L[y[t]]

field=VectorPlot[Evaluate[{1,Last[L[y]]}],{t,0,20},{y,-5,15},VectorScale->{Tiny,Small,None}]

Let's overlay this with a plot of the important equilibrium solutions:

lines = Plot[{7,3,0},{t,0,20},PlotStyle->Thick,
PlotRange -> {-2,10}, ImageSize -> {400,300},AxesLabel-> {"t",""},Ticks-> {Automatic,{{7,"K "},{3,"T "},{0,"0 "}}}]

Show[lines,field,BaseStyle->16]

We can see that for y(0) < T the function will approach zero as time approaches infinity; there is not a large enough population to begin with to sustain itself. For T<y(0) < K  The solution increases to K, and for y(0) > K the solution decreases to K. K is called the critical threshold of the system, and is the stable population that will result if the population does not die out. y(t) = K and y(t)  = 0 are stable equilibrium solutions, and y(t) = T is an unstable equilibrium solution.

Explicit Solution

This analysis is important to understand since Mathematica produces a very nasty looking solution to the initial value problem.

L[y_]:=  D[y,t]== -r (1 - y/K) ( 1 - y/T) y

sol = DSolve[{L[y[t]],y[0]==y0},y[t],t]

Implicit Solution

We can get a better looking solution is we separate and integrate--the problem is that going from the implicit solution to an explicit solution makes it look complicated. Notice I include a constant of integration:

Integrate[(1/(r y(1 - y/K) ( 1 - y/T)  )),y]==Integrate[-1,t]+c

Let's determine the value of the constant if we have the initial condition y(0)=y0 (notice I had to put in the absolute values for the integral):

Solve[((K T (((Log[Abs[y]])/(K T))+((Log[Abs[-K+y]])/(K (K-T)))+((Log[Abs[-T+y]])/(T (-K+T)))))/r)==c-t/.{t->0,y->y0},c]

So here is the implicit solution:

implicitsol=((K T (((Log[Abs[y]])/(K T))+((Log[Abs[-K+y]])/(K (K-T)))+((Log[Abs[-T+y]])/(T (-K+T)))))/r)==((K Log[Abs[y0]]-T Log[Abs[y0]]+T Log[Abs[-K+y0]]-K Log[Abs[-T+y0]])/(r (K-T)))-t

It's still pretty nasty! I can actually get a plot of these, but it is a bit a pain to make Mathematica show the functions. Also since the functions are implicit, we get some extra bits when we plot them.

ImplicitSolToPlot=implicitsol/.{K->7,T->3,r->0.5}

plot1=ContourPlot[Evaluate[ImplicitSolToPlot/.y0->2],{t,0,20},{y,0,10},PlotPoints->50,ContourStyle->{Thick,Orange}];
Show[plot1,lines,field,BaseStyle->20]

In the above sketch, we should only have the lower part that begins at y(0)=2. The other two curves are extraneous due to the implcit function we are plotting.

plot1=ContourPlot[Evaluate[ImplicitSolToPlot/.y0->3.1],{t,0,20},{y,0,10},PlotPoints->50,ContourStyle->{Thick,Orange}];
Show[plot1,lines,field,BaseStyle->20]

Again, notice the extraneous parts of the solution.