Weekly Assignment #5

Instructions:  integrations in problems 1 and 2 should be evaluated by hand. Some numerical integrations should also be done by hand (whenever n ≤  4); if n>4, you may use technology. Show your work!

1. Average value

Let f(x)=x((e^(2x))+(e^(3x))).

Define and plot f(x) on the interval [0,1/2].

f[x_]:=x(Exp[2 x]+Exp[3 x])
a=0;
b=1/2;
Plot[f[x],{x,a,b}]

Demonstrate integration by parts to find the true average value of f(x) on [0,1/2]. Give the name “true” to this average value.

true=1/(b-a)Integrate[f[x],{x,a,b}]
N[true]

Use the Trapezoidal and Midpoint Rules with n=2 to approximate the average value over the interval [0,1/2] (calls these values “trap” and “mid”). Combine them in an appropriate way to give the Simpson’s Rule (“simp”) approximation, with n=4 (note: this n=4 value (Subscript[S, 4]) uses a weighted average of (Subscript[T, 2])and (Subscript[M, 2])).

n=2;
deltaX=(b-a)/n;
trap=1.0/(b-a)*deltaX/2*(f[a]+2f[a+deltaX]+f[b])
mid=1.0/(b-a)*deltaX*(f[a+deltaX/2]+f[a+3/2deltaX])
simp=(2*mid + trap)/3

Compute the errors (true - estimate) for the Trapezoidal and Midpoint rules. By plotting and bounding the absolute value of the second derivative of f(x), demonstrate that the approximations are within the error bounds for each method. Demonstrate that the errors are of opposite sign, as is typical.

Plot[Abs[f''[x]],{x,a,b}]
K=N[f''[b]]/(b-a) (* looks like we can take the right endpoint value for K; notice that I divide by (b-a), since that's actually the integral that is being approximated! The integrand can be considered f(x)/(b-a).   *)
trapBound=K(b-a)^3/(12 n^2);
trapError=true-trap;
{trapError,trapBound}
midBound=K(b-a)^3/(24 n^2);
midError=true-mid;
{midError,midBound}

Compute the absolute error (|true - estimate|) for the Simpson’s Rule approximation. By plotting the absolute value of the fourth derivative of f(x) (choose K as big as, or bigger than, the largest value), demonstrate that the approximations are within the error bounds for Simpson’s Rule.

Plot[Abs[f''''[x]],{x,a,b}]
K=N[f''''[b]]/(b-a) (* Looks like it's greatest on the right endpoint *)
bound=K(b-a)^5/(180 (2n)^4);
simpError=true-simp;
{simpError,bound}

2. Area and Volume

Let R be the region in the first quadrant of the plane between the curve y=f(x)=(Sqrt[4-(x^2)]) and the x-axis.

Define and plot f(x) on the interval [0,2].

f[x_]:=Sqrt[4-x^2]
a=0;
b=2;
Plot[f[x],{x,0,2}]

Find the area of R, by computing the appropriate integral using trigonometric substitution.

true=Integrate[f[x],{x,a,b}]

Estimate the area using (Subscript[S, 10]).

trapezoidal[exp_,{x_,a_,b_,n_}]:=Module[{dx=N[(b-a)/n],f},
f[u_]:=exp/.x->u;
Return[0.5dx(f[a]+f[b]+2Sum[f[a+k dx],{k,1,n-1}])]
]
midpoint[exp_,{x_,a_,b_,n_}]:=Module[{dx=N[(b-a)/n],f},
f[u_]:=exp/.x->u;
Return[dx Sum[f[a+k dx - 1/2dx],{k,1,n}]]
]
simpson[exp_,{x_,a_,b_,n_}]:=If[EvenQ[n],(trapezoidal[exp,{x,a,b,n/2}]+2midpoint[exp,{x,a,b,n/2}])/3,
Print["Sorry: Simpson's Rule only works on an even number of sub-intervals."]
]
simp=simpson[f[x],{x,a,b,10}]

Compute the absolute error in the approximation.

simpError=Abs[true-simp]

Compute the fourth derivative of f(x). Why can we not use the error estimate for Simpson’s rule in this case? (Is it possible to bound the fourth derivative on this interval?)

Simplify[f''''[x]]
Plot[Abs[f''''[x]],{x,a,b}]

3. Application

When we die, our bodies become rigid (rigor mortis sets in). Niderkorn’s (1872) observations on 113 bodies provides the main reference database for the development of rigor mortis.

The data: 

(({{hour, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13}, {proportion in rigor mortis, (2/123), (16/123), (47/123), (61/123), (81/123), (92/123), (99/123), (103/123), (110/123), (111/123), (111/123), (113/123)}}))

bodies={2,14,31,14,20,11,7,4,7,1,0,2};
(* bodies transitioning to rigor mortis at the associated hour: *)hours=Table[i,{i,2,13}];
cumulative=Accumulate[bodies];
data=Table[{hours[[i]],cumulative[[i]]/113},{i,1,12}];
MatrixForm[Transpose[data],TableHeadings->{{"hour","proportion in rigor mortis"},None}]
dataplot= ListPlot[data,PlotRange->{{0,15},{0,1.1}},PlotMarkers->{"x"},
AxesLabel->{"Hour","Number of Bodies"},
PlotLabel->"Number of bodies in rigor mortis by a given hour"]

One can fit a lovely model to this somewhat unlovely data, for the proportion p(t) of bodies in complete rigor mortis after t hours. It is illustrated in the graph below:

The model is (e^(-(22.47/(t^2.216)))): that is,

p[t_]:=E^(-22.47t^(-2.216))

a. Compute the average proportion of bodies in rigor mortis in the time interval from 6 to 10 hours after death, based on this model (write the integral, but you may use your calculator or Mathematica to produce your answer!).

NIntegrate[p[t]/(10-6),{t,6,10}]
(* So the average proportion of bodies you'd expect to see in the 6-10 hour mark passed the time of death is about 79%. *)

b. Now compute an approximation using only the data, and Simpson’s rule with n=4 (by hand). How do they compare?

a=6;
b=10;
n=4;
deltaX=(b-a)/n;
simp=1/(b-a)( deltaX /3.0(81+4*92+2*99+4*103+110)/123)