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>