Nth-order derivative of an arbitrary function

ta0kira · 12-31-2008, 06:55 PM

I've been trying to use Maxima to calculate the Nth-order derivative of a numerically-calculated function, which it obviously has a huge problem doing. This isn't any surprise since it would probably need to take N^2 samples of the function to differentiate N times, but there really isn't another way to do it.

Code:

 n
d    /   1     x        m-1           \
---  |  ---- · ∫ (x - y)    · f(y) dy |
  n  \  Γ(m)   0                      /
dx

You might recognize this as the fractional derivative of f(x) of order n-m. I've actually calculated this symbolically for the f(x) in question (with some very liberal interpretations,) but I'm trying to calculate it the formal way to see how close the symbolic version is. The symbolic version matches perfectly for integer orders, but I used a set of fit curves to find the in-between orders, which probably isn't an actual fractional derivative. Also, it only worked symbolically because I strategically avoided constants, but I'd like to be able to calculate this with an arbitrary f(x). I'm exporting it to gnuplot in the form of a pm3d splot with x on the x-axis, n-m on the y-axis, and the resulting value on the z-axis.

Has anyone been able to calculate something like this? I have Matlab sitting around ready to install, but I haven't ever used it. My other thought is to write a loop and/or use a matrix. Thanks.
ta0kira

Sergei Steshenko · 12-31-2008, 10:50 PM

You can't, generally speaking, calculate the n-the derivative because, generally speaking, not every function is that smooth.

ta0kira · 01-01-2009, 01:07 AM

Supposing it's smooth for x+[0,δ] where δ is some sample range, what would be the simplest way to calculate the Nth derivative? I'll mostly be dealing with harmonic functions and natural exponents, all in the real domain. My first thought is to sample N+1 times initially, then recursively difference pairs of two until the Nth order:

Code:

2^x | x=1, N=3, δ=0.1

d \ x	1.0	1.033	1.067	1.1	        →
0	2.0	2.047	2.095	2.144	Δ/(δ/3) →
1	1.402	1.435	1.469		Δ/(δ/3) →
2	0.983	1.006			Δ/(δ/3) →
3	0.690    (actual: ln(2)^3*2^x=0.666)

The only problem I see with this is that the entire process would need to be repeated for every point on the plot, not to mention the fact that using the function in my initial post would also require numerical calculation of the integral, giving me at least x*y*(N+1)/2 calculations of that integral. With my plot as it is, 80x50 & N=[0,8] that would be 18,000 numerical integral calculations! Of course, programatically I could carry over points and might be able to get away with ~9,000 calculations of the integral, but I'm not in the mood to write my own fractional-derivative mapper. I suppose that might be in order if I decide to continue this path of research. If not, it's still the sort of thing I would occupy myself with just to prove I haven't been beaten! It really only needs to be within about 10% at this point since I'm doing mostly visual analyses right now. Thanks.
ta0kira

PS Something I just thought of is that the functions I'll deal with will almost always be symbolically differentiable to any integer order (without converging to a constant,) so for real-order n I could calculate symbolically to order ⌈n⌉, then calculate the fractional integral ⌈n⌉-n of that, which should get me down to 4,000 numerical integrations. I don't know who's idea it was to differentiate after the fractional integral! Now I need to figure out how to get Maxima to numerically integrate for a plot (could actually be a gnuplot problem; I'll try exporting to a data file.)

PPS I understand that the issue stated initially isn't exactly what it actually is. Sorry.