Derivatives in C/C++?

Question

I have some expressions such as x^2+y^2 that I'd like to use for some math calculations. One of the things I'd like to do is to take partial derivatives of the expressions.

So if f(x,y) = x^2 + y^2 then the partial of f with respect to x would be 2x, the partial with respect to y would be 2y.

I wrote a dinky function using a finite differences method but I'm running into lots of problems with floating point precision. For example, I end up with 1.99234 instead of 2. Are there any libraries that support symbolic differentiation? Any other suggestions?

Solution 1

I've implemented such libraries in several different languages, but unfortunately not C. If you are dealing only with polynomials (sums, products, variables, constants, and powers) it is pretty easy. Trig functions are also not too bad. Anything more complicated and you will probably be better off taking the time to master somebody else's library.

If you decide to roll your own, I have a few suggestions that will simplify your life:

• Use immutable data structures (purely functional data structures) to represent expressions.

• Use Hans Boehm's garbage collector to manage the memory for you.

• To represent a linear sum, use a finite map (e.g., a binary search tree) to map each variable to its coefficient.

If you're willing to embed Lua into your C code and do your computations there, I have put my Lua code at http://www.cs.tufts.edu/~nr/drop/lua. One of the nicer features is that it can take a symbolic expression, differentiate it, and compile the results into Lua. You will of course find no documentation whatever :-(

Solution 2

If you're doing numerical differentiation ("evaluate the derivative of f(x) at x = x0") and you know you're equations in advance (ie, not user input), then I'd recommend FADBAD++. It's a C++ template library for solving numeric derivatives using Automatic differentiation. It's very quick, and accurate.

Solution 3

Getting numerical differentiation "right" (in the sense of minimizing errors) can be quite tricky. To get started, you may want to take a look at the Numerical Recipies' section on numerical derivatives.

For free symbolic math packages, you should look at GiNaC. You could also check out SymPy, a self-contained, pure-python symbolic math package. You will find that SymPy is much easier to explore, since you can use it interactively from the Python command line.

On the commercial end, both Mathematica and Maple have C APIs. You need an installed/licensed version of the program to use the libraries, but both can easily do the type of symbolic differentiation you're after.

Solution 4

You can improve the accuracy of your numerical differentiation in two simple ways

1. Use a smaller delta. You appear to have used a value of around 1e-2. Start with 1e-8, and test if getting any smaller hurts or helps. Obviously you can't get too close to the machine precision - about 1e-16 for double.

2. Use central differences rather than forward (or backwards) differences. i.e. df_dx =(f(x+delta) - f(x-delta)) / (2.0*delta) For reasons to do with cancellation of higher truncation terms, the error in the central differences estimate is of the order delta^2 rather than the delta of forward differences. See http://en.wikipedia.org/wiki/Finite_difference

Solution 5

Sorry for bringing this up after 6 years. However, I was looking for such a library to my project and I've seen @eduffy suggests FADBAD++. I've read the documentation and came back to your question. I feel my answer will be beneficial, therefore, the following code is for your case.

#include <iostream>
#include "fadiff.h"

using namespace fadbad;

F<double> func(const F<double>& x, const F<double>& y)
{
    return x*x + y*y;
}

int main()
{
    F<double> x,y,f;     // Declare variables x,y,f
    x=1;                 // Initialize variable x
    x.diff(0,2);         // Differentiate with respect to x (index 0 of 2)
    y=1;                 // Initialize variable y
    y.diff(1,2);         // Differentiate with respect to y (index 1 of 2)
    f=func(x,y);         // Evaluate function and derivatives

    double fval=f.x();   // Value of function
    double dfdx=f.d(0);  // Value of df/dx (index 0 of 2)
    double dfdy=f.d(1);  // Value of df/dy (index 1 of 2)

    std::cout << "    f(x,y) = " << fval << std::endl;
    std::cout << "df/dx(x,y) = " << dfdx << std::endl;
    std::cout << "df/dy(x,y) = " << dfdy << std::endl;

    return 0;
}

The output is

    f(x,y) = 2
df/dx(x,y) = 2
df/dy(x,y) = 2

Another example, let's say that we are interested in the first derivative of sin(). Analytically, it is cos. This is great because we need to compare the true derivative of a given function and its numerical counterpart to compute the true error.

#include <iostream>
#include "fadiff.h"

using namespace fadbad;

F<double> func(const F<double>& x)
{
    return sin(x);
}



int main()
{
    F<double> f,x;
    double dfdx;
    x = 0.0;
    x.diff(0,1);
    f = func(x);
    dfdx=f.d(0);


    for (int i(0); i < 8; ++i ){
        std::cout << "       x: " << x.val()        << "\n"
                  << "    f(x): " << f.x()          << "\n" 
                  << " fadDfdx: " << dfdx           << "\n"
                  << "trueDfdx: " << cos(x.val())   << std::endl;
        std::cout << "=========================="   << std::endl;

        x += 0.1;
        f = func(x);
        dfdx=f.d(0);
    }


    return 0;
}

The result is

       x: 0
    f(x): 0
 fadDfdx: 1
trueDfdx: 1
==========================
       x: 0.1
    f(x): 0.0998334
 fadDfdx: 0.995004
trueDfdx: 0.995004
==========================
       x: 0.2
    f(x): 0.198669
 fadDfdx: 0.980067
trueDfdx: 0.980067
==========================
       x: 0.3
    f(x): 0.29552
 fadDfdx: 0.955336
trueDfdx: 0.955336
==========================
       x: 0.4
    f(x): 0.389418
 fadDfdx: 0.921061
trueDfdx: 0.921061
==========================
       x: 0.5
    f(x): 0.479426
 fadDfdx: 0.877583
trueDfdx: 0.877583
==========================
       x: 0.6
    f(x): 0.564642
 fadDfdx: 0.825336
trueDfdx: 0.825336
==========================
       x: 0.7
    f(x): 0.644218
 fadDfdx: 0.764842
trueDfdx: 0.764842
==========================

View more solutions

Comments

[Jules]

Exactly. Symbolic math would be a nightmare in C. Just getting the mallocs and frees right would take longer than using a high level functional language.

[nlucaroni]

I initially didn't even think of that (too used to using functional languages with GCs I guess). But you are additionally right. It would be a nightmare.

[alfC]

The core of Mathematica is written in C. So, it is possible.

[Admin]

It is 4 years after I've heard first about boost-proto and I'm still unclear about what value is provided in boost-proto.

[Ankit Roy]

@user678269: 2 reasons: (1) expressiveness. You can use it to easily create an AST data structures, e.g. one representing the expression x * x + y * y just by typing out that C++ expression verbatim in your code; if you have written code for manipulating ASTs (e.g. for evaluating, compiling, pretty-printing or (for mathematical expressions) differentiating/integrating/expanding/simplifying them you can then apply this. ...

[Ankit Roy]

... (2) efficiency. If you just want to evaluate the expression, the usual C++ computations will create and destroy many temporaries unnecessarily; you can often do better by building and immediately evaluating an AST, which behaves a bit like a custom-built constructor that only performs one memory allocation (see the first google result for more details).

[CroCo]

Thank you for suggesting this library. I've double checked and posted an answer for the OP. It seems very good library. Thanks again.