First Derivative in C++?

[CroCo]

deleted

[CroCo]

In Matlab, I have this function
t = 0:0.2:10
x = sin(t)
x' = diff(x)

x          x'
-------------------------
0       0.1987
0.1987	0.1907
0.3894	0.1752
0.5646	0.1527
0.7174	0.1241
0.8415	0.0906
0.932	0.0534
0.9854	0.0141
0.9996	-0.0257
0.9738	-0.0646
0.9093	-0.1008
0.8085	-0.133
0.6755	-0.16
0.5155	-0.1805
0.335	-0.1939
0.1411	-0.1995
-0.0584	-0.1972
-0.2555	-0.187
-0.4425	-0.1693
-0.6119	-0.1449
-0.7568	-0.1148
-0.8716	-0.08
-0.9516	-0.0421
-0.9937	-0.0025
-0.9962	0.0372
-0.9589	0.0755
-0.8835	0.1107
-0.7728	0.1415
-0.6313	0.1667
-0.4646	0.1852
-0.2794	0.1963
-0.0831	0.1996
0.1165	0.195
0.3115	0.1826
0.4941	0.1629
0.657	0.1367
0.7937	0.105
0.8987	0.0692
0.9679	0.0306
0.9985	-0.0092
0.9894	-0.0486
0.9407	-0.0861
0.8546	-0.1202
0.7344	-0.1495
0.5849	-0.1728
0.4121	-0.1892
0.2229	-0.1981
0.0248	-0.1991
-0.1743	-0.1922
-0.3665	-0.1775
-0.544	

Now I want to test that in C++. But with the function in my code for each 5 data there is one data for the derivative. How that possible?

[CroCo]

This is my code

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <fstream>


double First_Derivative(double a[5])
{
    double firstDer(0.0);
    double t(0.001);
    
    // no loop here as 5 points is just enough to calculate the derivative
    // of one elem (the middle one). If you do use a loop, it has to start
    // at 2 and end 2 elems before the end (so i - 2 and i + 2 are always
    // valid indexes for the array.)
    
    size_t i(2);
    
    firstDer = (   4.0 / 3.0 * (a[i + 1] - a[i - 1]) / (2.0 * t)
                - 1.0 / 3.0 * (a[i + 2] - a[i - 2]) / (4.0 * t) );
    // Rather than use (double)4 you can just use 4.0. The .0 tells the
    // compiler you mean a double (if you need a float instead, append
    // an f, like 3.0f )
    // std::cout << firstDer << std::endl;
    return firstDer;
}


int main(void)
{
    std::ifstream xdata;
    std::ofstream xdot;
    double temp[5] = {0};
    std::vector<double> x, xd;
   
    
    xdata.open("/Users/bandarhd/Desktop/TestX.txt");
    xdot.open("/Users/bandarhd/Desktop/TestXdot.txt");
    
    if (xdata.fail())
    {
        std::cout << "Error in open xdata file";
    }
    if (xdot.fail())
    {
        std::cout << "Error in Opening xdot file";
    }
    
    double line;
    while ( !xdata.eof() )
    {
        xdata >> line;
        x.push_back(line);
        
    }
        
//    for (int i = 0; i < x.size(); ++i)
//    {
//        std::cout << x[i] << std::endl;
//    }
    
    int xdata_counter = 0;
    while ( xdata_counter < x.size() )
    {
        
        for (int i = 0; i < 5; ++i)
        {
            temp[i] = x[xdata_counter];
            xdata_counter++;
        }
        
        xdot << First_Derivative(temp) << std::endl;
        std::cout << First_Derivative(temp) << std::endl;
    }
    
    
    xdata.close();
    xdot.close();
    std::cin.get();
    return 0;
}

and this is what I got

    x'
----------
184.15
34.0417
-147.475
-193.317
-61.45
126.942
198.608
87.675
-103.825
-199.917
317.333

Weird results!

[naraku9333]

while ( !xd.empty() )I don't see where you use that vector.

[CroCo]

