c++ - 在C++中计算雅可比矩阵(符号数学)

介绍

假设以下ODE集需要Jacobian矩阵：

dxdt[ 0 ] = -90.0 * x[0] - 50.0 * x[1];
dxdt[ 1 ] = x[0] + 3*x[1];
dxdt[ 2 ] = x[1] + 50*x[2];

在Matlab / Octave中，这将非常简单：

syms x0 x1 x2;
f = [-90.0*x0-50.0*x1, x0+3*x1, x1+50*x2]
v=[x0, x1, x2]
fp = jacobian(f,v)

这将产生以下输出矩阵：

[-90  -50  0 ]
[ 1    3   0 ]
[ 0    1   50]

我需要的

现在，我想在C ++中重现相同的结果。我无法先计算雅可比行列式并对其进行硬编码，因为它取决于例如用户输入和时间。所以我的问题是：该怎么做？通常对于数学运算，我使用Boost库，但是在这种情况下，我找不到任何解决方案。 implicit systems中对此只有简短的注释，但是以下代码不起作用：

sys.second( x , jacobi , t )

它还要求时间（t），因此它可能不会生成解析形式的解决方案。我会误解文档吗？还是应该使用其他功能？我宁愿留在Boost中，因为我需要将雅可比行列化为ublas::matrix，并且我想避免转换。

编辑：

更具体地说，我将在rosenbrock4 ODE求解器中使用Jacobian。示例here-第47-52行。我需要自动生成此结构，因为ODE集可能会在以后更改，并且我想避免每次手动重写Jacobian。另外，ODE定义中的某些变量在时间上不是恒定的。

最佳答案

雅可比行列式基于函数的导数。如果函数f仅在运行时已知（并且没有线性等约束），则必须自动执行微分。如果您希望准确地做到这一点（与数字估计相反），则需要使用symbolic computation。查找例如here和here以获得支持此功能的库。

请注意，雅可比行列式通常取决于状态和时间，因此不可能将其表示为常数矩阵（例如在您的示例中），除非您的问题很无聊以至于您可以通过分析来解决。