Re: Re: Re: Finding derivatives of a list?
- To: mathgroup at smc.vnet.net
- Subject: [mg40919] Re: [mg40888] Re: [mg40854] Re: [mg40816] Finding derivatives of a list?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 24 Apr 2003 05:25:51 -0400 (EDT)
- References: <200304230517.BAA05828@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Selwyn Hollis wrote:
>
> I feel compelled to point out that using interpolation for this purpose
> is, in general, a Bad Idea.
>
> Please have a look at the following, which uses an example attributed
> to Runge. The Do loop plots the difference between f''[x] and the
> second derivative of the InterpolatingFunction based on 21 nodes -1,
> -.9, ..., .9, 1, for InterpolationOrders 3 through 20. You'll notice
> that the difference is quite large in all cases.
>
> runge[x_] := 1/(1 + 20 x^2);
> Plot[runge[x], {x, -1, 1}, Frame->True, PlotLabel->"Runge's example"];
>
> pts = Table[{x, runge[x]}, {x, -1., 1, 0.1}];
>
> Do[ interp = Interpolation[pts, InterpolationOrder -> n];
> Plot[Evaluate[{runge''[x] - interp''[x]}], {x, -1, 1},
> PlotRange -> {-10, 10},
> PlotLabel -> "InterpolationOrder --> " <> ToString[n]],
> {n, 3, 20}]
>
> -----
> Selwyn Hollis
> http://www.math.armstrong.edu/faculty/hollis
>
> On Tuesday, April 22, 2003, at 06:44 AM, Daniel Lichtblau wrote:
>
> > AES/newspost wrote:
> >>
> >> Specific problem is how to generate a list of values of the second
> >> derivative of a relatively smooth function at a set of equally spaced
> >> points, when the function itself is known only as a list of numerical
> >> values at those same points?
> >>
> >> --
> >> "Power tends to corrupt. Absolute power corrupts absolutely."
> >> Lord Acton (1834-1902)
> >> "Dependence on advertising tends to corrupt. Total dependence on
> >> advertising corrupts totally." (today's equivalent)
> >
> > Here are some possibilities.
> >
> > (i) Form an interpolation of relatively high order (say 6 or so). Take
> > second derivatives.
> >
> > (ii) Use finite differences to approximate the second derivatives.
> >
> > (iii) Use Fourier to get the approximated derivatives. See for example
> >
> > Wang, Jing. B (2002). Numerical differentiation using Fourier. The
> > Mathematica Journal 8:3. 383-388.
> >
> > I believe there was a small error in the code provided; you might want
> > to contact the author at wang at physics.uwa.edu.au
> >
> >
> > Daniel Lichtblau
> > Wolfram Research
It depends on how it is used, on the underlying assumptions of
smoothness of the data set, and on relative merits of alternative
approaches. Your example will serve to illustrate.
First let me say why I want to use "relatively high order" for the
interpolation (I'm sure you know this, but others may not). For low
order one simply does not get smooth second derivatives in the
interpolation. An order of five or six should suffice to get this much.
But now we run into another issue, as seen in your example. One must
have alot more points in order for the interpolation to be a useful
approximation at most points. For example, with interpolation order of
15 the function and its approximation do not agree at all well, and, not
surprisingly, the second derivatives disagree all the more. I think a
general rule of thumb is to use an order no higher than Sqrt[n] for
interpolating n points. Actually this is for polynomial interpolations
(to keep "wiggles" down), but I think the rule is also often reasonable
for piecewise interpolating functions such as those provided by
Interpolation.
Hence to get an order of 6, I'd generally want twice as many points as
you have. But that's not the real issue for your example (or so I
believe). What matters more is the amount of variation of that second
derivative relative to the point spacing. It's quite large. Hence in
some sense that point spacing violates the "smoothness" assumption. To
see what I mean, look at
Plot[runge'[x], {x, -.9, .9}]
near the origin. A spacing of 1/10 simply cannot capture this variation.
I'll illustrate by finding those derivatives using an alternative
approach of finite differencing. Before I take this further I should
point out that the comparison is mildly unfair insofar as I used a
reasonably high order interpolation, but only the most basic discrete
difference approximation. So you may want to try more careful discrete
approximations to see if they yield significantly better results than
those below.
dx = .1;
derivapprox = ListConvolve[{1,-2,1},pts[[All,2]]]/dx^2
innerpts = Take[pts[[All,1]],{2,-2}];
innerptderivs = Map[runge''[#]&, innerpts]
We now have the second derivatives and their approximations evaluated at
the interior points -.9,-.8,...,.9. If we plot them together we see they
do not agree well near the middle.
<<Graphics`MultipleListPlot`;
MultipleListPlot[innerptderivs,derivapprox]
Perhaps a better way is to make a plot akin to yours, where we first do
an interpolation (of default order) of the discretely approximated
second derivatives.
derivinterp = Interpolation[Transpose[{innerpts,derivapprox}]]
Plot[Evaluate[runge''[x] - derivinterp[x]], {x, -.9, .9}]
You will see that this also give significant variation, just as did the
second derivative of the sixth order interpolation.
Another point to make is that one really should check relative rather
than absolute error of these second derivatives. For that you might
compare
Plot[(runge''[x]-interp[6]''[x])/Sqrt[runge''[x]^2+interp[6]''[x]^2],
{x, -1, 1}, PlotRange -> {-10, 10}]
Plot[(runge''[x]-derivinterp[x])/Sqrt[runge''[x]^2+derivinterp[x]^2],
{x, -.9, .9}, PlotRange -> {-10, 10}]
These seem to give fairly similar pictures.
If you cut the spacing in half then both methods will give substantially
better results. Below is the code I used, in its entirety (I hope). My
guess, which I have not tried to verify, is that we really only need the
closer spacing near the center where the second derivative of the
underlying function is varying most rapidly.
dx = .05;
pts = Table[{x,runge[x]}, {x,-1.,1,dx}];
interp[6] = Interpolation[pts, InterpolationOrder->6];
Plot[Evaluate[runge[x] - interp[6][x]], {x, -1, 1}]
Plot[runge''[x]-interp[6]''[x], {x,-1,1}]
derivapprox = ListConvolve[{1,-2,1}, pts[[All,2]]] / dx^2;
innerpts = Take[pts[[All,1]], {2,-2}];
innerptderivs = Map[runge''[#]&, innerpts];
MultipleListPlot[innerptderivs, derivapprox]
derivinterp = Interpolation[Transpose[{innerpts,derivapprox}]];
Plot[Evaluate[runge''[x] - derivinterp[x]], {x,-1+dx,1-dx}];
Plot[(runge''[x]-interp[6]''[x])/Sqrt[runge''[x]^2+interp[6]''[x]^2],
{x,-1,1}, PlotRange -> {-10,10}];
Plot[(runge''[x]-derivinterp[x])/Sqrt[runge''[x]^2+derivinterp[x]^2],
{x,-1+dx,1-dx}, PlotRange -> {-10,10}];
I think in this case the prior interpolation followed by differentiation
gave the better results. A last way to assess the error is by
integration.
NIntegrate[Abs[runge''[x]-interp[6]''[x]]/Sqrt[runge''[x]^2+interp[6]''[x]^2],
{x,-1,1}]
For the differentiated interpolation function I get around .011.
NIntegrate[Abs[runge''[x]-derivinterp[x]]/Sqrt[runge''[x]^2+derivinterp[x]^2],
{x,-1+dx,1-dx}]
For the finite difference approximation function I get around .054. So
in this example the differentiated interpolation seems to perform better
(again, subject to the caveat that I did not try the anything beyond the
simplest finite differencing).
In conclusion, as to what method to use for approximating second
derivatives from a table of data, I'd have to say that it depends
heavily on underlying assumptions about smoothness, variation relative
to spacing, etc. And this does not even touch upon the assumptions (or
coding) necessary in order to make good use of a Fourier-based approach.
Daniel Lichtblau
Wolfram Research
- References:
- Re: Re: Finding derivatives of a list?
- From: Selwyn Hollis <selwynh@earthlink.net>
- Re: Re: Finding derivatives of a list?