Second derivative

[gdargaud]

Hello all,

this is rather elementary numberical question, but I haven't done this in a while... Say I have 2 arrays: V[] and I[]. The 1st order derivative dI/dV can be computed for each point with:

Der1[i]=(I[i+1]-I[i])/(V[i+1]-V[i])

But what's the formula for second order derivatives ? i.e.: d2I/dV2 ?

I see the advanced library function difference(), but it only deals with dt. Is there a function that I can use ?

[RobertoBozzolo · Accepted solution]

Maybe too simplistic, but why not calulating the derivative of the derivative?

Der2[i]=(Der1[i+1]-Der1[i])/(V[i+1]-V[i])

[Wolfgang]

I am not using the simple form... but instead I adapted the function from the GNU scientific library

[gdargaud]

I should have added that it was my 1st idea, but I wasn't sure the divisor was correct. Also, shouldn't there be a /2 for each computation ?

[gdargaud]

Generally speaking, how would you want to filter the data before derivation ? There's little noise in it, but what little there is makes it hard on the 2nd derivative...So I was thinking of a lowpass filter.

And worse, since I derivate according to V, if V is not strictly growing (or decreasing), I get divide by zero. So maybe I also need a highpass filter!

Is there some generic data massaging protocol better than a large bandpass filter before doing a derivative ?

[gdargaud]

---

@Wolfgang wrote:

I am not using the simple form... but instead I adapted the function from the GNU scientific library

---

Ha, I'd missed the part about the 5-point vs 3-point and error estimate. But still, I don't see how to use gsl_deriv_central in the case of dI/dV where I and V are 2 arrays of values. Sorry, my signal processing skill have gone rusty in recent decades.

[Wolfgang]

well - I have slightly rewritten the routine because of some constraints; one is that I am also dealing with arrays - but the main idea remains the same, i.e., calculate more data points to allow for some error estimate which then in turn is used to find the optimum step size.

[gdargaud]

///////////////////////////////////////////////////////////////////////////////
/// HIFN	Compute dI/dV by finite difference method
/// HIFN	See https://en.wikipedia.org/wiki/Finite_difference_coefficient
/// HIPAR	I / Dividand input arrays
/// HIPAR	V / Divisor input array
/// HIPAR	n / Size of I and V
/// HIPAR	D / Output array, of size n-1
/// HIPAR	Method / 2 (forward)  or 3,5,7,9 (central)
/// OUT		D
///////////////////////////////////////////////////////////////////////////////
void Derivate(double I[], double V[], int n, double D[], int Points) {
// Accuracy 1 - Forward
	#define D2(A) (A[j]-A[j-1])
// Accuracy 2 - Central
	#define D3(A) (-A[j+1]+A[j-1])

// Accuracy 4 - Forward
//	#define EVAL5 (-25*I[j]+48*I[j-1]-36*I[j-2]+16*I[j-3]-3*I[j-4])/(-12*(V[j]-V[j-1] ??? ))
// Accuracy 4 - Central
	#define D5(A) (A[j+2]-8*A[j+1]+8*A[j-1]-A[j-2])
// Accuracy 6 - Central
	#define D7(A) (-A[j+3]+9*A[j+2]-45*A[j+1]+45*A[j-1]-9*A[j-2]+A[j-3])
// Accuracy 8 - Central
	#define D9(A) (3*A[j+4]-32*A[j+3]+168*A[j+2]-672*A[j+1]+672*A[j-1]-168*A[j-2]+32*A[j-3]-3*A[j-4])

	switch (Points) {
		case 2: for (int j=1; j<n  ; j++) D[j-1]=D2(I)/D2(V);                                          D[n-1]=0; break;
		case 3: for (int j=1; j<n-1; j++) D[j]  =D3(I)/D3(V); D[0]=                                    D[n-1]=0; break;
		case 5: for (int j=2; j<n-2; j++) D[j]  =D5(I)/D5(V); D[0]=D[1]=                        D[n-2]=D[n-1]=0; break;
		case 7: for (int j=3; j<n-3; j++) D[j]  =D7(I)/D7(V); D[0]=D[1]=D[2]=            D[n-3]=D[n-2]=D[n-1]=0; break;
		case 9: for (int j=4; j<n-4; j++) D[j]  =D9(I)/D9(V); D[0]=D[1]=D[2]=D[3]=D[n-4]=D[n-3]=D[n-2]=D[n-1]=0; break;
		default: break;
	}
}

