numdifftools.Jacobian 的时间分量不应为零

Question

我正在使用 solve_ivp 对微分方程进行数值求解，以便从初始条件对位置和速度进行积分。这在 Python 中很简单，也很容易做到。但我也在计算最终位置和速度相对于初始 positions/velocities 和使用 numdifftools.Jacobian 的时间的部分。这个我比较难掌握。让我解释一下。

我有一组微分方程（参见 rhs_CR3BP 函数），我在给定初始状态的情况下使用 Python 求解器 solve_ivp 在时间跨度 [0,T] 上对其进行数值积分矢量 X_0 在 t=0（位置和速度）。在集成结束时，我将所有从 t=0 到 t=T 的状态向量存储在 solve_ivp.y 中，最终状态向量 X_f 在 t=T 存储在 solve_ivp.y[:, -1] 中.

我正在尝试计算最终状态向量 X_f 相对于初始状态向量 X_0 和周期 T 的数值偏导数，它们实际上是雅可比矩阵。

考虑以下向量：

X = [x_0, y_0, z_0, vx_0, vy_0, vz_0, T] # this is X_0 and T

X_f 相对于 X 的向量偏导数由以下雅可比矩阵给出：

为了计算这个矩阵，我创建了一个函数（参见 propagation 函数），它将一个由 X_0 和 T 组成的向量作为变量，称为 X，和 returns X_f。然后，此函数与 numdifftools.Jacobian 结合计算雅可比矩阵。当我运行我的代码（见下文）时，雅可比矩阵被正确计算为 X_f 相对于 X_0 的偏导数（前 6 列），但最后一列带有偏导数X_f 相对于 T 的导数等于零。我知道这不可能，因为这些导数实际上是 t=T 时的最终速度和加速度。

为什么在X_f的演化中考虑了初始条件X_0而没有考虑时间T？那没有任何意义。在 solve_ivp 中，时间跨度 (tspan=(0,T)) 似乎被完全忽略了，但初始条件 (y0=X_0) 却没有。这是为什么？我没有正确使用 solve_ivp 和 numdifftools.Jacobian 吗？

如果我的问题需要更多细节，请告诉我，我很乐意在一天结束前提供更多解释。

import numpy as np
from scipy.integrate import solve_ivp
import numdifftools as nd

def rhs_CR3BP(t, X0, mu):
    """Integrates the CR3BP equations of motion"""
    
    x, y, z, v_x, v_y, v_z = X0
    
    r = np.sqrt((x-1+mu)**2+y**2+z**2)
    d = np.sqrt((x+mu)**2+y**2+z**2)
    
    a_x = 2*v_y + x - (1-mu)*(x+mu)/d**3 - mu/r**3*(x-1+mu)
    a_y = -2*v_x + y - (1-mu)*y/d**3 - mu*y/r**3
    a_z = -(1-mu)*z/d**3 - mu*z/r**3
    
    return np.array([v_x, v_y, v_z, a_x, a_y, a_z])

def propagation(X, mu):
    
    X_0 = X[0:6]
    T = X[6].real
    
    sol = solve_ivp(fun=rhs_CR3BP, t_span=(0, T), y0=X_0, method='RK45', rtol=1e-10, atol=1e-10, args=(mu,), dense_output=True)
    X_f = sol.y[:, -1]
    
    return X_f

T = 0.732 # period
x, y, z, vx, vy, vz = 1.085, 0, 0, 0, -0.464, 0 # position and velocity components
X0 = [x, y, z, vx, vy, vz] # initial state vector
X = np.array([x, y, z, vx, vy, vz, T]) # X_0 and T

mu = 0.012 # arbitrary parameter
f = lambda X: propagation(X, mu)
f = nd.Jacobian(f, method='complex')
DF = np.real(f(X))

Answer 1

用 nd.Jacobian(f, method='central') 替换 nd.Jacobian(f, method='complex') 产生以下矩阵：

array([[ -1.809,   0.609,   0.   ,  -0.16 ,   1.148,   0.   ,   0.027],
       [ -5.119,   3.386,   0.   ,  -1.052,   1.34 ,   0.   ,   0.468],
       [  0.   ,   0.   ,  -0.849,   0.   ,   0.   ,   0.153,   0.   ],
       [-22.765,  10.802,   0.   ,  -3.568,   8.253,   0.   ,   1.907],
       [  0.836,  -0.096,   0.   ,   0.103,   0.621,   0.   ,  -0.156],
       [  0.   ,   0.   ,  -1.532,   0.   ,   0.   ,  -0.902,   0.   ]])

最后一栏中 df/dT 的这些条目看起来更合理。

您的所有变量似乎都是实值，但您在调用 Jacobian 时请求了 method=complex。这是自找麻烦，原因我将在下面展开。

实数可微与复数可微

如 numdifftools 文档中所述：

Complex methods are usually the most accurate provided the function to differentiate is analytic

解析函数是可以用幂级数表示的函数，解析函数与complex-differentiable函数：这两个类函数是相同的（跳过一些技术细节）。相比之下，有许多函数沿实数线可微但不是解析的。

重要的是，仅仅因为一个函数对于实数是可微的，并不意味着它也是复可微的（即解析的）。您的函数保证相对于 T 是实可微的，但不能保证它对于复数 T.

是复可微的

复可微函数比实可微函数多 select 组，并且它们具有许多并非所有实可微函数都具有的属性（除了解析）。例如，如果一个复函数可以微分一次，那么它就可以微分无数次。 f(x) = |x|**3.

是非复可微的实可微函数示例

numdifftools.Jacobian 的时间分量不应为零

Time components of numdifftools.Jacobian should not be equal to zero

python

python-3.x

math

ode

solver

实数可微与复数可微