python - 求解 Frenet 框架的非线性 ODE 系统

Question

我查过python non linear ODE with 2 variables ，这不是我的情况。也许我的情况不被称为非线性ODE，请纠正我。

问题实际上是Frenet Frame，其中有3个向量T(s)、N(s)和B(s)；参数s>=0。有 2 个标量，其数学公式表达式 t(s) 和 k(s) 已知。我有初始值 T(0)、N(0) 和 B(0)。

diff(T(s), s) =             k(s)*N(s)
diff(N(s), s) = -k(s)*T(s)            + t(s)*B(s)
diff(B(s), s) =            -t(s)*N(s)

那么如何以数字或符号方式得到 T(s)、N(s) 和 B(s)？

我已经检查了 scipy.integrate.ode 但我根本不知道如何将 k(s)*N(s) 传递到其第一个参数

def model (z, tspan):
   T = z[0]
   N = z[1]
   B = z[2]

   dTds =            k(s) * N           # how to express function k(s)?      
   dNds = -k(s) * T          + t(s) * B
   dBds =           -t(s)* N

   return [dTds, dNds, dBds]


z = scipy.integrate.ode(model, [T0, N0, B0]

Answer 1

这是使用 solve_ivp 的代码来自 Scipy 的接口(interface)(而不是 odeint)来获取数值解:

import numpy as np
from scipy.integrate import solve_ivp

from scipy.integrate import cumtrapz
import matplotlib.pylab as plt

# Define the parameters as regular Python function:
def k(s):
    return 1

def t(s):
    return 0

# The equations: dz/dt = model(s, z):
def model(s, z):
    T = z[:3]   # z is a (9, ) shaped array, the concatenation of T, N and B 
    N = z[3:6]
    B = z[6:]

    dTds =            k(s) * N  
    dNds = -k(s) * T          + t(s) * B
    dBds =           -t(s)* N

    return np.hstack([dTds, dNds, dBds])


T0, N0, B0 = [1, 0, 0], [0, 1, 0], [0, 0, 1] 

z0 = np.hstack([T0, N0, B0])

s_span = (0, 6) # start and final "time"
t_eval = np.linspace(*s_span, 100)  # define the number of point wanted in-between,
                                    # It is not necessary as the solver automatically
                                    # define the number of points.
                                    # It is used here to obtain a relatively correct 
                                    # integration of the coordinates, see the graph

# Solve:
sol = solve_ivp(model, s_span, z0, t_eval=t_eval, method='RK45')
print(sol.message)
# >> The solver successfully reached the end of the integration interval.

# Unpack the solution:
T, N, B = np.split(sol.y, 3)  # another way to unpack the z array
s = sol.t

# Bonus: integration of the normal vector in order to get the coordinates
#        to plot the curve  (there is certainly better way to do this)
coords = cumtrapz(T, x=s)

plt.plot(coords[0, :], coords[1, :]);
plt.axis('equal'); plt.xlabel('x'); plt.xlabel('y');

T、N 和 B 是向量。因此，需要求解 9 个方程:z 是一个 (9,) 数组。

对于恒定曲率且无扭转，结果是一个圆: