11-18-2016 09:10 AM
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 ?
Solved! Go to Solution.
11-18-2016 09:25 AM
Maybe too simplistic, but why not calulating the derivative of the derivative?
Der2[i]=(Der1[i+1]-Der1[i])/(V[i+1]-V[i])
11-18-2016 09:54 AM
I am not using the simple form... but instead I adapted the function from the GNU scientific library
11-18-2016 09:54 AM
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 ?
11-22-2016 03:41 AM
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 ?
11-22-2016 06:50 AM
@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.
11-22-2016 07:25 AM
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.
12-01-2016 05:35 AM - edited 12-01-2016 05:38 AM
/////////////////////////////////////////////////////////////////////////////// /// 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.