python - 求解 Frenet 框架的非线性 ODE 系统
问题描述
我已经用 2 个变量检查了 python 非线性 ODE,这不是我的情况。也许我的案子没有被称为nonlinear 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]
解决方案
这是使用solve_ivp
Scipy 的接口(而不是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,) 数组。
对于恒定曲率和无扭转,结果是一个圆:
推荐阅读
- python - 使用 Python C 扩展获取当前工作目录
- mysql - My SQL 5.7.33 错误代码 - 1126。无法打开共享库
- firebase - 无法在firebase实时数据库中的一个节点中获取一个孩子的所有孩子
- sql - 需要支持将 SQL Server 查询转换为 Oracle
- c# - 如何使用 NPOI C# 包创建数据透视表并保存它?
- python - 使用 split 函数应用:'expand' 是 split() 的无效关键字参数
- python - python/dataframe - 合并重复的行
- sql - Laravel 中的 SQL 查询没有返回正确的结果
- javascript - 标准 DOM 元素事件是否在 DOMContentLoaded 之前触发?
- celery - KeyError 收到未注册的“芹菜”类型的任务,而任务已注册