Shuttlecocks and the Akaike Information Criterion

Kurt Bryan and SIMIODE

A notebook to help explore the project in Section 3.5.4.

The Data: First, the data for the shuttlecock's fall, in time (seconds)/distance (meters) pairs:

shuttledata = {{0, 0}, {0.347, 0.61}, {0.47, 1.00}, {0.519, 1.22}, {0.582, 1.52}, {0.650, 1.83}, {0.674, 2.00}, {0.717, 2.13}, {0.766, 2.44}, {0.823, 2.74}, {0.870, 3.00}, {1.031, 4.00}, {1.193, 5.00}, {1.354, 6.00}, {1.501, 7.00}, {1.726, 8.50}, {1.873, 9.50}}

Number of data points

n = Length[shuttledata]

A plot

plt1 = ListPlot[shuttledata,AxesLabel->{"Time (seconds)","Distance (meters)"}]

The Model: We might posit a model of the form v'(t) = g (no air resistance) and consider g as an unknown, to be estimated. Then the governing ODE is (from equation (3.68) in the text)

de = v'[t] == g

The solution with v(0) = 0 is

sol = DSolve[{de,v[0]==0},v,t]

Define this as a function of t

vsol = v/.sol[[1]]

Integrate to find position, using x(0) = 0 where x(t) is the distance the shuttlecock has fallen:

dis= Integrate[vsol[tau],{tau,0,t}]

Make this into a function of t

x[t_] = dis

Estimating Parameters: Form a sum of squares

SS = Sum[(x[shuttledata[[i,1]]]-shuttledata[[i,2]])^2,{i,1,n}]

Minimize in g. First, a plot

Plot[SS,{g,0,15}]

Solve SS'(g) = 0 to find the least-squares estimate for gravitational acceleration

eq = D[SS,g] == 0;
bestg = Solve[eq,g]

The residual is

SS/.bestg

A plot to compare the fit of this model to the data:

plt2 = Plot[x[t]/.bestg, {t,0,1.873},PlotStyle->{Red}];
Show[plt1,plt2]