Derivados en C/C++?

Question

Tengo algunas expresiones como x^2+y^2 que me gustaría usar para algunos cálculos matemáticos. Una de las cosas que me gustaría hacer es tomar derivados parciales de las expresiones.

Así que si entonces el f(x,y) = x^2 + y^2 parcial de f con respecto al x habría 2x, el parcial con respecto a y habría 2y.

Escribí una función dinky usando un método de diferencias finitas, pero estoy teniendo muchos problemas con la precisión de punto flotante. Por ejemplo, termino con 1.99234 en lugar de 2. ¿Hay bibliotecas que respalden la diferenciación simbólica? ¿Cualquier otra sugerencia?

Answer 1

Obtener la diferenciación numérica "a la derecha" (en el sentido de minimizar los errores) puede ser bastante complicado. Para comenzar, es posible que desee echar un vistazo a la sección de Recetas Numéricas en numerical derivatives.

Para paquetes matemáticos simbólicos gratuitos, debe consultar GiNaC. También puede consultar SymPy, un paquete matemático simbólico autónomo de pura pitón. Descubrirá que SymPy es mucho más fácil de explorar, ya que puede usarlo de forma interactiva desde la línea de comando de Python.

En el extremo comercial, tanto Mathematica como Maple tienen C API. Necesita una versión instalada/licenciada del programa para usar las bibliotecas, pero ambas pueden hacer fácilmente el tipo de diferenciación simbólica que busca.

Answer 2

Esto es algo aparte ya que se aplica a Lisp y no a C/C++, pero podría ayudar a otros a buscar tareas similares o podría obtener algunas ideas sobre cómo implementar algo similar en C/C++ por su cuenta. SICP tiene algunas conferencias sobre el tema de Lisp:

1. reglas derivadas 3b
2. reglas algebraicas 4a

En Lisp, es bastante sencillo (y en otros lenguajes funcionales con una potente coincidencia de patrones y polimórfica tipos). En C, es probable que tenga que utilizar enums y estructuras para obtener el mismo poder (sin mencionar asignación/desasignación). Uno podría codificar definitivamente lo que necesita en ocaml en menos de una hora - Diría que la velocidad de tipeo es el factor limitante. Si necesita C, puede llamar a ocaml desde C (y viceversa).

Answer 3

Implementé estas bibliotecas en varios idiomas, pero desafortunadamente no C. Si solo se trata de polinomios (sumas, productos, variables, constantes y potencias), es bastante fácil. Las funciones Trig tampoco son tan malas. Algo más complicado y probablemente será mejor que se tome el tiempo para dominar la biblioteca de otra persona.

Si decide rodar su propia, tengo algunas sugerencias que simplificarán su vida:

• utilizar las estructuras de datos inmutables (estructuras de datos puramente funcionales) para representar expresiones.

• Utilice Hans Boehm's garbage collector para administrar la memoria para usted.

• Para representar una suma lineal, utilice un mapa finito (p., un árbol de búsqueda binario) para asignar cada variable a su coeficiente.

Si usted está dispuesto a incrustar Lua en el código C y hace sus cálculos allí, he puesto mi código Lua en http://www.cs.tufts.edu/~nr/drop/lua. Una de las características más agradables es que puede tomar una expresión simbólica, diferenciarla y compilar los resultados en Lua. Por supuesto, no encontrará documentación :-(

Answer 4

Si este es realmente el tipo de función que desea usar, sería bastante fácil escribir una biblioteca de clases. Comience con un solo Término, con un coeficiente y un exponente. Tener un polinomio que consistiría en una colección de términos.

Si define una interfaz para los métodos matemáticos de interés (por ejemplo, add/sub/mul/div/differentiate/integrate), usted está buscando en un patrón Compuesto GoF. Tanto Término como Polinomio implementarían esa interfaz. El Polinomio simplemente iteraría sobre cada Término en su Colección.

Answer 5

Sin duda sería más fácil aprovechar un paquete existente que para escribir el suyo, pero si está decidido a escribir el suyo, y está dispuesto a pasar algún tiempo aprendiendo sobre algunos rincones oscuros de C++, puede usar el Boost.Proto del Boost para diseñar su propia biblioteca.

