23
$\begingroup$

I have conducted a chemistry experiment and obtained a list of concentrations of a reactant against time. I plot a graph of this using Mathematica's ListPlot functionality, as follows:

ListPlot[{{0, 0.562}, {10, 0.523}, {20, 0.480}, {30, 0.438}, {40, 0.398}, {50, 0.357},
 {60, 0.320}, {70, 0.285}, {80, 0.255}, {90, 0.230}, {100, 0.220}}, Joined -> True,
 InterpolationOrder -> 2, Mesh -> Full, AxesLabel -> {"Time (s)", "[A] (moldm^-3)"}]

This produces the graph below:

enter image description here

However, I wish to calculate the rate of reaction, for which I need the first derivative at $t=0$. I want to be able to get the interpolating function used by the ListPlot function when I specify InterpolationOrder -> 2, is this possible?

Thanks in advance!

$\endgroup$
3

4 Answers 4

22
$\begingroup$

Like Chris says:

data = {{0, 0.562}, {10, 0.523}, {20, 0.480}, {30, 0.438}, {40, 0.398}, {50, 0.357}, 
        {60, 0.320}, {70, 0.285}, {80, 0.255}, {90, 0.230}, {100, 0.220}}
f = Interpolation[data, InterpolationOrder -> 2]

then

f'[0]

returns

-0.0037
$\endgroup$
4
  • $\begingroup$ Thank you to both you and Chris for your help! :) $\endgroup$ Sep 23, 2012 at 10:21
  • $\begingroup$ Is there anyway I can reclaim a function in terms of $t$ from the InterpolatingFunction? Thanks! $\endgroup$ Sep 23, 2012 at 13:24
  • $\begingroup$ @Shaktal - Have a look at Fit and FindFit, which will give you a function for the given data, in any form you want. If you think you have a quadratic equation then f := Fit[data, {x^2, x, 1}, x] $\endgroup$
    – stevenvh
    Sep 23, 2012 at 13:31
  • 1
    $\begingroup$ @Shaktal: you'll probably want to use FindFit[]; it looks to me here that your reaction is following first-order kinetics, and linearization (by taking logarithms) isn't the best way to get good values for your parameters. $\endgroup$ Sep 23, 2012 at 13:42
43
$\begingroup$

Let me elaborate on @stevenvh's answer using Splines instead of Interpolation.

The danger of using f'[0] is that the built-in interpolation requires that the (Hermite) polynomials go through each data points. Now if you data is noiseless that's fine, but if your data is noisy, the derivative of the interpolation will be all the more noisy (as a rule its always a bad idea to differentiate noise!).

In that case, you can either use a parametric fit (as suggested by @stevenvh) or a non-parametric fit, say, using Splines (or any parametrisation for which you can weight the relative degree of smoothness of your fit). Let me illustrate this, starting with your dataset:

 pts = {{0, 0.562}, {10, 0.523}, {20, 0.480}, {30, 0.438}, {40, 
0.398}, {50, 0.357}, {60, 0.320}, {70, 0.285}, {80, 0.255}, {90, 
0.230}, {100, 0.220}};

Let us define a noisy set:

 pts2 = pts + RandomReal[{-0.1, 0.1}, Dimensions[pts]];

Now the original solution was to use (Hermite) interpolation;

 f0 = Interpolation[pts2];

but let's also define a Spline interpolation

 f = BSplineFunction[pts2];

Comparing the two curves

 Show[Graphics[{Red, Point[pts2], Green, Line[pts2]}, Axes -> True], 
 ParametricPlot[f[t], {t, 0, 1}], 
 Plot[f0[x], {x, 0, 100}, PlotStyle -> {Dashed, Red}], 
 AspectRatio -> 1/2] // Quiet;

Mathematica graphics

we can visually see that the grey curve (the Spline) provides a better representation of the derivative at the origin. Indeed

  f'[0] // #[[2]]/#[[1]] &
 (* -0.00365569 *)

is closer to the underlying value of -0.0037, than the value returned by the original method (red dashed curve):

 f0'[0]
 (* 0.0120881 *).

