从曲率与挠率重建3D曲线：solve_ivp求解ODE系统报错问题

问题背景 #

已有3D曲线坐标，通过splipy计算曲率与挠率：

from splipy import curve_factory
import numpy as np

points = np.load('my_coordinates.npy')
curve = curve_factory.cubic_curve(points)
sample_points = 1000 
t = np.linspace(curve.start()[0], curve.end()[0], sample_points)
curvature = curve.curvature(t)
torsion = curve.torsion(t)

尝试通过求解Frenet-Serret方程组重建曲线，已将离散的曲率、挠率插值为连续函数，并定义了ODE右端函数，但使用scipy.integrate.solve_ivp时报错，疑问是否需要将向量形式的ODE拆分为12个标量方程。

解决方案 #

scipy.integrate.solve_ivp不支持直接传入二维数组作为初始条件或处理向量形式的状态变量，必须将所有状态变量扁平化为一维数组，在右端函数内部再恢复为原结构计算。具体步骤如下：

1. 修改初始条件为一维数组 #

将gamma(0)（起点）、T(0)（切向量）、N(0)（法向量）、B(0)（副法向量）的4个3维向量展平为12个元素的一维数组：

import numpy as np
from scipy import interpolate
from scipy.integrate import solve_ivp

# 加载并插值曲率、挠率
kappa = np.load("curvature.npy")
tao = np.load("torsion.npy")
length = len(kappa)
s_vals = np.linspace(0, 1, length)
kappa_interp = interpolate.interp1d(s_vals, kappa)
tao_interp = interpolate.interp1d(s_vals, tao)

# 初始条件展平：[gamma_x, gamma_y, gamma_z, T_x, T_y, T_z, N_x, N_y, N_z, B_x, B_y, B_z]
initial_cond = np.array([0, 0, 0,  # gamma(0)：原点
                         0, 1, 0,  # T(0)：y方向单位向量
                         1, 0, 0,  # N(0)：x方向单位向量
                         0, 0, 1]) # B(0)：z方向单位向量

2. 修改右端函数，适配一维状态变量 #

在rhs函数中，先将输入的一维数组reshape为(4,3)的矩阵，对应gamma, T, N, B，计算Frenet-Serret方程组后，再将结果展平为一维数组返回：

def rhs(s, x):
    # 将一维状态变量恢复为4个3维向量
    gamma, T, N, B = x.reshape(4, 3)
    k = kappa_interp(s)
    t = tao_interp(s)
    
    # Frenet-Serret方程组
    dgamma_dt = T
    dT_dt = k * N
    dN_dt = -k * T + t * B
    dB_dt = -t * N
    
    # 将结果展平为一维数组返回
    return np.concatenate([dgamma_dt, dT_dt, dN_dt, dB_dt])

3. 调用solve_ivp求解 #

此时传入扁平化的初始条件即可正常求解：

# 求解区间s∈[0,1]
res = solve_ivp(rhs, (0, 1), initial_cond, t_eval=np.linspace(0,1,sample_points))

# 提取重建的曲线点：gamma的x,y,z分量
reconstructed_points = res.y[:3].T

关键说明 #

Frenet-Serret方程组本质是12个耦合的标量ODE，solve_ivp仅能处理一维状态向量，因此必须通过扁平化/重构的方式适配API。求解完成后，从结果中提取前3个元素即可得到重建的3D曲线坐标。

内容的提问来源于stack exchange，提问作者MangoMat