Re: Numerical Differentiation

*To*: mathgroup at smc.vnet.net*Subject*: [mg4281] Re: Numerical Differentiation*From*: tlm at ameslab.gov (Dr. T. L. Marchioro II)*Date*: Sat, 29 Jun 1996 03:53:47 -0400*Organization*: Iowa State University, Ames, Iowa*Sender*: owner-wri-mathgroup at wolfram.com

Mark James wrote: > Does anyone know of a function that calculates the derivative of > a function (that can't be differentiated symbolically) at a given > point by numerical means? I can't find it as a built-in or in the > standard packages. Thanks. This question arises a lot, either in this form, or in the more common one "how do I take the derivative of a series of data points." The answer is that there is no perfect way, unless you happen to know a fair amount about the data set, such as the applicability of polynomial fits, the frequency spectrum, etc. But overall, the best way to accomplish this is with the use of Distributed Approximating Functionals (usually called DAFs) which are a localized, numerically efficient way of representing L^2 integrable functions and linear operations (such as derivatives) on them. They are the single most useful "representation" I have ever encountered P and no, I did not invent them, although I've since discovering them I've done lots of work advancing them. DAFs provide an accurate representation of the function over a continuous interval. The error in the representation, whether on a grid point or off, is essentially constant throughout the interval, except at the endpoints (where edge effects occur. These can be obviated if need be.) Similarly, the derivative is represented equally well throughout the interval, with none of the anomalies that can occur with inappopriately applied polynomial fits or finite difference formulae. Below is a session demonstrating some of these features you might want to look through. I'm in the process of writing a couple of papers on DAFs for the Mma Journal which will present these and other aspects of DAFs, and which will include my standard DAF package. Unfortunately, due to the press of my "day job", I've been working on this paper for about a year now, sigh..... Anyway, the short answer is: generate a set of data points from your function (even spaced is easiest, but not necessary) and then apply the DAF derivative operator as given in my package (which I will be happy to send to you). Hope this proves helpful --- Tom -- Dr. Thomas L. Marchioro II Two-wheeled theoretical physicist Applied Mathematical Sciences 515-294-9779 Ames Laboratory 515-432-9142 (home) Ames, Iowa 50011 tlm at ameslab.gov Project Coordinator: Undergraduate Computational Engineering and Sciences http://uces.ameslab.gov/ >>>>>>>> dirac> math Mathematica 2.2 for NeXT Copyright 1988-93 Wolfram Research, Inc. -- NeXT graphics initialized -- (* Make up a "generic" function *) In[1]:= fun=Table[ E^(-x^2/3)*(.3x^4-.2x^3-x^2+x+.001)*Cos[x], {x, -7, 7, .25}]//N; (* Here is one possible form of DAF, based on Hermite Polynomials *) In[2]:= gdafH[ x_, m_, s_]:= 1/Sqrt[2*Pi*s^2]*E^(-x^2/(2*s^2))*Sum[HermiteH[2*j,x/Sqrt[2.*s^2]]* HermiteH[2*j,0]*(1/(2^(2*j)))*1/((2*j)!), {j, 0, m}] (* Load my standard package *) In[3]:= <<DAF.m (* Find optimal value of the width parameter *) In[4]:= s=optsig[18]*.25 Out[4]= 0.500025 (* Make a DAF for fitting the function, just to show how it works *) In[5]:= fitH= Table[gdafH[i, 18, s], {i, -5.5, 5.5, .25}]//N; In[6]:= fitH -16 -14 -13 -11 Out[6]= {7.59084 10 , -8.84531 10 , -1.22763 10 , 1.32137 10 , -11 -10 -9 -8 > -7.98243 10 , -4.71594 10 , 9.81288 10 , -6.69121 10 , -7 -7 -6 > 2.03673 10 , 4.10018 10 , -8.58167 10 , 0.0000566282, -0.00025614, > 0.000906552, -0.0026461, 0.00655072, -0.0139979, 0.0261361, -0.0430334, > 0.0629337, -0.0822126, 0.0963393, 3.89846, 0.0963393, -0.0822126, > 0.0629337, -0.0430334, 0.0261361, -0.0139979, 0.00655072, -0.0026461, -6 -7 > 0.000906552, -0.00025614, 0.0000566282, -8.58167 10 , 4.10018 10 , -7 -8 -9 -10 > 2.03673 10 , -6.69121 10 , 9.81288 10 , -4.71594 10 , -11 -11 -13 -14 > -7.98243 10 , 1.32137 10 , -1.22763 10 , -8.84531 10 , -16 > 7.59084 10 } In[7]:= SetParameters[s,.25] (* Apply the DAF to the Function *) In[8]:= ans=ToeplitzMatTimesVec[fitH, fun]; (* See how well it fits it. The relative error, except at the boundaries, is essentially constant, but non-zero, as it would be with interpolation formulae. Evaluation between grid points will be equally accurate. *) In[9]:= (ans-fun)/fun Out[9]= {-0.00536289, 0.00144985, -0.000398982, 0.000109871, -0.0000300339, -6 -6 -7 -7 -7 > 8.14213 10 , -2.20934 10 , 6.19937 10 , -2.02905 10 , 2.41878 10 , -9 -10 -11 -11 > 5.47186 10 , -3.04859 10 , 2.27914 10 , 4.82129 10 , -11 -11 -11 -11 > 3.44883 10 , -1.04593 10 , -6.19808 10 , -7.35209 10 , -11 -10 -9 -10 -9 > 2.57312 10 , 3.91886 10 , 2.8484 10 , 5.09101 10 , -4.05693 10 , -10 -11 -10 -10 > -6.12523 10 , 9.9343 10 , 4.23474 10 , 3.97645 10 , -10 -7 -9 -10 > -3.18627 10 , 3.19527 10 , 1.14463 10 , -5.85762 10 , -9 -9 -8 -8 -9 > -2.32052 10 , -3.59793 10 , 1.907 10 , 6.18858 10 , -2.81078 10 , -10 -10 -11 -11 > 5.67949 10 , 3.59834 10 , 7.3603 10 , -6.17788 10 , -11 -11 -11 -11 > -7.41647 10 , -2.47583 10 , 3.01013 10 , 5.53946 10 , -11 -10 -9 -7 -7 > 3.58665 10 , -3.2438 10 , 5.99097 10 , 2.61984 10 , -2.17367 10 , -7 -6 -6 > 6.57658 10 , -2.32213 10 , 8.48354 10 , -0.0000310387, 0.00011268, > -0.00040625, 0.00146632, -0.00538943} (* Okay, now do the derivatives *) In[10]:= vector=Table[i, {i, -7, 7, .25}]; (* Find the analytic derivative, available in this case *) In[11]:= dfun[x_]=D[E^(-x^2/3)*(.3x^4-.2x^3-x^2+x)*Cos[x], x]; (* Here is the DAF derivative operator. Because it is not an interpolation formula, it gives derivatives equally well throughout the interval *) In[12]:= dfitH=BandedDAFDeriv[20,18,1,optsig[18]*.25, 0,.25]; In[13]:= dfitH -12 -10 -9 -9 Out[13]= {-6.98418 10 , -1.07717 10 , 1.66257 10 , -6.28633 10 , -8 -7 -6 > -5.08179 10 , 8.3523 10 , -5.84927 10 , 0.0000255138, -0.0000647072, > -0.0000109105, 0.00108739, -0.00676197, 0.0281798, -0.0935711, 0.263526, > -0.651742, 1.45487, -3.02748, 6.2247, -15.0287, 0, 15.0287, -6.2247, > 3.02748, -1.45487, 0.651742, -0.263526, 0.0935711, -0.0281798, > 0.00676197, -0.00108739, 0.0000109105, 0.0000647072, -0.0000255138, -6 -7 -8 -9 -9 > 5.84927 10 , -8.3523 10 , 5.08179 10 , 6.28633 10 , -1.66257 10 , -10 -12 > 1.07717 10 , 6.98418 10 } (* Apply the Derivative Operator and check it. Think this might be good enough for you? *) In[14]:= check=ToeplitzMatTimesVec[dfitH, fun]; In[15]:= (check-dfun[vector])//N -6 -6 Out[15]= {0.0000407922, -0.0000166356, 8.10154 10 , -3.90769 10 , -6 -7 -7 -8 -8 > 1.79204 10 , -6.77261 10 , 3.38611 10 , 1.63462 10 , 1.60503 10 , -7 -6 -6 > -4.80046 10 , -1.88495 10 , -5.24074 10 , -0.000012069, > -0.0000241565, -0.0000427325, -0.0000669032, -0.000091552, -0.000105546, > -0.0000917361, -0.0000303707, 0.0000934259, 0.000279596, 0.000504596, > 0.000719816, 0.000861042, 0.000868396, 0.000710226, 0.000400452, -9 > -1.38499 10 , -0.000400454, -0.000710225, -0.000868393, -0.00086104, > -0.000719816, -0.000504597, -0.000279597, -0.0000934256, 0.0000303715, > 0.0000917366, 0.000105546, 0.0000915517, 0.000066903, 0.0000427325, -6 -6 -7 > 0.0000241566, 0.000012069, 5.24078 10 , 1.88495 10 , 4.79511 10 , -8 -8 -7 -7 -6 > -1.29019 10 , -2.9278 10 , -2.96181 10 , 5.59041 10 , -1.5025 10 , -6 -6 > 3.26685 10 , -6.77636 10 , 0.0000139131, -0.00003412} ==== [MESSAGE SEPARATOR] ====