///////////////////////////////////////////////////////////////////////////////
/// HIFN	Compute d2I/dV2 by finite difference method
/// HIPAR	Points / 3,5,7,9 (central), -2,-3,-5,-7,-9, iterate Derivate twice with |Points|
///////////////////////////////////////////////////////////////////////////////
void Derivate2(double I[], double V[], int n, double D[], int Points) {
// Accuracy 2 - Central
	#define DD3(A) (A[j+1]-2*A[j]+A[j-1])
// Accuracy 4 - Central
	#define DD5(A) (-A[j+2]+16*A[j+1]-30*A[j]+16*A[j-1]-A[j-2])
// Accuracy 6 - Central
	#define DD7(A) (2*A[j+3]-27*A[j+2]+270*A[j+1]-490*A[j]+270*A[j-1]-27*A[j-2]+2*A[j-3])
// Accuracy 8 - Central
	#define DD9(A) (-9*A[j+4]+128*A[j+3]-1008*A[j+2]+8064*A[j+1]-14350*A[j]+8064*A[j-1]-1008*A[j-2]+128*A[j-3]-9*A[j-4])

	switch (Points) {
		case -2: 
		case -3: 
		case -5: 
		case -7: 
		case -9: { double d[n]; Derivate(I, V, n, d, -Points); Derivate(d, V, n, D, -Points); } break; 

//		Those are wrong
//		case 3: for (int j=1; j<n-1; j++) D[j]=DD3(I)/DD3(V); D[0]=D[n-1]=0; break;
//		case 5: for (int j=2; j<n-2; j++) D[j]=DD5(I)/DD5(V); D[0]=D[1]=D[n-2]=D[n-1]=0; break;
//		case 7: for (int j=3; j<n-3; j++) D[j]=DD7(I)/DD7(V); D[0]=D[1]=D[2]=D[n-3]=D[n-2]=D[n-1]=0; break;
//		case 9: for (int j=4; j<n-4; j++) D[j]=DD9(I)/DD9(V); D[0]=D[1]=D[2]=D[3]=D[n-4]=D[n-3]=D[n-2]=D[n-1]=0; break;

		case 3: for (int j=1; j<n-1; j++) D[j]=4*DD3(I)/(     D3(V)*D3(V)); D[0]=                                    D[n-1]=0; break;
		case 5: for (int j=2; j<n-2; j++) D[j]=4*DD5(I)/(12*  D3(V)*D3(V)); D[0]=D[1]=                        D[n-2]=D[n-1]=0; break;
		case 7: for (int j=3; j<n-3; j++) D[j]=4*DD7(I)/(180* D3(V)*D3(V)); D[0]=D[1]=D[2]=            D[n-3]=D[n-2]=D[n-1]=0; break;
		case 9: for (int j=4; j<n-4; j++) D[j]=4*DD9(I)/(5040*D3(V)*D3(V)); D[0]=D[1]=D[2]=D[3]=D[n-4]=D[n-3]=D[n-2]=D[n-1]=0; break;
		default: break;
	}
}

#if 0
#include <stdio.h>

#define N 20
double f(double x) { 
//	return 10;	// 0 everywhere 
//	return x; 	// 1 for D, 0 for DD
//	return x*x; // 2x+1 for D2, 2x for Dn, 2 for DD
//	return x*x*x; 
	return exp(x); 
}

int main (int argc, char *argv[]) {
	double I[N], V[N], D2[N], D3[N], D5[N], D7[N], D9[N], DD2[N], DD3[N], DD5[N], DD7[N], DD9[N];
	for (int i=0; i<N; i++) I[i]=f(V[i]=i);
	Derivate (I, V, N, D2, 2);
	Derivate (I, V, N, D3, 3);
	Derivate (I, V, N, D5, 5);
	Derivate (I, V, N, D7, 7);
	Derivate (I, V, N, D9, 9);
	Derivate2(I, V, N, DD2, 2);
	Derivate2(I, V, N, DD3, 3);
	Derivate2(I, V, N, DD5, 5);
	Derivate2(I, V, N, DD7, 7);
	Derivate2(I, V, N, DD9, 9);
		printf("I\tV=x\tForw 3 5 7 9\tForw 3 5 7 9\n");
	for (int i=0; i<N; i++) 
		printf("%.4g %g\t%.4g %.4g %.4g %.4g %.4g\t%.4g %.4g %.4g %.4g %.4g\n",
			   I[i], V[i], D2[i], D3[i], D5[i], D7[i], D9[i], DD2[i], DD3[i], DD5[i], DD7[i], DD9[i]);
	return 0;
}
#endif

I've written the above based on the finite difference method. I'm pretty sure it's correct for the 1st derivative.

But for the 2nd derivative, what do I divide with ? The various formulas I see are always for dI/dt where dt is constant (actually 1), and in that case they divide by dt^2. But here since my dV is (slightly) variable, I don't know what to use.

For a 3-point central double derivative, it seems straightforward to use a 3-point central difference squared, but I'm not too sure about it or about the higher orders.