Lab 1:  First Order Differential Equations

Part 1 - Analyzing equilibrium behavior

For each of the following differential equations: graph the slope field over a reasonable domain, find all equilibrium values, explain what happens to solutions that start between these values (i.e., do they head towards a particular equilibrium values, away from an equilibrium value, both?), and graph a solution or two on the slope field that demonstrates each of these different qualitative behaviors.

So, for example, if your equation were:  y' = 0.5  y(4 - y), your answer would look something like:

<< Graphics`

Clear[y, y0, t, slopeField, solutionCurves]

slopeField = 
      PlotVectorField[{1, 0.5  y (4 - y)}, {x, -5, 5}, {y, -5, 
          6}, Axes -> True]

Notice, there are 2 equilibrium values:  0 and 4.  For simpler cases like this, you can probably tell what is going on by just looking at the slope field, but if the equation is an autonomous differential equation, you could also graph y vs. y' on the same axis (so here, you are graphing 0.5  y(4 - y)):

Plot[0.5  y (4 - y), {y, -5, 5}, 
      AxesLabel -> {y, y'}]

The equilibrium values are where this graph crosses the horizontal axis.  Notice, that if y < 0, the slopes are all negative (making y head towards -∞).  If 0 < y < 4, the slope is positive, so y must head upwards (towards 4).  And if y > 4, the slope is negative, so y must head back down (towards 4 again).

Thus, all solutions starting below 0 head toward -∞, all solutions starting between 0 and 4 head away from 0 and towards 4, and all solutions that start above 4 head towards ∞.  So, it might make sense to choose 3 solutions to demonstrate this:  y(0) = -2,  y(0) = 1,  and y(0) = 5

solution = 
      DSolve[{y'[t] == 0.5  y[t] (4 - y[t]), y[0] == y0}, y[t], 
        t]

ySol[t_] = y[t] /. solution

y0 = {-2, 1, 5}

solutionCurves = 
      Plot[Evaluate[ySol[t]], {t, -5, 5}, 
        PlotStyle -> Blue]

(Notice that you got a few solutions with a vertical asymptote.  Look back up at the slope field and see if this makes any sense.)  Now, combine these with the slope field:

Show[slopeField, solutionCurves, 
      PlotRange -> {-5, 6}]

For each of the following differential equations, perform a similar type of analysis, showing your work in Mathematica, along with appropriate explanations to go with it.  Warning:  For some of these, you may want to use NDSolve to plot the solutions instead of DSolve.

y' = 4  y - (y^3)

(hint: only plot and solve for t ≥ 0)

y' = y(1 - ln y) (y - 4)

(Hint: use NDSolve on this one.  Also, remember that Mathematica uses Log for natural logarithm.)

Part 2 - Sensitivity to initial conditions and NDSolve

A differential equation is said to be insensitive to initial conditions (or "stable") if a small change in the initial condition causes a (relatively) small change in the behavior of the solution (over time).  It is sensitive to initial conditions if a small change in the initial conditions causes a (relatively) significant change in the behavior of the solution.  So, as a simplistic example, consider y' = 0.4  y.  This has a solution of y = k (e^(0.4  t)), and if we graph this for several closely spaced values of k, we get:

k = Range[1, 1.1, 0.01]

Plot[Evaluate[k (E^(0.4  t))], {t, 0, 5}]

As you can see, these solutions draw further and further apart over time.  If you graph over a larger time period, this becomes even more evident:

Plot[Evaluate[k (E^(0.4  t))], {t, 99, 100}]

(Notice, these solutions that started out only covering a range of 0.1, now spread over a range of about 0.25*(10^17) = 2.5*(10^16).)  The opposite case can be shown by the equation y' = -0.4 y:

k = Range[-10, 10, 1]

Plot[Evaluate[k (E^(-0.4 t))], {t, 0, 50}, 
      PlotRange -> All]

Here, all of the solutions converge towards 0, so this equation clearly is pretty insensitive to initial conditions (stable).  

A useful fact (which I'm not even going to try to prove here) is that if you have a first order differential equation of the form y' = f(x, y), then the equation is stable if ∂(f/∂)y < 0 and unstable if ∂(f/∂)y > 0.  (Things become a bit more complicated if ∂(f/∂)yis sometimes greater and sometimes less than 0.)  Thus, in the first equation, ∂(f/∂)y = 0.4, which means it is unstable (as we saw).

Now, this can cause a problem if you are using a numerical method to approximate the solution to a differential equation (i.e., NDSolve).  If the equation is stable, then life is good.  Your approximation should be pretty safe.  However, if the equation is unstable, then each small error in the approximation has the effect of kicking the solution over onto a slightly different solution curve; since the equation is sensitive to initial conditions, that means that when the error moves it over to a new curve, that curve diverges from the original solution curve more and more over time, potentially greatly increasing your total error.  By default, Mathematica tries to make the error for approximating differential equations ~(10^(-6)) for each step of the approximation.  If the solutions are stable, then the overall error should remain pretty small.  However, if the solutions are unstable, then those small errors can cause a much larger error, especially over longer time intervals.  So, in our example, if we use NDSolve to approximate the solution and evaluate that at t = 1000, compared to the actual solution:

approxSolution = y[t] /. 
          NDSolve[{y'[t] == 0.4  y[t], y[0] == 1}, 
            y[t], {t, 0, 1000}] /. t -> 1000

exactSolution = (E^(0.4  t)) /. 
        t -> 1000

approxSolution - exactSolution

I think most people will agree that this is a pretty significant difference (though notice that it will actually be quite accurate for small t values).  You can correct for this by increasing the accuracy goals for NDSolve.  The details of this are a bit complicated, but generally this is pretty easy to do:

approxSolution = y[t] /. 
          NDSolve[{y'[t] == 0.4  y[t], y[0] == 1}, 
            y[t], {t, 0, 1000}, WorkingPrecision -> 16] /. 
        t -> 1000

approxSolution - exactSolution

This still isn't great, but it IS better (by several orders of magnitude).  For our purposes, WorkingPrecision should generally be about 10 greater than the number of digits you plan to trust (the default value for WorkingPrecision is about 16).  The error message here is complaining that the 0.4 in the equation doesn't even come close to the requested WorkingPrecision.  You can fix that by replacing all constants in the equation with exact values:

approxSolution = y[t] /. 
          NDSolve[{y'[t] == (4/10) y[t], y[0] == 1}, 
            y[t], {t, 0, 1000}, WorkingPrecision -> 25] /. 
        t -> 1000

approxSolution - exactSolution

(Notice that this also makes your result noticeably more accurate as well.)  Don't go too crazy with this, however.  If you decide to shoot for 25 digits accuracy, for example, you run into the following problem:

approxSolution = y[t] /. 
          NDSolve[{y'[t] == (4/10) y[t], y[0] == 1}, 
            y[t], {t, 0, 1000}, WorkingPrecision -> 75] /. 
        t -> 1000

The error message here is warning you that it isn't able to maintain the requested accuracy over the entire interval.  This means the solution is not trustworthy (as you can see by the pretty wildly inaccurate solution it gives).  There are ways to deal with this, but I don't plan on covering this (look it up in the Help system if you are curious).

Problem

Consider the differential equation y' = y - 3  (e^(-2 t)).  Answer the following questions about it:

a)  Is this equation sensitive to initial conditions or not?  Why/why not?

b)  Use DSolve to plot solutions with initial conditions near x = 1 (say, from 0.9 to 1.1, perhaps with a stepsize of 0.01).  Does this support or contradict your answer in part a?  (Hint:  Plot all solutions for this problem out to at least t = 15.)

c)  Use NDSolve to find an approximate solution with initial condition y(0) = 1 and plot that on a graph with the solution obtained with DSolve (same initial condition).  Explain why NDSolve gives a different answer (hint: "because it's an approximation" won't get you any credit here).  Is this difference "significant"?  Why/why not?

d)  Increase the WorkingPrecision for NDSolve gradually (remember, the default value starts at 16) until it agrees with DSolve.  Using that WorkingPrecision, increase the domain for t to be twice as long and plot it against DSolve again.  Find a new WorkingPrecision that will fix this one.

e)  Plot NDSolve with WorkingPrecision of 25 (0 ≤ t ≤ 30).  What is going on with this?  (It should look qualitatively different than most of the solutions you plotted in part d.)

Part 3 - An application to population growth

As we have discussed in class, there are several varieties of equations we can use to model population growth.  In this part, we will investigate some further extensions to this.  In this situation, we will consider the population of fish in a pond.  We can use a logistic growth equation:

	Q' = r Q (1 - (Q/K))

Here:
	Q(t) = Quantity of fish at time t
	r = growth rate
	K = carrying capacity of the lake (or "Karrying Kapacity")
	
However, we need to also add in a term to take care of any fish "harvesting" that might occur (people do get hungry...).  We will use H(t) to represent the rate at which the fish are harvested over time:

	Q' = r Q (1 - (Q/K)) - H(t)
	
In this case, assume that r = 4,  K = 100,  and T = 10.  For each of the following types of harvesting, answer the following questions (hint: make liberal use of copy and paste on this...):

a)  If the eqaution is autonomous, find the equilibrium points and describe the behavior between them.  If it isn't autonomous, graph the slope field to pick some good initial values.  (Alternately, you could just plug in a wide range of values, covering the region of interest, perhaps 10 ≤ y ≤ 150 or so.  Also, generally you only need to plot solutions for about 0 ≤ t ≤ 4.  If NDSolve give you a lot of errors, you could either use fewer initial points, use a smaller domain, or try a higher WorkingPrecision.)

b)  Pick several initial values to demonstrate the different possible behaviors of the population (you only need to consider cases where y ≥ 0).

c)  Clearly explain (using actual words and sentences) what will happen to the fish population in this lake over time, given the different possible types of initial conditions.  Explain what effect the harvesting term has (compared to no harvesting).

Constant harvesting:  H(t) = 90

Constant harvesting:  H(t) = 110

"Constant effort" harvesting:  H(t) = 3  Q

Periodic harvesting:  H(t) = 51 (1 - sin π t)

Periodic harvesting:  H(t) = 101 (1 - sin π t)