Of course, statistically they both yield something reasonable

 Mean@ Table[
 pts2 = pts + RandomReal[{-0.1, 0.1}, Dimensions[pts]];
 f = BSplineFunction[pts2 ];
 f0 = Interpolation[pts2];
 {f0'[0], f'[0] // #[[2]]/#[[1]] &}
 , {15000}] // Quiet
(* {-0.00364452,-0.00393355}  *)

but the RMS on the Spline solution is typically $\sim$3 times smaller;

 (StandardDeviation@ Table[
 pts2 = pts + RandomReal[{-0.1, 0.1}, Dimensions[pts]];
 f = BSplineFunction[pts2 ];
 f0 = Interpolation[pts2];
 {f0'[0], f'[0] // #[[2]]/#[[1]] &}
 , {1500}] // Quiet) // #[[1]]/#[[2]] &
 (* 2.70277 *)

More generally, let us build on the ref/BSplineCurve Least square fitting example in the documentation. Consider the set of points:

 pts = Table[{i, Exp[-i] Sin[10 Pi i] + RandomReal[.3]}, {i, 0, 1, .003}];

Use uniform parametrization:

 uparam[pts_] := N[Range[0, 1, 1/(Length[pts] - 1)]];

Define a function to generate clamped knots for a given number of control points and degrees:

 kfun[n_, d_] := 
 Join[ConstantArray[0, d], Range[0, 1, 1/(n - d)], ConstantArray[1, d]];

Define the basis matrix for least squares:

 mbasis[pts_, n_, d_] := 
With[{param = uparam[pts]}, 
Table[BSplineBasis[{d, kfun[n, d]}, j - 1, param[[i]]], {i, 
 Length[param]}, {j, n}]];

A cubic B-Spline curve with 25 control points seems reasonable here for fitting that particular curve:

 ctrlpts = LeastSquares[mbasis[pts, 25, 3], pts];

Let me show the data with the curve:

 ListPlot[pts, Epilog -> {Red, BSplineCurve[ctrlpts, SplineDegree -> 3]}]

Mathematica graphics

Now the derivative at the origin is

 BSplineFunction[ctrlpts, SplineDegree -> 3]'[0] // #[[2]]/#[[1]] &
 (*  39.7595 *),

to be compared to the exact value:

 D[Exp[-x] Sin[10 Pi x], x] /. x -> 0 // N
 (* 31.4159 *)

which, given the amount of noise, is a reasonable estimate. On average, the non-parametric estimate is fine:

  {Mean[#], StandardDeviation[#]/Sqrt[Length[#]]} & @Table[
  pts = Table[{i, Exp[-i] Sin[10 Pi i] + 
  RandomVariate[NormalDistribution[0, 0.1]]}, {i, 0, 1, .01}];
  ctrlpts = LeastSquares[mbasis[pts, 25, 3], pts];
   BSplineFunction[ctrlpts, SplineDegree -> 3]'[0] // #[[2]]/#[[1]] &, {100}]

(* {32.378,1.25732} *)

The main virtue of this second example is that we now have a control parameter, namely the number of control points, which we can adapt to the level of noise in order to get an accurate estimate of the derivative. More generally, fitting with a variable amount of smoothing can be formalized in the context of Maximum a Posteriori analysis. There are even methods (such as cross validation) to estimate automatically what the proper amount of smoothing should be.

EDIT

As an alternative to using the number of Splines basis to control the smoothness of the curve, let us consider an explicit penalty function. The idea here is that instead of solving for the best (spline) weights in the least square sense (a maximum likelihood solution solution while assuming Gaussian statistics), we find these, subject to a prior corresponding to a 'roughness' penalty (which allows us to tune how smooth the spine function should be, which involves adding a tunable cost to unsmooth spline).

Clear[ctrlpts];
ctrlpts[lambda_:0]:= 
With[{mat = mbasis[pts, 50, 3], 
reg = SparseArray[{{i_, i_} -> 
    2., {i_, j_} /; Abs[i - j] == 1 -> -1.}, {50, 50}, 0.]}, 
LinearSolve[Transpose[mat].mat + 10^(lambda) Transpose[reg].reg, 
Transpose[mat].(Last /@ pts)]];

Let us start with an under-smooth solution (lambda=-4.25)

Show[ListPlot[pts, AxesLabel -> font[] /@ {x, y}], 
ListLinePlot[{First /@ pts, mbasis[pts, 50, 3].ctrlpts[-4.25]} // Transpose,
PlotStyle -> Red]]

Mathematica graphics

Note that the red curve seems to somewhat over fit the noise for this value of lambda.

Now let us put more weights on our penalty using ctrlpts[0.25]

 Show[ListPlot[pts, AxesLabel -> font[] /@ {x, y}], 
 ListLinePlot[{First /@ pts, mbasis[pts, 50, 3].ctrlpts[0.25]} // Transpose,
 PlotStyle -> Red]]

Mathematica graphics

Better?

As far as the estimated derivative is concerned, we can define with mathematica the derivative of our spline basis

Clear[dbasis];
dbasis[pts_, n_, d_] := 
With[{param = uparam[pts]}, 
Table[Release[D[BSplineBasis[{d, kfun[n, d]}, j - 1, x], x]] /. 
 x -> param[[i]], {i, Length[param]}, {j, n}]];

and watch how changing the weight of the penalty changes the roughness of the derivative

{ListLinePlot[{First /@ pts, dbasis[pts, 50, 3].ctrlpts[-6]}//Transpose, PlotStyle -> Blue],
 ListLinePlot[{First /@ pts, dbasis[pts, 50, 3].ctrlpts[0.15]}//Transpose,PlotStyle -> Red],
 Plot[D[Exp[-x] Sin[10 Pi x], x] // Release, {x, 0, 1},PlotStyle -> {Dashed,Green}]} // Show

Mathematica graphics

Note that the derivative at the origin is still available via

 BSplineFunction[ctrlpts[0.25], SplineDegree -> 3]'[0]

and more generally e.g. the second derivative

  df[x_] = BSplineFunction[ctrlpts[0.25], SplineDegree -> 3]''[x]

In closing, Wikipedia page on this topic is rather good.

UPDATE 2021

Since version 12 Mathematica has built in regularisation which is implemented for Interpolation here

$\endgroup$
1
  • $\begingroup$ Oh man great answer! Learned a bunch, thanks $\endgroup$
    – Gabriel
    Oct 24, 2012 at 6:10
7
$\begingroup$

To see the relevance of the comments by @chris for your data, plot the derivatives for two interpolations of your data:

f2 = Interpolation[data, InterpolationOrder -> 2];
f3 = Interpolation[data, InterpolationOrder -> 3];
p = Plot[{f2'[x], f3'[x]}, {x, 0, 100}, PlotStyle -> {Red, Blue}, 
AxesLabel -> {"Time", "Rate of Reaction"}]

Mathematica graphics

$\endgroup$
1
$\begingroup$

I think that the first problem is to know the physical law between the measured parameter and the time.

For the kinetic analysis of the A vitamin via the Carr-Price reaction (absorbance vs. time), according to the Beer-Lambert law, I use LinearModelFit. The intercept is the absorbance at t0.

model = LinearModelFit[data, x, x];
equation = Normal[model];
slope = Select[equation[[2]], # Epsilon Reals &, 1];
intercept = equation[[1]];
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.