Calculus: The Language of Change

Exercise Set 9.5
Graphing Differential Equations

by John Robeson

Exercise 1) Exponential Decay Without Formulas

A fundamental "decay law" says that the rate of decrease of a quantity is proportional to the amount that is left.  Radioactive substances have this property.  You can sketch the amount of such a substance as a function of time without any formula.

(a)  Given that a quantity q satisfies q[0] = 4 and (dq/dt) = -(q/2) , sketch the graph of q vs. t.

(b)  What happens to the quantity as it gets close to q = 0?

(c)  Can you find an analytical solution to the conditions q[0] = 4 and  (dq/dt) = -(q/2) ?  (See Chapter 8.  This is not needed for the sketch.)

(a)

Given that a quantity q satisfies q[0] = 4 and (dq/dt) = -(q/2) , sketch the graph of q vs. t.

Solution to (a)

We know that the graph of q vs. t starts at (0, 4) since q[0] = 4.  We also know that the slope is negative for all positive values of q, so the graph curves downward.  Finally, the second derivative is the derivative of (dq/dt), which is


	(((d^2) q)/(dt^2)) =(d/dt)((dq/dt)) = (d/dt)(((-q)/2)) =-(1/2) (((-q)/2)) = (q/4)

So the second derivative is positive (a smile) for all positive values of q.  Evaluate the cell below for the graph.

Clear[q, t]q  =  4 E^(-t)Plot[q, {t, 0, 10}, PlotRange -> {0, 5}, 
      AxesLabel -> {"t", "q"}]

(b)

What happens to the quantity as it gets close to q = 0?

Solution to (b)

The quantity is still decreasing, but the slope is getting very small because q is very small.  So the graph is "flattening."

(c)

Find an analytical solution to the conditions q[0] = 4 and

(dq/dt) =  -(q/2).  (See Chapter 8.  This is not needed for the sketch.)

Solution to (c)

From Section 8.2 (Exercise Set 8.2 number 1(d)) we know that our solution is q = 4 (e^(k t)) where k = -(1/2)

We can check this solution.

	q[0] = 4 (e^0)= 4
	
	(dq/dt) = 4 k  (e^(k t))	 = -(q/2)

Exercise 2) Sketch Solutions

Sketch the solutions of

		(dy/dt) = y(4 - y)
		
and		(dy/dt) = 10y(4 - y)
that begin with

(a)	y[0] = 1

(b)	y[0] = 2

(c)	y[0] = 6

(d)	y[0] = 4

(e)	y[0] = 0

Solution for (dy/dt) = y(4 - y)

First let's compute the second derivative.

		(((d^2) y)/(dt^2))= ((dy (4 - y))/dt) + ((y (- dy))/dt)
		
		
		= ((dy (4 -  2  y))/dt) = y (4-y) (4-2 y)

(a)	For y[0] = 1 we know the graph starts at (0, 1) and the inital slope is 1(4 - 1) = 3 which means that y is increasing.  As y approaches 4, the slope is still positive but becomes very small, so our graph flattens out as t gets as large as we want.  The second derivative is initially 1(4 - 1)(4 - 2) = 6.  It remains positive (a smile) until y = 2.  For all y between 2 and 4 the second derivative is negative (frown).

\(Needs["Graphics`Colors`"]; 

    eulerApprox[f_, {t0_, x0_}, dt_, final_]\  := 

      Module[{steps, G}, 
steps\  = \ Floor[final/dt]; 

        G[{tt_, xx_}]\  := \ {tt + dt, xx\  + \ f[{tt, xx}]\ dt}; 

        NestList[\ G\ , \ {t0, x0}\ , \ steps\ ]]\)\(Clear[t0, x0, G, f]; 
G[{t_, x_}] := {t + dt, x + f[t, x]\ dt}; 

    Clear[f]\)\(f[t_, y_] := y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 1.0; 
dt = 0.1; 

    finalTime = 2; 
solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 
    
plot1 = 
      ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; \)

(b)	For y[0] = 2 we have the same situation as (a) except that our graph starts at (0, 2) and the second derivative is initially zero, but since y is increasing, the second derivative quickly becomes negative (frown).

Clear[f]\(f[t_, y_] := y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 2.0; 
dt = 0.1; 

    finalTime = 2; 
solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 
    
plot2 = 
      ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; \)

