Corn Seedling example ( Curve Fitting via the Criterion of Least Squares, by John Alexander):

In which we use non-linear regression on corn growth data to estimate the parameters of a differential equation. 

That differential equation then provides a solution in line with the corn growth data, and allows us estimate its growth at any moment in time.

Andy Long
Spring, 2018

The Data: 

You might remember this as a dataset that I thought Alexander should not have fit with an exponential (because of its S-shape). 

Take a look:

weeks=Table[2i,{i,1,10}];
weight={21,28, 58, 75, 170, 422, 706, 853, 924, 966};
corndata=Transpose[{weeks, weight}];
p1=ListPlot[corndata,PlotMarkers->{{▲,10}},PlotStyle->{Blue}]

We’re going to try to fit a non-linear model to the corn weight data (cw),

cw(t)=((a (e^(r t)) K)/(1-a+a (e^(r t))))

cw[t_]:=((a (E^(r t)) K)/(1-a+a (E^(r t))))

The reason why we chose this form of the logistic equation is that this is solution to the logistic differential equation

y’(t)= r (1 - ((y(t))/K)) y(t)
y(0)=a K

Simplify[D[cw[t],t]]
Simplify[r (1-((cw[t])/K)) cw[t]]

Here r is the growth parameter, K is the carrying capacity of the environment, and a is the fraction of the carrying capacity population at the outset.

One can rewrite cw(t) (using α=Log[a]) as

(( (e^(α+r t)) K)/(1-(e^α)+ (e^(α+r t))))

which can be written as

((  K)/((1-(e^α))(e^(-(α+r t)))+ 1))

which can be written as

((  K)/(((e^(-α))-1)(e^(-r t))+ 1))

which can be written as

((  K)/(((ce)^(-r t))+ 1))

From this form one can see that K is the "carrying capacity" (horizontal asymptote), that r is a rate constant, and that when t=0, the initial population is at 
(Ke^α) =Ka
(that is,the initial population). 

Usually we're interested in when ((ce)^(-r t)) is equal to 1, since that is when the thing modeled hits half the carrying capacity (in this context it means the corn is half its ultimate weight). So r's job is to determine this.

From the data we can guess that K is around 1000, so that leads us to the approximate model

cw(t)=((  1000)/(((ce)^(-r t))+ 1))

Let’s linearize the model to get a good guess to the parameters as starting values for the non-linear regression:

((ce)^(-r (Subscript[t, i])))=((  1000)/(Subscript[cw, i]))-1

or

ln(c)-(Subscript[rt, i])=ln(((  1000)/(Subscript[cw, i]))-1)

transdata=Transpose[{weeks,Log[Abs[1000.0/weight-1]]}]
p2=ListPlot[transdata];
tm=LinearModelFit[transdata,x,x]
Show[p2,Plot[tm[x],{x,0,20}]]
{a,b}=tm["BestFitParameters"]
Show[p1,Plot[1000/(1+ Exp[a+b x]),{x,0,20}]]
{astart,bstart}={a,b}

Perform non-linear regression:

Clear[a,b,c]
nlm=NonlinearModelFit[corndata,c /(1+(1-Exp[-a])Exp[(a+b x)]),{{a ,astart},{b,bstart},{c,1000}},x]
nlm["ParameterTable"]

Show[p1,Plot[nlm[x],{x,0,20}]]

Match parameters to ODE for solution:

Clear[r,y,t,a,K]
LogisticEquation = y'[t]== r(1-(y[t]/K))*y[t];
sol=DSolve[{LogisticEquation, y[0]==a*K}, y,t]

numSoln=Evaluate[Table[y[t]/.sol[[1]]/.{K-> 970, a-> Exp[-7.21249], r-> 0.579798}, {t,0,20}]]
p2=ListPlot[Transpose[{Table[t, {t,0,20}],numSoln}],PlotMarkers->{{■,10}},PlotStyle->{Red}];
Show[p2,p1]