最近項(xiàng)目里總是要對(duì)很龐大的公式求導(dǎo)，很煩人，手工求導(dǎo)容易出錯(cuò)。
當(dāng)然MATLAB是個(gè)好選擇，不過當(dāng)它要錢的時(shí)候，您可能就不這么認(rèn)為了。
于是，實(shí)現(xiàn)了一個(gè)可以
編譯期求導(dǎo)(不用擔(dān)心運(yùn)行時(shí)負(fù)擔(dān))的小型庫(kù)，還不完全，僅支持多項(xiàng)式，sin，cos，pow，exp，log等函數(shù)求導(dǎo)。
后期的表達(dá)式優(yōu)化做的不是很好。
下面是一些測(cè)試代碼，完整的源碼在
http://www.boostpro.com/vault/index.php?action=downloadfile&filename=[math]AD.zip實(shí)現(xiàn)部分很復(fù)雜，請(qǐng)多多指教。
只有1個(gè)函數(shù)， d(...)
支持高階，多元求導(dǎo)。
d(exp, var)(value1, value2, ...)
exp內(nèi)可以有多個(gè)變量，var表示要對(duì)其求導(dǎo)的變量，value表示求導(dǎo)以后用于計(jì)算表達(dá)式的變量的值。
比如:
d(d(x*x*x, x),x)(3.0) 表示對(duì)x*x*x求二階導(dǎo)數(shù)在x=3.0時(shí)候的值。
d(d(x*x*y, x), y)(3.0, 4.0) 表示d(x*x*y)/(dxdy)在x=3.0，y=4.0的值。
d(d(x*x*x, x) +d(y*x, y), y) (2.0) 則表示 (d(x*x*x)/dx + d(y*x)/dy)/dy == 0。
可以直接用cout把求導(dǎo)后的表達(dá)式輸出，不用給變量給值。
cout<<d(x*x, x) // 結(jié)果是：2*x
這里沒有用任何迭代，是直接對(duì)表達(dá)式求導(dǎo)的。返回值是求導(dǎo)后的表達(dá)式，本質(zhì)是一個(gè)仿函數(shù)。可以用boost::function保存起來使用。
例如：
boost::function<double (double)> df = d(pow(x, const_<10>::type()), x); //df 參數(shù)為1個(gè)double，返回double
然后就可以在任何地方使用 df 了：
double res = df(3.0) // res == pow(3, 9)
1
#include "ad.h"
2#include <iostream>
3#include <iterator>
4
5using namespace std;
6
7int main()
8
{
9
variable<0>::type x;
10 variable<1>::type y;
11
12 double res[14];
13
14 res[0] = d(pow(x, const_<10>::type()), x)(2.0);
15
16 res[1] = d(x * x * x, x)(2.0);
17 res[2] = d(x + x + x, x)(2.0);
18 res[3] = d(x - x - x, x)(2.0);
19 res[4] = d(x / x, x)(2.0);
20
21 res[5] = d(pow(x, var(3.0)), x)(2.0);
22 res[6] = d(pow(var(3.0), x), x)(2.0);
23 res[7] = d(pow(x, x), x)(2.0);
24
25 res[8] = d(log(x), x)(2.0);
26 res[9] = d(exp(x), x)(2.0);
27
28 res[10] = d(sin(x), x)(2.0);
29 res[11] = d(cos(x), x)(2.0);
30
31 res[12] = d(d(sin(x) * cos(y), x), y)(2.0, 3.0);
32
33 res[13] = (d(log(x) + x, x) * x)(2.0);
34
35 copy(res, res + 14, ostream_iterator<double>(cout, "\n"));
36
37 cout<<d(pow(x, const_<10>::type()), x)<<endl;
38
39 cout<<d(x * x * x, x)<<endl;
40 cout<<d(x + x + x, x)<<endl;
41 cout<<d(x - x - x, x)<<endl;
42 cout<<d(x / x / x, x)<<endl;
43
44 cout<<d(pow(x, var(3.0)), x)<<endl;
45 cout<<d(pow(var(3.0), x), x)<<endl;
46 cout<<d(pow(x, x), x)<<endl;
47
48 cout<<d(log(x), x)<<endl;
49 cout<<d(exp(x), x)<<endl;
50
51 cout<<d(sin(x), x)<<endl;
52 cout<<d(cos(x), x)<<endl;
53
54 cout<<d(d(sin(x) * cos(y), x), y)<<endl;
55
56 cout<<(d(log(x) + x, x) * x)<<endl;
57
58 return 0;
59
}
60
輸出結(jié)果如下：
1512
212
33
4-1
50
612
79.88751
86.77259
90.5
107.38906
11-0.416147
120.909297
13-0.0587266
143
15pow(x,9)
16(((x+x)*x)+(x*x))
173
18-1
19(-1/(x*x))
20(pow(x,3)*(3*(1/x)))
21(pow(3,x)*log(3))
22(pow(x,x)*(log(x)+1))
23(1/x)
24exp(x)
25cos(x)
26sin(x)
27(cos(x)*sin(y))
28(((1/x)+1)*x)
29
posted on 2009-05-01 23:50
尹東斐 閱讀(2625)
評(píng)論(6) 編輯 收藏 引用