Re: Finding derivatives of a list?
- To: mathgroup at smc.vnet.net
- Subject: [mg40968] Re: [mg40919] Finding derivatives of a list?
- From: Selwyn Hollis <selwynh at earthlink.net>
- Date: Fri, 25 Apr 2003 08:08:32 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
My previous "Bad Idea" comment was definitely an exaggeration ... In fact, there's really no difference between differentiating interpolating polynomials and using finite differences. (So why do the extra work of computing the interpolation?) However, with "black box" piecewise interpolation --- even with a smart choice of degree --- it's not clear to me that one can avoid discontinuities in its derivatives. All this begs the question: Just what are we getting from Interpolation[ pts, InterpolationOrder -> n]? It's clearly not a spline, yet it does often tend to be smooth. I've also noticed that it's common for the first derivative to have jump discontinuities while the second derivative is "smooth." Could we get some clarification on how Interpolation works? Thanks. ----- Selwyn Hollis http://www.math.armstrong.edu/faculty/hollis On Thursday, April 24, 2003, at 05:25 AM, Daniel Lichtblau wrote: > 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 > >