@naraku9333,
I think you've posted after I changed it to
while ( xdata_counter < x.size() )

Now I guarantee that I will not exceed the size of my vector data.

[andywestken]

The MATLAB results you pasted above look a bit odd to me...

As the derivative of sin(x) is cos(x) it's easy enough to check the results: while cos(0.2) does equal 0.1987, as expected, sin(0.2) is 0.9801 not 0.1907.

So what the matlab output is showing is not the derivative. Going by the MATLAB help, diff just calculates the changes between values, which needs to be divided by the increment (0.2 in your case) to provide an approximation of the derivative. i.e. x' = diff(x)/0.2 (or whatever the MATLAB syntax is.)

See the comments (1#, 2#, 3#) I added to the code below. After I made these changes, and added std::setprecision(4) to fix the number of decimal places, I got:

0.9207
0.8255
0.6969
0.54
0.3621
0.1702
...

which is not quite cos(0.2), cos(0.4), cos(0.6), ... but is closer than you get using MATLAB's diff() divided by 0.2.

Edit: it looks like MATLAB should have been probably used like this:

t = 0:0.2:10
x = sin(t)
x' = diff(x)./diff(t)

(from example given here:
diff - Differences and approximate derivatives
http://www.mathworks.co.uk/help/matlab/ref/diff.html )

Andy

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <fstream>


double First_Derivative(double a[5])
{
    double firstDer(0.0);
    double t(0.001); // 1# Why is 0.001 used here when i/p data uses 0.2 ?
    
    // no loop here as 5 points is just enough to calculate the derivative
    // of one elem (the middle one). If you do use a loop, it has to start
    // at 2 and end 2 elems before the end (so i - 2 and i + 2 are always
    // valid indexes for the array.)
    
    size_t i(2);
    
    firstDer = (   4.0 / 3.0 * (a[i + 1] - a[i - 1]) / (2.0 * t)
                - 1.0 / 3.0 * (a[i + 2] - a[i - 2]) / (4.0 * t) );
    // Rather than use (double)4 you can just use 4.0. The .0 tells the
    // compiler you mean a double (if you need a float instead, append
    // an f, like 3.0f )
    // std::cout << firstDer << std::endl;
    return firstDer;
}


int main(void)
{
    std::ifstream xdata;
    std::ofstream xdot;
    double temp[5] = {0};
    std::vector<double> x, xd;
   
    
    xdata.open("/Users/bandarhd/Desktop/TestX.txt");
    xdot.open("/Users/bandarhd/Desktop/TestXdot.txt");
    
    if (xdata.fail())
    {
        std::cout << "Error in open xdata file";
    }
    if (xdot.fail())
    {
        std::cout << "Error in Opening xdot file";
    }
    
    double line;
    while ( !xdata.eof() )
    {
        xdata >> line;
        x.push_back(line);
        
    }
        
//    for (int i = 0; i < x.size(); ++i)
//    {
//        std::cout << x[i] << std::endl;
//    }
    
    int xdata_counter = 0;
    while ( xdata_counter < x.size() )
    {
        
        for (int i = 0; i < 5; ++i)
        {
            temp[i] = x[xdata_counter];
            xdata_counter++;
            // #2 Why use each point only once? Each point of input
            // should be used as the center point, and four other
            // times as the other elems in the temp array. Of course,
            // as you need points from -2 to 2, you can't calculate
            // the derivative for the first two or the last
            // two points in vector x, so while loop needs its
            // start and termination conditions changes. But then
            // you can calculate the deriv for each other point
            // by making it the middle one in temp.
        }
        
        // #3 why run the calc twice?
        xdot << First_Derivative(temp) << std::endl;
        std::cout << First_Derivative(temp) << std::endl;
    }
    
    
    xdata.close();
    xdot.close();
    std::cin.get();
    return 0;
}

[CroCo]

@andywestken,
Wow man, I really owe you. I did follow you suggestions and this is my modifications

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <fstream>


double First_Derivative(double a[5])
{
    double firstDer(0.0);
    double t(0.2); // 1# Why is 0.001 used here when i/p data uses 0.2 ?
    
    // no loop here as 5 points is just enough to calculate the derivative
    // of one elem (the middle one). If you do use a loop, it has to start
    // at 2 and end 2 elems before the end (so i - 2 and i + 2 are always
    // valid indexes for the array.)
    
    size_t i(2);
    
    firstDer = (   4.0 / 3.0 * (a[i + 1] - a[i - 1]) / (2.0 * t)
                - 1.0 / 3.0 * (a[i + 2] - a[i - 2]) / (4.0 * t) );
    // Rather than use (double)4 you can just use 4.0. The .0 tells the
    // compiler you mean a double (if you need a float instead, append
    // an f, like 3.0f )
    std::cout << firstDer << std::endl;
    return firstDer;
}


int main(void)
{
    std::ifstream xdata;
    std::ofstream xdot;
    double temp[5] = {0};
    std::vector<double> x, xd;
    
    
    xdata.open("/Users/bandarhd/Desktop/TestX.txt");
    xdot.open("/Users/bandarhd/Desktop/TestXdot.txt");
    
    if (xdata.fail())
    {
        std::cout << "Error in open xdata file";
    }
    if (xdot.fail())
    {
        std::cout << "Error in Opening xdot file";
    }
    
    double line;
    while ( !xdata.eof() )
    {
        xdata >> line;
        x.push_back(line);
        
    }
    
    int xdata_counter = 0;
    
    while ( xdata_counter < (x.size() - 2) )
    {
        temp[2] = x[xdata_counter];
        temp[3] = x[xdata_counter+1];
        temp[4] = x[xdata_counter+2];

        
        xdot << First_Derivative(temp) << std::endl;
        temp[1] = temp[2];
        temp[0] = temp[1];
        xdata_counter++;
    }
    
    
    xdata.close();
    xdot.close();
    std::cin.get();
    return 0;
}

and this is what I've got

0.500083
1.06275
0.92075
0.8255
0.696917
0.54
0.362083
0.170208
-0.0292083
-0.227292
-0.415958
-0.588375
-0.737375
-0.856917
-0.942208
-0.990083
-0.998042
-0.966583
-0.897
-0.790958
-0.653542
-0.49025
-0.30725
-0.112292
0.087625
0.283625
0.468292
0.634708
0.775625
0.885625
0.960083
0.996292
0.993042
0.950292
0.8695
0.754
0.60825
0.438375
0.251125
0.0541667
-0.145458
-0.339292
-0.519125
-0.67875
-0.811125
-0.911
-0.974667
-0.999583
-0.984792
-0.995333

which are perfect results. I'm dealing with a sensitive device. Is it possible to increase the accuracy of this method? Also, is there any remedy to smooth the beginning? In my case, the centre point of my temp data is the first element of my actual data but this is only at the beginning. This is why the beginning is odd. Also, I'm afraid if this could affect the stability of the system. And by the way, do yo suggest books in this field I mean C++ and math for real applications? Again thank you so much for your invaluable time.

By drawing it in Matlab, this is what I've got

http://s13.postimg.org/74gml77jr/Untitled.png

[andywestken]

1# results

which are perfect results

line 68 and 68 are the wrong way around: you are copying elem 2 over elem 1 and then copying the new elem 1 over elem 0. So elem 0 and 1 are always the same.

If you swap these two lines then you're results will be better.

#2 initial two and final two points

When is comes to the first two and last two points, you could special case them. That is, use one of the other formula discussed earlier in the thread.
- use the forward delta version of Alg #1 for the first point
- use the reverse delta version of Alg #1 for the last point
- use version of Alg #2 for the second and one before the last points

But to get improved accuracy, you will need to look for a way of interpolating which uses more than 2 or 3 points. (I've always just dropped the two first points.)

#3 Books

I am unsure about what you mean "maths for real applications". Some areas of computing, like games and fincancial applications, use quite a lot of maths. So books on these subjects should also be mathematical.

One book which might be relevant is "Numerical Recipes in C++: The Art of Scientific Computing".

Andy

[CroCo]

@Andy,
Thank you so much for being helpful and informative. Actually, the results that I've posted based on the following code which exactly what you suggested. I mistakenly did post the old code.

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <fstream>


double First_Derivative(double a[5])
{
    double firstDer(0.0);
    double t(0.2); // 1# Why is 0.001 used here when i/p data uses 0.2 ?
    
    // no loop here as 5 points is just enough to calculate the derivative
    // of one elem (the middle one). If you do use a loop, it has to start
    // at 2 and end 2 elems before the end (so i - 2 and i + 2 are always
    // valid indexes for the array.)
    
    size_t i(2);
    
    firstDer = (   4.0 / 3.0 * (a[i + 1] - a[i - 1]) / (2.0 * t)
                - 1.0 / 3.0 * (a[i + 2] - a[i - 2]) / (4.0 * t) );
    // Rather than use (double)4 you can just use 4.0. The .0 tells the
    // compiler you mean a double (if you need a float instead, append
    // an f, like 3.0f )
    std::cout << firstDer << std::endl;
    return firstDer;
}


int main(void)
{
    std::ifstream xdata;
    std::ofstream xdot;
    double temp[5] = {0};
    std::vector<double> x, xd;
    
    
    xdata.open("/Users/bandarhd/Desktop/TestX.txt");
    xdot.open("/Users/bandarhd/Desktop/TestXdot.txt");
    
    if (xdata.fail())
    {
        std::cout << "Error in open xdata file";
    }
    if (xdot.fail())
    {
        std::cout << "Error in Opening xdot file";
    }
    
    double line;
    while ( !xdata.eof() )
    {
        xdata >> line;
        x.push_back(line);
        
    }
    
    int xdata_counter = 0;
    
    while ( xdata_counter < (x.size() - 2) )
    {

        temp[2] = x[xdata_counter];
        temp[3] = x[xdata_counter+1];
        temp[4] = x[xdata_counter+2];
        
        xdot << First_Derivative(temp) << std::endl;
        temp[0] = temp[1];
        temp[1] = temp[2];
        
        xdata_counter++;
    }
    
    
    xdata.close();
    xdot.close();
    std::cin.get();
    return 0;
}

I'm still confused with delta. Now in my case the delta t = 0.2. Correct me if I'm wrong. In my case (sin and cos), for the first time calling derivative function, I need to set the delta t = 0.4e-5 (for example). Then after that I need to reset delta t = 0.2 ? Am I right? These are my modifications

long double First_Derivative(double a[5])
{
    long double firstDer(0.0);
    long double t(0.2); 
    static int delta_counter(0);
    
    // indicate to the first calling ( first points have been called )
    if ( delta_counter == 0 )
    {
        t = -0.4e2;
    }
    
    // indicate to the rest of points
    else 
    {
        t = 0.2;
    }
    
    size_t i(2);
    
    firstDer = (   4.0 / 3.0 * (a[i + 1] - a[i - 1]) / (2.0 * t)
                - 1.0 / 3.0 * (a[i + 2] - a[i - 2]) / (4.0 * t) );
    
    std::cout << firstDer << std::endl;
    delta_counter++;
    return firstDer;
}

But this didn't smooth the beginning. Also, the data table you've posted in the pervious post, why delta is not constant?

These what I got

-0.00250042 <- first element
1.06275
0.92075
0.8255
0.696917
0.54
.
.
.

[andywestken]

I need to set the delta t = 0.4e-5 (for example). Then after that I need to reset delta t = 0.2 ?

Why do you think you need to change the sampling period? Isn't this fixed??

What changes it the number of points available for the calculation. As you do not have two preceeding points available when you calculate the first point, you cannot use

Alg#3 = 4/3 (f(x_n+1) - f(x_n-1)) / 2t - 1/3 (f(x_n+2) - f(x_n-2)) / 4t

You have to use the n/n+1 version of

Alg#1 = ( f(x_n) - f(x_n+1) ) / t

Similarly for the last point, except then you have to use the n-1/n form.

For the second and one-before-the-last points you can use

Alg#2 = ( f(x_n+1) - f(x_n-1) ) / 2t

If this isn't enough, you will have to look for a formula which works with more of the available points.

Andy

[CroCo]

@Andy,

Thank so much again. You are very helpful and informative person. I'm going now to test it with the device.

[CroCo]

Hello Andy,

These what I've got after implementing your method.

This is x_position

http://s1.postimg.org/f2jrvay3j/xposition.png

This is the results of the velocities by using Matlab

http://s18.postimg.org/qv2xetxcp/vxmatlab.png

This is the results of the velocities by using your method

http://s17.postimg.org/dsh1hof0f/image.png

[andywestken]

Well, from what I can see, they look quite similar? So, are you happy with the results??

But it's not my method, it's from Karsten Ahnert and Markus Abel's paper!

Andy

PS For the record:

Computer Physics Communications 177 (2007) 764–774
Numerical differentiation of experimental data: local versus global methods
http://www.stat.physik.uni-potsdam.de/~kahnert/publications/COMPHY3395.pdf

[CroCo]

@Andy,
Yes, I'm very happy with these results. Thank you so much again. In programming, I'm facing new problem. I'm dealing with a lot of variables.
x, y, z & x', y', z' & x'', y'', z'' These are only for one device. The problem is how can I deal with velocities the way I did with positions? This is the code

  while (1)
    {
        /* Perform a synchronous call to copy the most current device state.
           This synchronous scheduler call ensures that the device state
           is obtained in a thread-safe manner. */
        hdScheduleSynchronous(copyDeviceDataCallback, &currentData, HD_MIN_SCHEDULER_PRIORITY);           
		
		
		x_Position[2] =  currentData.m_devicePosition[0];
		x_Position[3] =  currentData.m_devicePosition[0];
		x_Position[4] =  currentData.m_devicePosition[0];

		y_Position[2] =  currentData.m_devicePosition[1];
		y_Position[3] =  currentData.m_devicePosition[1];
		y_Position[4] =  currentData.m_devicePosition[1];

		z_Position[2] =  currentData.m_devicePosition[2];
		z_Position[3] =  currentData.m_devicePosition[2];
		z_Position[4] =  currentData.m_devicePosition[2];

		std::cout << First_Derivative(x_Position) << " " << First_Derivative(y_Position) << " " << First_Derivative(z_Position) << std::endl;

		x_Position[0] =  x_Position[1];
		x_Position[1] =  x_Position[2];

		y_Position[0] =  y_Position[1];
		y_Position[1] =  y_Position[2];

		z_Position[0] =  z_Position[1];
		z_Position[1] =  z_Position[2];

	}

I've tried to create x_Velocity[5] and

	
                x_Position[2] =  currentData.m_devicePosition[0];
		x_Position[3] =  currentData.m_devicePosition[0];
		x_Position[4] =  currentData.m_devicePosition[0];

                x_Velocity[2] = First_Derivative(x_Position);
		x_Velocity[3] = First_Derivative(x_Position);
		x_Velocity[4] = First_Derivative(x_Position);

		
		x_Velocity[1] = x_Velocity[2];
		x_Velocity[0] = x_Velocity[1];

		x_Position[0] =  x_Position[1];  // <----- is affected by calling First_Derivative()
		x_Position[1] =  x_Position[2];  // <----- is affected by calling First_Derivative()

I'm trying also to store x_Velocity in vector. Every time I deal with whether array or vector, the first two points of x position are affected because I'm calling the function for the velocity. Any suggestions?

[andywestken]

If it's a different problem, it would better if you shifted this new question to a new thread.

It's late where I am, so my brain has shut down for the day. If no one esle has responded, I'll look out for you new thread tomorrow some time.

Andy

[CroCo]

@Andy,
Thank you again. I really appreciate your help.

I will post a new post for that.