(c)	For y[0] = 6 we know the graph starts at (0, 6) and the inital slope is 6(4 - 6) = -12 which means that y is decreasing.  As y approaches 4, the slope is still negative but becomes very small, so our graph flattens out as t gets as large as we want.  The second derivative is initally 6(4 - 6)(4 - 12) = -96 and remains negative (frown) for all y between 6 and 4.

Clear[f]\(f[t_, y_] := y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 6.0; 
dt = 0.1; 

    finalTime = 2; 
solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 
    
plot3 = 
      ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; \)

(d)	For y[0] = 4 we know the graph starts at (0, 4) and the inital slope is 4(4 - 4) = 0 which means that y starts at 4.  The second derivative is initially also 0 which means that the first derivative does not change so our graph is a horizontal line, y = 4.

Clear[f]\(f[t_, y_] := y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 4.0; 
dt = 0.1; 

    finalTime = 2; 
solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 
    
plot4 = 
      ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; \)

(e)	For y[0] = 0 we know the graph starts at (0, 0) and the inital slope is 0(4 - 0) = 0 which means that y is initially flat.  The second derivative is initially also 0 which means that the first derivative does not change so our graph is a horizontal line, y = 0.

Clear[f]\(f[t_, y_] := y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 0.0; 
dt = 0.1; 

    finalTime = 2; 
solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 
    
plot5 = 
      ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; \)

Now all 5 graphs together.

\(\(\(Show[plot1, plot2, plot3, plot4, plot5]; \)
\)\)

Solution for (dy/dt) = 10 y(4 - y)

First let's compute the second derivative.

		(((d^2) y)/(dt^2))= 10[((dy (4 - y))/dt) + ((y (- dy))/dt)]
		
		
		= 10[((dy (4 -  2  y))/dt)] = 10 y (4-y) (4-2 y)

(a)	For y[0] = 1 we know the graph starts at (0, 1) and the inital slope is 10(4 - 1) = 30 which means that y is increasing even faster than in the first part.  As y approaches 4, the slope is still positive but becomes very small, so our graph flattens out as t gets as large as we want.  The second derivative is initially 10(4 - 1)(4 - 2) = 60.  It remains positive (a smile) until y = 2.  For all y between 2 and 4 the second derivative is negative (frown).

Clear[f]\(f[t_, y_] := 10\ y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 1.0; 
dt = 0.01; 
    
finalTime = 2; 

    solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 

    plot6 = ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Blue}]; 
    \)

(b)	For y[0] = 2 we have the same situation as (a) except that our graph starts at (0, 2) and the second derivative is initially zero, but since y is increasing, the second derivative quickly becomes negative (frown).

Clear[f]\(f[t_, y_] := 10\ y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 2.0; 
dt = 0.01; 
    
finalTime = 2; 

    solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 

    plot7 = ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Blue}]; 
    \)

(c)	For y[0] = 6 we know the graph starts at (0, 6) and the inital slope is 60(4 - 6) = -120 which means that y is decreasing.  As y approaches 4, the slope is still negative but becomes very small, so our graph flattens out as t gets as large as we want.  The second derivative is initally 60(4 - 6)(4 - 12) = -960 and remains negative (frown) for all y between 6 and 4.

Clear[f]\(f[t_, y_] := 10\ y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 6.0; 
dt = 0.01; 
    
finalTime = 2; 

    solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 

    plot8 = ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Blue}]; 
    
\)

(d)	For y[0] = 4 we know the graph starts at (0, 4) and the inital slope is 40(4 - 4) = 0 which means that y starts at 4.  The second derivative is initially also 0 which means that the first derivative does not change so our graph is a horizontal line, y = 4.

Clear[f]\(f[t_, y_] := 10\ y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 4.0; 
dt = 0.01; 
    
finalTime = 2; 

    solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 

    plot9 = ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Blue}]; 
    
\)

(e)	For y[0] = 0 we know the graph starts at (0, 0) and the inital slope is 0(4 - 0) = 0 which means that y is initially flat.  The second derivative is initially also 0 which means that the first derivative does not change so our graph is a horizontal line, y = 0.

Clear[f]\(f[t_, y_] := 10\ y\ \((4\  - \ y)\); 
t0 = 0; 
y0 = 0.0; 
dt = 0.01; 
    
finalTime = 2; 

    solnList = NestList[G, {t0, y0}, Floor[finalTime\/dt]]; 

    plot10 = ListPlot[solnList, PlotJoined -> True, 
        PlotStyle -> {Blue}]; \)

