[Date Index] [Thread Index] [Author Index]

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] ====