2009-01-07 93 views
21

我有一些表達式,例如x^2+y^2,我想用於一些數學計算。我想要做的一件事就是採用表達式的偏導數。C/C++中的衍生物?

所以如果f(x,y) = x^2 + y^2則部分的f相對於x2x,部分相對於y2y

我寫了一個使用有限差分法的dinky函數,但我遇到了很多浮點精度問題。例如,我最終以1.99234而不是2。有沒有支持符號區分的圖書館?還有其他建議嗎?

+0

滾動你自己將是一個壞主意。我在「計算機代數」中選擇了一門交叉列出的畢業課程,而且這些課程的功能還不完整。 簡單的區分可以通過應用maty規則來完成,然後嘗試用C/C++替代和評估...很難做到。 – Calyth 2009-01-08 00:42:11

+0

另請參閱http://stackoverflow.com/questions/627055/compute-a-derivative-using-discrete-methods,也許標記你的問題的衍生物 – 2009-07-22 07:47:12

回答

5

獲得數值區分「正確」(從最小化錯誤的意義上說)可能非常棘手。要開始,您可能需要查看numerical derivatives上的「數字收錄」部分。

對於免費的符號數學包,你應該看看GiNaC。您還可以查看SymPy,這是一款自包含的純python符號數學軟件包。您會發現SymPy更容易探索,因爲您可以從Python命令行以交互方式使用SymPy。

在商業上,Mathematica和Maple都有C API。您需要安裝/許可的程序版本才能使用這些庫,但兩者都可以輕鬆完成您所追求的符號差異類型。

0

由於它適用於Lisp而不適用於C/C++,但它可以幫助其他人尋找類似的任務,或者您可能會自行獲得有關在C/C++中實現類似功能的一些想法。 SICP有對口齒不清的問題的一些講座:

  1. 衍生規則3b
  2. 代數規則4a

在Lisp中,這是很簡單的(並與強大的模式匹配和多態等功能性語言類型)。在C語言中,您可能不得不大量使用枚舉和結構來獲得相同的權力(更不用說分配/解除分配)。人們可以在一小時內明確地編寫你需要的東西 - 我認爲打字速度是限制因素。如果你需要C,你實際上可以從C中調用ocaml(反之亦然)。

+0

準確。符號數學將是C中的噩夢只需要使用高級功能語言就可以花費更長的時間。 – Jules 2009-01-08 00:09:29

+0

我起初甚至沒有想到(我習慣使用功能語言和GC)。但你也是對的。這將是一場噩夢。 – nlucaroni 2009-01-08 00:38:00

11

我已經用幾種不同的語言實現了這樣的庫,但不幸的是不是C.如果你只處理多項式(和,產品,變量,常量和冪),這很容易。 Trig函數也不錯。任何更復雜的事情,你可能會花更多的時間去掌握別人的圖書館。

如果您決定推出自己的,我有幾個建議,可以簡化你的生活:

  • 使用一成不變的數據結構(純功能性數據結構)來表示表達式。使用Hans Boehm's garbage collector爲您管理內存。

  • 爲了表示線性和,使用有限圖(例如,,二叉搜索樹)將每個變量映射到它的係數。

如果你願意嵌入Lua到你的C代碼,做你的計算有,我已經把我的Lua代碼在http://www.cs.tufts.edu/~nr/drop/lua。其中一個更好的特點是它可以採用符號表達,區分它,並將結果編譯成Lua。你當然會發現沒有任何文檔:-(

1

如果這真的是你想要使用的函數,編寫一個類庫會很容易,從一個單詞開始,有一個係數和一個指數有一個多項式將包含一個術語集合

如果您爲感興趣的數學方法(例如add/sub/mul/div/differentiate/integrate)定義接口,那麼您正在尋找在GoF的複合模式。這兩種期限和多項式將實現該接口。多項式只會在其收集各期限迭代。

1

它肯定會更容易利用現有的包比編寫自己的,但如果你決定寫自己的,並且你準備花一些時間學習一些C++的黑暗角落,你可以使用BoostBoost.Proto來設計你自己的庫。

基本上,Boost.Proto允許你任何有效的C++表達式轉換,如x * x + y * yexpression template - 即表達的解析樹的基本上是一個表示使用嵌套struct秒 - 然後經執行任意的計算稍後通過調用proto::eval()解析樹。默認情況下,proto::eval()用於評估樹,就好像它已經直接運行一樣,儘管沒有理由不修改每個函數或運算符的行爲來代替採用符號導數。

雖然這將是一個非常複雜的解決方案,但它會比使用C++模板元編程技術展示自己的表達式模板容易得多。

5

如果你正在做數值微分(「評估F(X)在X = X0的衍生物」),你知道你預先方程(即,不是用戶輸入),那麼我建議FADBAD++。它是一個C++模板庫,用於解決使用Automatic differentiation的數字衍生產品。它非常快速,準確。

0

只計算一階導數是相當簡單的實現。 但它是一種讓藝術變得快速的藝術。 你需要一個類,其中包含

  • 階導數值與獨立變量

的陣列然後將寫爲加法和減法運算符等和功能等罪()實現了這個操作的基本和衆所周知的規則。

要計算更高階的衍生物應該使用截斷泰勒系列。 您也可以將上面提到的類應用於自身 - 值和派生值的類型應該是模板參數。 但這意味着不止一次計算和存儲衍生物。

截斷泰勒 - 有兩個庫可用於此:

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

http://www.foelsche.com/ctaylor

4

您可以提高數值微分的準確性,兩種簡單的方法

  1. 使用較小的增量。您似乎使用了大約1e-2的值。從1e-8開始,並測試是否獲得更小的傷害或幫助。顯然你不能太接近機器的精度 - 約爲1e-16雙倍。

  2. 使用中心差異而不是前向(或後向)差異。 即df_dx =(f(x+delta) - f(x-delta))/(2.0*delta) 由於與取消較高截斷項有關的原因,中心差異估計的誤差爲delta^2的順序,而不是前向差的增量。請參見http://en.wikipedia.org/wiki/Finite_difference

0

看看Theano,它支持符號分化(在神經網絡中)。該項目是開源的,所以你應該能夠看到他們是如何做到的。

0

對不起,6年後提起這件事。但是,我正在爲我的項目尋找這樣一個圖書館,我看到@eduffy建議FADBAD++。我已閱讀文檔並回到您的問題。我覺得我的回答會有好處,因此,下面的代碼適合您的情況。

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

輸出是

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

另一個例子,讓我們說,我們感興趣的sin()一階導數。分析上,它是cos。這很好,因爲我們需要比較一個給定函數的真實導數和它的數值副本來計算真實誤差。

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

結果是

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