Chapter 14

Needs["Graphics`MultipleListPlot`"]
Needs["Graphics`Legend`"]
Needs["Graphics`PlotField`"]

Question 1

k = ((Log[6] - Log[5.8])/10)

estp = 5.8  (E^(k (t - 1701)))

pointsestp = Table[{t, estp}, {t, 1701, 1791, 10}]

pointsactp = {{1701, 5.8}, {1711, 6.0}, {1721, 6.0}, {1731, 
          6.1}, {1741, 6.2}, {1751, 6.5}, {1761, 6.7}, {1771, 7.2}, {1781, 
          7.5}, {1791, 8.3}}

\(\(MultipleListPlot[pointsactp, pointsestp, PlotJoined -> True, 
        PlotRange -> {5.5, 8.5}, 
        AxesLabel -> {"year", "population (million)"}, 
        PlotLabel -> "Population of England and Wales, 1701-1791

", 
        PlotLegend -> {"actual", "estimated"}, 

        LegendPosition -> {\(-0.35\), \(-1.2\)}, 

        LegendOrientation -> Horizontal, 

        LegendSize -> {1,  .5}];\)\)

Question 2

We have the following information for each country after the epidemic:
	Country A: b = 3%, d = 2% with (b-d) = 1%	λ = 0.01
	Country B: b = 5%, d = 3% with (b-d) = 2%	λ = 0.02

Solve[2  p0 (E^(0.01  t)) == p0 (E^(0.02  t)), t]

Question 3

We have
	p(t) = p0 (E^([(b - d) - m] t))

FullSimplify[
      Solve[p0 (E^(((b - d) - m) t1)) == (1/2) 
            p0 (E^(((b - d) - m) t0)), t1]]

Apart[((Log[2] - Log[(E^((b - d - m) t0))])/(-b + d + m))]

