我正在尝试修改此处的示例：
# include <cppad/cppad.hpp>
namespace { // ---------------------------------------------------------
// define the template function JacobianCases<Vector> in empty namespace
template <typename Vector>
bool JacobianCases()
{ bool ok = true;
using CppAD::AD;
using CppAD::NearEqual;
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
using CppAD::exp;
using CppAD::sin;
using CppAD::cos;
// domain space vector
size_t n = 2;
CPPAD_TESTVECTOR(AD<double>) X(n);
X[0] = 1.;
X[1] = 2.;
// declare independent variables and starting recording
CppAD::Independent(X);
// a calculation between the domain and range values
AD<double> Square = X[0] * X[0];
// range space vector
size_t m = 3;
CPPAD_TESTVECTOR(AD<double>) Y(m);
Y[0] = Square * exp( X[1] );
Y[1] = Square * sin( X[1] );
Y[2] = Square * cos( X[1] );
// create f: X -> Y and stop tape recording
CppAD::ADFun<double> f(X, Y);
// new value for the independent variable vector
Vector x(n);
x[0] = 2.;
x[1] = 1.;
// compute the derivative at this x
Vector jac( m * n );
jac = f.Jacobian(x);
/*
F'(x) = [ 2 * x[0] * exp(x[1]) , x[0] * x[0] * exp(x[1]) ]
[ 2 * x[0] * sin(x[1]) , x[0] * x[0] * cos(x[1]) ]
[ 2 * x[0] * cos(x[1]) , -x[0] * x[0] * sin(x[i]) ]
*/
ok &= NearEqual( 2.*x[0]*exp(x[1]), jac[0*n+0], eps99, eps99);
ok &= NearEqual( 2.*x[0]*sin(x[1]), jac[1*n+0], eps99, eps99);
ok &= NearEqual( 2.*x[0]*cos(x[1]), jac[2*n+0], eps99, eps99);
ok &= NearEqual( x[0] * x[0] *exp(x[1]), jac[0*n+1], eps99, eps99);
ok &= NearEqual( x[0] * x[0] *cos(x[1]), jac[1*n+1], eps99, eps99);
ok &= NearEqual(-x[0] * x[0] *sin(x[1]), jac[2*n+1], eps99, eps99);
return ok;
}
} // End empty namespace
# include <vector>
# include <valarray>
bool Jacobian(void)
{ bool ok = true;
// Run with Vector equal to three different cases
// all of which are Simple Vectors with elements of type double.
ok &= JacobianCases< CppAD::vector <double> >();
ok &= JacobianCases< std::vector <double> >();
ok &= JacobianCases< std::valarray <double> >();
return ok;
}
我正在尝试通过以下方式对其进行修改：
jac令 G 为在本例中计算的雅可比行列式，如下所示：
jac = f.Jacobian(x);
并且，如示例中所示，设X为自变量。我想构造一个新函数 ，H它是 的函数jac，即H(jacobian(X))= 某物，使得 H 是可自微的。一个例子可以是H(X) = jacobian( jacobian(X)[0])，即jacobian(X)wrt的第一个元素的雅可比X（某种二阶导数）。
问题是jac这里写的是 type Vector，它是 raw 上的参数化类型double，而不是AD<double>. 据我所知，这意味着输出不可自微。
我正在寻找一些关于是否可以在更大的运算中使用雅可比行列式的建议，并采用该较大运算的雅可比行列式（与任何算术运算符不同），或者这是否不可能。
编辑：这已经提出了一次赏金，但我再次提出来看看是否有更好的解决方案，因为我认为这很重要。更清楚一点，“正确”答案需要的元素是：
a) 一种计算任意阶导数的方法。
b) 一种不必先验地指定导数顺序的智能方法。如果在编译时必须知道最大阶导数，则无法通过算法确定导数的阶数。此外，在给出的当前答案中指定一个非常大的订单将导致内存分配问题，我想还会导致性能问题。
c) 将衍生订单的模板从最终用户中抽象出来。这很重要，因为很难跟踪所需的衍生品的顺序。如果 b) 得到解决，这可能是“免费”提供的。
如果有人能破解这个，那将是一个了不起的贡献和一个非常有用的操作。