我有一些表达式，例如x^2+y^2我想用于一些数学计算的表达式。我想做的一件事是取表达式的偏导数。
因此，如果f(x,y) = x^2 + y^2then 的部分f相对于x将是2x，则相对于的部分y将是2y。
我使用有限差分法编写了一个极小的函数，但我遇到了很多浮点精度问题。例如，我最终得到1.99234而不是2. 有没有支持符号微分的库？还有其他建议吗？
idontwanttortfm
问问题
34369 次
10 回答
13
我已经用几种不同的语言实现了这样的库，但不幸的是不是 C。如果你只处理多项式（和、乘积、变量、常数和幂），这很容易。三角函数也不错。任何更复杂的事情，你可能最好花时间掌握别人的图书馆。
如果您决定自己动手，我有一些建议可以简化您的生活：
使用不可变数据结构（纯函数数据结构）来表示表达式。
使用Hans Boehm 的垃圾收集器为您管理内存。
为了表示线性和，使用有限映射（例如，二叉搜索树）将每个变量映射到其系数。
如果您愿意将Lua嵌入到您的 C 代码中并在那里进行计算，我已将 Lua 代码放在http://www.cs.tufts.edu/~nr/drop/lua上。一个更好的特性是它可以接受一个符号表达式，区分它，并将结果编译到 Lua 中。您当然不会找到任何文档:-(
于 2009-01-07T23:57:05.530 回答
7
如果您正在进行数值微分（“评估f(x)在x = x0处的导数”）并且您事先知道自己是方程式（即，不是用户输入），那么我建议您使用FADBAD++。它是一个 C++ 模板库，用于使用自动微分求解数值导数。它非常快速，准确。
于 2009-06-12T12:42:00.200 回答
4
您可以通过两种简单的方式提高数值微分的准确性
使用较小的增量。您似乎使用了 around 的值
1e-2。从 开始1e-8，并测试是否有任何较小的伤害或帮助。显然你不能太接近机器精度——大约1e-16是两倍。使用中心差异而不是前向（或后向）差异。即
df_dx =(f(x+delta) - f(x-delta)) / (2.0*delta)由于取消较高截断项的原因，中心差异估计中的误差是阶数delta^2而不是前向差异的增量。见http://en.wikipedia.org/wiki/Finite_difference
于 2013-03-12T21:51:35.523 回答
1
如果这确实是您想要使用的那种函数，那么编写一个类库就很容易了。从一个 Term 开始，带有一个系数和一个指数。有一个由一组项组成的多项式。
如果您为感兴趣的数学方法定义了一个接口（例如，add/sub/mul/div/differentiate/integrate），那么您正在查看一个 GoF 复合模式。项和多项式都将实现该接口。多项式将简单地迭代其集合中的每个项。
于 2009-01-08T03:19:01.003 回答
1
利用现有的包肯定比编写自己的包更容易，但如果你决心编写自己的包，并且准备花一些时间了解 C++ 的一些黑暗角落，你可以使用Boost.Proto从Boost到设计您自己的库。
基本上，Boost.Proto 允许您将任何有效的 C++ 表达式转换为x * x + y * y表达式模板（基本上是使用嵌套structs 表示该表达式的解析树），然后稍后对该解析树执行任何任意计算时间通过调用proto::eval()它。默认情况下，proto::eval()用于评估树，就好像它已经直接运行一样，尽管没有理由不能修改每个函数或运算符的行为以取而代之的是符号导数。
尽管这对于您的问题来说是一个极其复杂的解决方案，但它仍然比尝试使用 C++ 模板元编程技术滚动您自己的表达式模板要容易得多。
于 2009-01-08T04:52:50.427 回答
1
很抱歉在 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
==========================
于 2015-12-02T08:30:12.527 回答
0
这是一种旁白，因为它适用于 Lisp 而不是 C/C++，但它可以帮助其他人寻找类似的任务，或者您可能会获得一些想法，以自己在 C/C++ 中实现类似的东西。SICP 有一些关于 lisp 主题的讲座：
在 Lisp 中，它非常简单（在其他具有强大模式匹配和多态类型的函数式语言中）。在 C 中，您可能不得不大量使用枚举和结构来获得相同的功能（更不用说分配/解除分配）。一个小时内就可以明确地在 ocaml 中编写你需要的代码——我想说打字速度是限制因素。如果你需要 C，你实际上可以从 C 调用 ocaml（反之亦然）。
于 2009-01-07T23:46:11.247 回答
0
看看 Theano，它支持符号微分（在神经网络的上下文中）。该项目是开源的，所以你应该能够看到他们是如何做到的。
于 2015-03-06T03:06:36.820 回答
0
只计算一阶导数是很容易实现的。但它是一门艺术，使其快速。你需要一个类，其中包含
- 价值
- 一组导数值与自变量
然后，您将编写用于加法和减法等的运算符以及诸如 sin() 之类的函数，它们实现了此操作的基本和众所周知的规则。
要计算高阶导数，应使用截断的泰勒级数。您也可以将上述类应用于自身——值和派生值的类型应该是模板参数。但这意味着不止一次地计算和存储衍生品。
截断的 taylor 系列——有两个库可用于此：
于 2012-06-01T16:15:53.257 回答