But
	\(TraditionalForm
      \`Ln(E\^\(\((b - d - m)\) t0\)) = \((b - d - m)\) t0\)
Hence,
	t1 = ((-Log[2])/(b - d - m)) + t0t1 - t0 = ((-Log[2])/(b - d - m))

To check, using the information in the previous question, let b = 5%, d = 3% and m = 3%, then b-d-m = -1%

-Log[2]/(-0.01)

Question 4

Given
	p(t) = p0 (E^(k t))

(i) Population trebbling in size

Solve[p0 (E^(0.03  t1)) == 3  p0 (E^(0.03 t0)), t1]

Simplify[Solve[t1 == 33.3333 (0.03 t0 + Log[3]), t1]]

i.e., t1-t0 = 36.6204

(ii) General result

Clear[p0, λ, k]

Solve[p0 (E^(k t1)) == λ p0 (E^(k t0)), t1]

Which can be expressed
	t1 - t0 = ((Log(λ))/k)

Question 5

To derive the Taylor expansion about the point a/b, we use the Series command

Series[p (a - b p), {p, a/b, 1}]

Neglecting the error term, then
	\(TraditionalForm\`\(p\& . \) = \(-\(a(p - a\/b)\)\)\)

DSolve[{p'[t] == -a (p[t] - (a/b)), p[0] == p0}, p[t], 
      t]

Which can be expressed
	\(TraditionalForm
      \`p(t)\  = \ a\/b + \((p0 - a\/b)\) E\^\(\(-a\)\ t\)\)

Mathematica cannot take the limit of (E^(-a t)) since it does not know whether a is positive or negative. Assuming a is positive, and given

Limit[(E^(-t)), t -> ∞]

then p(t) -> a/b, which implies that equilibrium is never achieved in a finite time period.

Question 6

(i)

The fixed points are found from

Solve[p (a + c p) == 0, p]

The turning point is

Solve[D[p (a + c p), p] == 0, p]

which is negative since a > 0 and c > 0. Furthermore,

D[D[p (a + c p), p], p]

which is positive, hence the turning point is a minimum.

(ii)

These observations can be verified by assuming for example a = 0.2 and c = 0.004.

Clear[a, c, p]

example = p (a + c p)  /. {a -> 0.2,  c -> 0.004}

Plot[example, {p, -100, 100},  
        AxesLabel -> {"p", "dp/dt"}];

Hence, for positive population p0, the population grows indefinitely.

(iii)

DSolve[p'[t] == p[t] (a + c p[t]), p[t], t]

If we let k = (E^C)[1], a constant, then this can be expressed
	p(t) = ((-a)/(c - k (E^(-a t))))

But p(t) = ∞ if c - k (E^(-a t)) = 0.

Solve[k (E^(-a t)) == c, t]

Hence, an infinite population is reached in a finite time period

Question 7

(i)

The two fixed points are

Solve[r (M - p) (p - m) == 0, p]

The turning point is

Solve[D[r (M - p) (p - m), p] == 0, p]

D[D[r (M - p) (p - m), p], p]  /. p -> ((m + M)/2)

which is negative, hence the turning point is a maximum.

D[r (M - p) (p - m), p]  /.  p -> m

which is positive since M > m, so the fixed point m is locally unstable.

D[r (M - p) (p - m), p]  /.  p -> M

which is negative since M > m, so the fixed point M is locally stable.

The parameter m is associated with such factors as:
	(a) reproduction system
	(b) density of species (i.e., the area over which it operates).
The parameter M is associated with such factors as:
	(a) food availability
	(b) predation
	(c) catch

(ii)

In the case of (b) species can begin to die out well before p = 0, i.e., at p = m. It is important, therefore, to establish the value of m, e.g., the blue whale.

Question 8

(i)

DSolve[{p'[t] == r p[t] (a - Log[p[t]]), p[0] == p0}, p[t], 
      t]

PowerExpand[
      (E^((E^(-r t)) (-a + a (E^(r t)) + Log[p0])))]

FullSimplify[
      (E^((E^(-r t)) (-a + a (E^(r t)) + Log[p0])))]

which can be written
	\(TraditionalForm
      \`E\^\(-\(a(E\^\(\(-r\)\ t\) - 1)\)\)\ 
          E\^\(\(E\^\(\(-r\)\ t\)\) Log[p0]\) = 
        E\^\(-\(a(E\^\(\(-r\)\ t\) - 1)\)\)\ p0\^\(E\^\(\(-r\)\ t\)\)\)

since

(E^Log)[(p0^((E^(-r t))))]

(ii) and (iii)

The fixed points are,

Solve[r p (a - Log[p]) == 0, p]

Limit[r p (a - Log[p]), p -> 0]

Hence, there are two fixed points one at p = 0 and the other at (E^a).

D[r p (a - Log[p]), p]

The slope of \(TraditionalForm\`\(p\& . \)\) at the point p = 0 cannot be determined since this involves Log[0]. On the other hand,

Simplify[D[r p (a - Log[p]), p]  /. p -> (E^a)]

which is equal to -r, since Log[(E^a)] = a Log[E] = a. Since r is positive, then the slope of \(TraditionalForm\`\(p\& . \)\) at the point (E^a) is negative, and so the point p = (E^a) is locally stable.

Also \(TraditionalForm\`\(p\& . \)\) is at a maximum at the value

Solve[D[r p (a - Log[p]), p] == 0, p]

(iv)

Limit[(E^(-a ((E^(-r t)) - 1))) 
        (p0^((E^(-r t)))), t -> ∞]

So Mathematica cannot solve for the limit directly. But we note that the limit of (E^(-r t)) as t -> ∞ is zero. Hence, the limit is (E^a) as t -> ∞.

Question 9

DSolve[{x'[t] == -3 y[t], y'[t] == -9 x[t]}, {x[t], 
        y[t]}, t]

Since
	\(TraditionalForm
      \`dy\/dx = 
        \(\(dy/dt\)\/\(dx/dt\) = 
          \(\(y\& . \)\/\(x\& . \) = 
            \(\(\(-9\) x\)\/\(3  y\) = \(3  x\)\/y\)\)\)\)
then we can integrate by parts for we have,
	y dy  =  3 x dx

∫y \[DifferentialD]y

∫3  x \[DifferentialD]x

Solve[((y^2)/2) == ((3 (x^2))/2) + (c/2), y]

where c/2 is the constant of integration.

Alternatively, we can solve directly

DSolve[y'[x] == 3  x/y[x], y[x], x]

Which can be expressed,
	y(x) = -(Sqrt[3  (x^2) + c]) and y(x) = (Sqrt[3  (x^2) + c]), where c = 3  C[1].

Question 10

The \(TraditionalForm\`\(T\& . \) = 0\) and \(TraditionalForm\`\(B\& . \) = 0\) phase lines are

Solve[a (1 - (T/k1)) T - b T B == 0, B]

Apart[((a (k1 - T))/(b k1))]

Solve[c (1 - (B/k2)) B - d T B == 0, B]

Apart[((k2 (c - d T))/c)]

The two phase lines are then,
	B = (a/b) - (((a )/(b k1))) T	for \(TraditionalForm\`T\& .  = 0\)
	
	B = k2 - (((d k2 )/c)) T	for \(TraditionalForm\`B\& .  = 0\)
which intersect at the point

Simplify[
      Solve[{B == (a/b) - ((a T)/(b k1)), B == k2 - ((d k2 T)/c)}, {B, 
          T}]]

which is as far as we can take the general result with Mathematica.

Question 11

(i)

Given
	\(TraditionalForm
      \`\(N(t)\)\/K = 
        \(N0/K\)\/\(\((N0/K)\) + \((1 - \((N0/K)\) E\^\(\(-r\)\ t\)\)\)\)
then N(2)/K is

(0.25/(0.25 + (1 - 0.25) (E^(-(0.71) 2))))

and hence N(2) is

(0.579662) (80.5 (10^6))

or 46.6628×(10^6) Kg

(ii)

In solving for t, we write N(t) as Nt.

Solve[
      (Nt/K) == ((N0/K)/((N0/K) + (1 - (N0/K)) (E^(-r t)))), t]

which can be expressed
	\(TraditionalForm
      \`t = \(-\((1\/r)\)\) 
          Log[\(\((N0/K)\) \((1 - Nt/K)\)\)\/\(Nt/\(K(1 - N0/K)\)\)]\)

(a) t at which N(t) = 0.5  K or N(t)/K = 0.5.

-((1/0.71)) 
      Log[((0.25 (1 - 0.5))/(0.5 (1 - 0.25)))]

(b) t at which N(t) = 0.75  K or N(t)/K = 0.75.

-((1/0.71)) 
      Log[((0.25 (1 - 0.75))/(0.75 (1 - 0.25)))]