POSSIBLE SOLUTION(S)
1. If p = 0 or p = K then dp/dt = 0. p = 0 is uninteresting as we do not believe in spontaneous generation and thus if the population got to be 0 there is nothing more to study. However, we see that if K = 0 there is not any growth or decline, and certainly that would qualify K for something special, although it would not necessarily mean it is a limiting value or carrying capacity. We turn to question 2 for more insight.
2. If the population were less than K then dp/dt would be positive and the population would then increase. If the population were greater than K then dp/dt would be negative and the population would tend to fall. Thus the population would grow until it got to K and fall back to K if it were over K. This says that the population would not stay above K long and would tend to rise to K if it falls below K, thus limiting value or carrying capacity seems to be a good name.
3. Looking at the function dp/dt = r p (K - p)/K we see that it is a quadratic in p and has a maximum when p = K/2. thus dp/dt is greatest when the population is at a level of K/2 or one-half its carrying capacity.
4. We enter the data for the number of protozoa in a culture at time t in days.
Input :=
data = {{0, 40}, {1, 53}, {2, 70}, {3, 91}, {4, 118}, {5, 152},
{6, 193}, {7, 240}, {8, 293}, {9, 351}, {10, 411}, {11, 470},
{12, 526}, {13, 577}, {14, 622}, {15, 660}, {16, 691},
{17, 716}, {18, 736}, {19, 752}, {20, 764}};
This recovers r and K using a parabolic fit to the differential equation and non-symmetric differences.
Here is the growth rate data as a function of population size.
Input :=
ddata = Table[{data[[i]][[2]],
(data[[i+1]][[2]]- data[[i]][[2]])/
(data[[i+1]][[1]]- data[[i]][[1]])},
{i,1, Length[data]-1}];
Input :=
dp1 = ListPlot[ddata,PlotStyle->{PointSize[.02]},
AxesLabel->{"P", "P'(t)"},PlotRange->{{0,800},{0,80}}]
Output =
-Graphics-
This appears to be a parabolic shape and we attempt to fit this data to a parabola of the form a p - b p^2.
Input :=
q[p_] = Fit[ddata,{p,p^2},p]
Output =
2
0.310469 p - 0.000396418 p
Input :=
qp = Plot[q[p],{p,0,800}]
Output =
-Graphics-
We plot the data over the parabola and it appears a good match.
Input :=
Show[dp1,qp]
Output =
-Graphics-
We now use information about the function a p - b p^2 for p'(t) at t = 0 to obtain r and k values.
Input :=
eq1 = q'[0] == r; eq2 = q''[0] == - 2 r/K;
Input :=
Solve[{eq1,eq2},{r,K}]
Output =
{{K -> 783.185, r -> 0.310469}}
These values are K = 783.185 and r = 0.3104.
This recovers r and K using a parabolic fit to the differential equation and symmetric differences.
We now take averages of the slope about each time interval to get an averaging of slopes.
Input :=
ddata = Table[{data[[i]][[2]],
((data[[i+1]][[2]]- data[[i]][[2]])/
(data[[i+1]][[1]]- data[[i]][[1]]) +
(data[[i]][[2]]- data[[i-1]][[2]])/
(data[[i]][[1]]- data[[i-1]][[1]]))/2},
{i,2, Length[data]-1}]//N;
Input :=
dp1 = ListPlot[ddata,PlotStyle->{PointSize[.02]},
AxesLabel->{"P", "P'(t)"}]
Output =
-Graphics-
This appears to be a parabolic shape and we attempt to fit this data to a parabola of the form a p - b p^2.
Input :=
dq[p_] = Fit[ddata,{p,p^2},p]
Output =
2
0.298494 p - 0.00037285 p
Input :=
dqp = Plot[dq[p],{p,0,800}]
Output =
-Graphics-
We plot the data over the parabola and it appears a good match.
Input :=
Show[dp1,dqp]
Output =
-Graphics-
We now use information about the function a p - b p^2 for p'(t) at t = 0 to obtain r and k values.
Input :=
eq1 = dq'[0] == r; eq2 = dq''[0] == - 2 r/K;
Input :=
Solve[{eq1,eq2},{r,K}]
Output =
{{K -> 800.575, r -> 0.298494}}
These values K = 800.575 and r = 0.2984 compared to the answers we get using non-symmetric difference K = 783.185 and r = 0.3104.
We have now recovered r and K using a linear fit to the differential equation divided by p(t) and non-symmetric differences.
We now take averages of the slope about each time interval and divide these by p(t) for each time t.
Input :=
rdata = Table[{data[[i]][[2]],
(data[[i+1]][[2]]- data[[i]][[2]])/
(data[[i+1]][[1]]- data[[i]][[1]])/data[[i]][[2]]},
{i,1, Length[data]-1}]//N;
Input :=
rp1 = ListPlot[rdata,PlotStyle->{PointSize[.02]},
AxesLabel->{"P(t)", "P'(t)/P(t)"}]
Output =
-Graphics-
This appears to be a straight line and we attempt to fit this data to a straight line of the form a - b p.
Input :=
rq[p_] = Fit[rdata,{1,p},p]
Output =
0.332564 - 0.000434191 p
Input :=
rqp = Plot[rq[p],{p,0,800}]
Output =
-Graphics-
Input :=
Show[rp1,rqp]
Output =
-Graphics-
We now use information about the function a - b p for p'(t) at t = 0 to obtain r and k values.
Input :=
eq1 = rq[0] == r; eq2 = rq'[0] == - r/K;
Input :=
Solve[{eq1,eq2},{r,K}]
Output =
{{K -> 765.941, r -> 0.332564}}
This compares favorably to the original toy values of K = 800, r = 0.3 but not as well as several previous results.
This recovers r and K using a linear fit to the differential equation divided by p(t) and symmetric differences.
We now take averages of the slope about each time interval and divide these by p(t) at each time t.
Input :=
sdata = Table[{data[[i]][[2]],
((data[[i+1]][[2]]- data[[i]][[2]])/
(data[[i+1]][[1]]- data[[i]][[1]]) +
(data[[i]][[2]]- data[[i-1]][[2]])/
(data[[i]][[1]]- data[[i-1]][[1]]))/2/data[[i]][[2]]},
{i,2, Length[data]-1}]//N;
Input :=
sp1 = ListPlot[sdata,PlotStyle->{PointSize[.02]},
AxesLabel->{"P(t)", "P'(t)/P(t)"}]
Output =
-Graphics-
This appears to be a straight line and we attempt to fit this data to a straight line of the form a - b p.
Input :=
sq[p_] = Fit[sdata,{1,p},p]
Output =
0.300345 - 0.000375982 p
Input :=
sqp = Plot[sq[p],{p,0,800}]
Output =
-Graphics-
Input :=
Show[sp1,sqp]
Output =
-Graphics-
We now use information about the function a - b p for p'(t) at t = 0 to obtain r and k values.
Input :=
eq1 = sq[0] == r; eq2 = sq'[0] == - r/K;
Input :=
Solve[{eq1,eq2},{r,K}]
Output =
{{K -> 798.827, r -> 0.300345}}
This compares favorably to the original toy values of K = 800, r = 0.3 but not as well as several previous results. And it is better than the non-symmetric difference approach which yielded K = 765.941 and
r = 0.332564.
This recovers r and K after solving the differential equation and directly doing a nonlinear least squares fit to determine r and K.
This approach to estimating the parameters uses solution of the differential equation followed by a least squares, non-linear fit of the closed-form solution with parameters to the data.
Here is our differential equation model for paramecium growth.
Input :=
eq = p'[t]== r p[t] (k - p[t])/k
Output =
r (k - p[t]) p[t]
p'[t] == -----------------
k
We solve with the parameters r and k unknown.
Input :=
sol = DSolve[eq,p[t],t]
Output =
r t
E k
{{p[t] -> -------------}, {p[t] -> 0}}
r t
E + k C[1]
We capture the solution.
Input :=
P[t_] = p[t]/.sol[[1]]
Output =
r t
E k
-------------
r t
E + k C[1]
We use the initial condition to determine the constant of integration and a "final" form for the population model.
Input :=
sc = Solve[P[0]==40, C[1]]
Output =
-40 + k
{{C[1] -> -------}}
40 k
Input :=
Pop[t_] = P[t]/.sc[[1]]
Output =
r t
E k
--------------
r t -40 + k
E + -------
40
Input :=
Pop[t_] = Simplify[Pop[t]]
Output =
r t
E k
--------------
r t k
-1 + E + --
40
Here is the data from the observations.
Input :=
data
Output =
{{0, 40}, {1, 53}, {2, 70}, {3, 91}, {4, 118}, {5, 152}, {6, 193},
{7, 240}, {8, 293}, {9, 351}, {10, 411}, {11, 470}, {12, 526},
{13, 577}, {14, 622}, {15, 660}, {16, 691}, {17, 716}, {18, 736},
{19, 752}, {20, 764}}
We offer a plot of the data to see the form we are looking for.
Input :=
p1 = ListPlot[data,PlotStyle->{PointSize[.02]},
AxesLabel->{"t", "p(t)"}]
Output =
-Graphics-
This is the sum of the squares for the differences between the observed data(i.e. the second coordinate of the ith data point - data[[i]] [[2]]) and the predicted data from the model (i.e. the value of the function at the first coordinate of the ith data point - Pop[data[[i]] [[1]] ]) .
Input :=
SS[r_,k_] = Sum[ (data[[i]][[2]] - Pop[data[[i]][[1]]])^2,
{i,1,Length[data]}];
Here we enter the function with all three variable, r, k , and t.
Input :=
ParPop[r_,k_,t_] = Pop[t]
Output =
r t
E k
--------------
r t k
-1 + E + --
40
We offer up a way to plot the data with a trial for any value of r and k one wishes to try.
Input :=
try[r_,k_] := Show[p1,
Plot[ParPop[r,k,t],
{t,0,20},PlotRange->{0,800},
DisplayFunction->Identity] ]
Input :=
try[.35,750]
Output =
-Graphics-
We seek to minimize this sum of the squares function SS[r, k] of two variables (the parameters r and k).
Input :=
FindMinimum[SS[r,k],{r,.6},{k,600}]
Output =
{2.45308, {r -> 0.29972, k -> 799.791}}
Input :=
vals = FindMinimum[SS[r,k],{r,.2},{k,900}]
Output =
{2.45308, {r -> 0.29972, k -> 799.791}}
These compare very favorably to the toy values of K = 800 and r = 0.30.
We pick off the values of r and K.
Input :=
r1 = r/.vals[[2]];
Input :=
k1 = k/.vals[[2]];
We substitute the values of the parameters in the model and plot data and model to see how good we are doing.
Input :=
Protozoa[t_] = Pop[t]/.vals[[2]]
Output =
0.29972 t
799.791 E
--------------------
0.29972 t
18.9948 + E
Input :=
ProtozoaPlot= Plot[Protozoa[t],{t,0,20}]
Output =
-Graphics-
Input :=
Show[p1, ProtozoaPlot]
Output =
-Graphics-
This looks very reasonable for prediction purposes and these values of
r = 0.2999 and K = 799.79 may be used to characterize this species under these conditions.
Summary of values obtained for r and K and comparisons.
After using a parabolic fit to the differential equation and non-symmetric differences we determine r = 0.3104 and K = 783.184 and a sum of squares error of 1880.42 in the actual model.
Input :=
SS[.3104,783.184]
Output =
1880.42
After using a parabolic fit to the differential equation and symmetric differences we determine r = 0.2984 and K = 800.574 and a sum of squares error of 51.4211 in the actual model.
Input :=
SS[0.2984,800.574]
Output =
51.4211
After using a linear fit to the differential equation divided by p(t) and non-symmetric differences we determine r = 0.332564 and K = 765.941 and a sum of squares error of 18636.2 in the actual model.
Input :=
SS[.332564,765.941]
Output =
18636.2
After using a linear fit to the differential equation divided by p(t) and symmetric differences we determine r = 0.3003 and K = 798.826 and
and a sum of squares error of 18636.2 in the actual model.
Input :=
SS[.3003,798.826]
Output =
8.23376
After solving the differential equation and directly doing a nonlinear least squares fit we determine r = 0.2999 and K = 799.79 and a sum of squares error of 22.562 in the actual model .
Input :=
SS[.299,799.79]//N
Output =
22.5652
Summary: It would appear that the linear fit of the differential equation divided by p(t) using symmetric difference gives the best least squares fit although the direct non-linear least squares on the solutions does quite well too. Indeed both symmetric numerical methods do quite well in comparison with the direct least squares method on the solution. The non-symmetric differencing does worst.