So the graphs of the second part are similar to those of the first part except that they flatten out faster.  Now all 10 graphs.

\(\(\(Show[plot1, plot2, plot3, plot4, plot5, plot6, plot7, plot8, plot9, 
      plot10]; \)

\)\)

Exercise 3) Euler's Approximation

Suppose y[0] = (1/4) and (dy/dt) = (Sqrt[y]).  Why is the tangent line to the exact solution y = y[t] below the curve at t = 0?  If Euler's approximate solution (Subscript[y, approx]) [δt] is low at the first step, why is the slope f[ (Subscript[y, approx]) [δt]] below the slope of y[t] at t = δt ?

Solution

First let's compute the second derivative.

	(((d^2) y)/(dt^2)) = (1/2) (1/(Sqrt[y])) (dy/dt) = (1/2) (1/(Sqrt[y])) (Sqrt[y]) = (1/2)
	
Since the second derivative is always positive, the function is smile shaped.  We know y[0] = (1/4) and (dy/dt)[0] = (Sqrt[0])  =   0.  Our function is smile shaped and the tangent at 0 is flat, so it is beneath the function.  In fact, since our function is always smile shaped, the tangent is always below the function.

Exercise 4) Solutions to the S-I-S Equations

An S-I-S disease is one that does not confer immunity; you are either susceptible or infectious.  If s is the fraction of the population that is susceptible, the spread of the disease in time t is given by
		(ds/dt) = b (1 - s) (1 - c s)
		
where b is the reciprocal of the infectious period and c is the contact number.  Assume that c ≻ 1 (so that each infectious person contacts more than one other person) and sketch the graphs of the solutions for various initial conditions between 0 and 1.

Solution

(ds/dt) = 0 when s = 1 and when s = (1/c).  For 0 ≺ s ≺ (1/c), the first derivative is positive and for  (1/c) ≺ s ≺ 1, it is negative.  The 

second derivative is (b^2)(s - 1)(1 - c s)((1 + c) - 2 c s ) so it is negative (frown) for 0 ≺ s ≺ (1/c), positive (smile) for (1/c)≺ s ≺ ((c + 1)/(2 c)) and negative for ((c + 1)/(2 c)) ≺ s ≺ 1.

Evaluate the cells below for graphs with b = 1/4 and c = 3.

Clear[f, b, c]\(b\  = \ 1/4; 
c\  = \ 3; 

    f[t_, s_] := b\ \((1 - s)\)\ \((1 - c\ s)\); 
t0 = 0; 
s0 = 0.1; 

    dt = 0.1; 
finalTime = 20; 

    solnList = NestList[G, {t0, s0}, Floor[finalTime\/dt]]; 

    plots1 = ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; 
    \)

Clear[f, b, c]\(b\  = \ 1/4; 
c\  = \ 3; 

    f[t_, s_] := b\ \((1 - s)\)\ \((1 - c\ s)\); 
t0 = 0; 
s0 = 0.3; 

    dt = 0.1; 
finalTime = 20; 

    solnList = NestList[G, {t0, s0}, Floor[finalTime\/dt]]; 

    plots2 = ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; 
    \)

Clear[f, b, c]\(b\  = \ 1/4; 
c\  = \ 3; 

    f[t_, s_] := b\ \((1 - s)\)\ \((1 - c\ s)\); 
t0 = 0; 
s0 = 0.7; 

    dt = 0.1; 
finalTime = 20; 

    solnList = NestList[G, {t0, s0}, Floor[finalTime\/dt]]; 

    plots3 = ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; 
    \)

Clear[f, b, c]\(b\  = \ 1/4; 
c\  = \ 3; 

    f[t_, s_] := b\ \((1 - s)\)\ \((1 - c\ s)\); 
t0 = 0; 
s0 = 0.9; 

    dt = 0.1; 
finalTime = 20; 

    solnList = NestList[G, {t0, s0}, Floor[finalTime\/dt]]; 

    plots4 = ListPlot[solnList, PlotJoined -> True, PlotStyle -> {Red}]; 
    \)

\(\(Show[plots1, plots2, plots3, plots4]
\)\)

Notice that all graphs approach a final value of (1/c).