Básicamente, Boost.Proto le permite convertir cualquier expresión válida de C++, como x * x + y * y a un expression template - básicamente una representación del árbol de análisis sintáctico de esa expresión utilizando anidados struct s - y luego realizar cualquier cálculo arbitrario sobre que analizar el árbol en otro momento llamando al proto::eval() en él. Por defecto, proto::eval() se usa para evaluar el árbol como si se hubiera ejecutado directamente, aunque no hay ninguna razón por la que no se pueda modificar el comportamiento de cada función u operador para tomar una derivada simbólica.

Aunque esta sería una solución compleja de extremadamente a su problema, sería mucho más fácil que intentar rodar sus propias plantillas de expresión usando técnicas de metaprogramación de plantillas C++.

Answer 6

Si estás haciendo la diferenciación numérica ("evaluar la derivada de f (x) en x = x0 ") y usted sabe que está ecuaciones de antelación (es decir, no la entrada del usuario), entonces yo recomendaría FADBAD++. Es una biblioteca de plantillas C++ para resolver derivados numéricos usando Automatic differentiation. Es muy rápido y preciso.

Answer 7

Para calcular solo el primer orden de derivados es bastante trivial de implementar. Pero es un arte hacerlo rápido. Se necesita una clase, que contiene

• el valor
• una matriz de valores derivados vs. variables independientes

A continuación, escribir los operadores de suma y resta y así sucesivamente y funciones como el pecado() que implementan las reglas básicas y conocidas para esta operación.

Para calcular las derivadas de orden superior se debe hacer utilizando series de taylor truncadas. También puede aplicar la clase mencionada anteriormente a sí mismo: el tipo para el valor y los valores derivados debe ser un argumento de plantilla. Pero esto significa el cálculo y el almacenamiento de derivados más de una vez.

truncada serie de Taylor - hay dos bibliotecas disponibles para esto:

http://code.google.com/p/libtaylor/

http://www.foelsche.com/ctaylor

Answer 8

se puede mejorar la exactitud de su diferenciación numérica de dos maneras simples

1. Usa un delta más pequeño. Parece que usaste un valor de alrededor de 1e-2. Comience con 1e-8 y pruebe si le duele algo más pequeño o lo ayuda. Obviamente, no puede acercarse demasiado a la precisión de la máquina - aproximadamente 1e-16 para el doble.

2. Use las diferencias centrales en lugar de las diferencias hacia delante (o hacia atrás). es decir, df_dx =(f(x+delta) - f(x-delta))/(2.0*delta) Por razones relacionadas con la cancelación de términos de truncamiento más altos, el error en la estimación de diferencias centrales es del orden delta^2 en lugar del delta de diferencias forward. Ver http://en.wikipedia.org/wiki/Finite_difference

Answer 9

Tenga una mirada en Teano, es compatible con la diferenciación simbólica (en el contexto de redes neuronales). El proyecto es de código abierto, por lo que debería poder ver cómo lo hacen.

Answer 10

Perdón por mencionar esto después de 6 años. Sin embargo, estaba buscando una biblioteca de este tipo para mi proyecto y he visto que @eduffy sugiere FADBAD++. He leído la documentación y volví a tu pregunta. Siento que mi respuesta será beneficiosa, por lo tanto, el siguiente código es para su caso.

#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; 
} 

La salida es

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

Otro ejemplo, digamos que estamos interesados ​​en la primera derivada de sin(). Analíticamente, es cos. Esto es genial porque necesitamos comparar la derivada verdadera de una función dada y su contraparte numérica para calcular el error verdadero.

#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; 
} 

El resultado es

 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 
========================== 

Answer 11

creé tal biblioteca en C++, se puede ver la muestra aquí: https://github.com/MartinKosicky/auto_diff/blob/master/sample/main.cpp

El uso es bastante simple, que también es compatible con matrices. Lo utilicé para crear una Red Neuronal Recurrente ... Tendría que agregar la multiplicación de la matriz GPU a través de