Решить нелинейную систему ОДУ для кадра Frenet

Question

Я проверил Python нелинейный ODE с 2 переменными , что не в моем случае. Возможно, мой случай не называется 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]

Answer 1

спасибо за ваш пример.И я снова подумал, обнаружив, что поскольку существует формула для dZ, где Z - матрица (T, N, B), мы можем вычислить Z[i] = Z[i-1] + dZ[i-1]*deltaS в соответствии с концепцией производной.Затем я пишу код и нахожу, что эта идея может решить пример круга.Так

1. подходит ли Z[i] = Z[i-1] + dZ[i-1]*deltaS для других ODE?в некоторых ситуациях произойдет сбой или scipy.integrate.solve_ivp / scipy.integrate.ode обеспечит преимущество перед прямым использованием Z[i] = Z[i-1] + dZ[i-1]*deltaS?
2. в моем коде, я должен нормализовать Z [i], потому что|| Z [I] ||не всегда 1. Почему это происходит?Ошибка числового вычисления с плавающей точкой?

мой ответ на мой вопрос, по крайней мере, он работает для круга

import numpy as np

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

def dZ(s, Z):
    return np.array(
        [k(s) * Z[1], -k(s) * Z[0] + t(s) * Z[2], -t(s)* Z[1]]
    )

T0, N0, B0 = np.array([1, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 1])

deltaS = 0.1    # step to calculate dZ/ds
num = int(2*np.pi*1/deltaS) + 1 # how many points on the curve we have to calculate

T = np.zeros([num, ], dtype=object)
N = np.zeros([num, ], dtype=object)
B = np.zeros([num, ], dtype=object)

T[0] = T0
N[0] = N0
B[0] = B0

for i in range(num-1):
    temp_dZ = dZ(i*deltaS, np.array([T[i], N[i], B[i]]))
    T[i+1] = T[i] + temp_dZ[0]*deltaS
    T[i+1] = T[i+1]/np.linalg.norm(T[i+1])  # have to do this

    N[i+1] = N[i] + temp_dZ[1]*deltaS
    N[i+1] = N[i+1]/np.linalg.norm(N[i+1])

    B[i+1] = B[i] + temp_dZ[2]*deltaS
    B[i+1] = B[i+1]/np.linalg.norm(B[i+1])


coords = cumtrapz(
    [
        [i[0] for i in T], [i[1] for i in T], [i[2] for i in T]
    ]
    , x=np.arange(num)*deltaS
)

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

plt.show()

Answer 2

Я обнаружил, что уравнение, которое я перечислил в первом посте, не работает для моей кривой.Поэтому я прочитал Gray A., Abbena E., Salamon S-Modern Differential Geometry of Curves and Surfaces with Mathematica. 2006 и обнаружил, что для произвольной кривой уравнение Френе должно быть записано в виде

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

, где ||r'|| (или ||r'(s)||) равно diff([x(s), y(s), z(s)], s).norm()

теперьпроблема изменилась, чтобы немного отличаться от той, что была в первом посте, потому что нет функции r '(s) или массива дискретных данных. Так что я думаю, что это подходит для нового ответа, отличного от комментария.

Я встретил 2 вопроса, пытаясь решить новое уравнение:

1. как мы можемзапрограммировать с помощью r '(s), если используется scipy solve_ivp?
2. Я пытаюсь изменить gaussian solution, но результат совершенно неверный.

еще раз спасибо

import numpy as np

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

def dZ(s, Z, r_norm):
    return np.array([
        r_norm * k(s) * Z[1],
        r_norm*(-k(s) * Z[0] + t(s) * Z[2]),
        r_norm*(-t(s)* Z[1])
    ])

T0, N0, B0 = np.array([1, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 1])

deltaS = 0.1    # step to calculate dZ/ds
num = int(2*np.pi*1/deltaS) + 1 # how many points on the curve we have to calculate

T = np.zeros([num, ], dtype=object)
N = np.zeros([num, ], dtype=object)
B = np.zeros([num, ], dtype=object)
R0 = N0

T[0] = T0
N[0] = N0
B[0] = B0

for i in range(num-1):
    r_norm = np.linalg.norm(R0)
    temp_dZ = dZ(i*deltaS, np.array([T[i], N[i], B[i]]), r_norm)
    T[i+1] = T[i] + temp_dZ[0]*deltaS
    T[i+1] = T[i+1]/np.linalg.norm(T[i+1])

    N[i+1] = N[i] + temp_dZ[1]*deltaS
    N[i+1] = N[i+1]/np.linalg.norm(N[i+1])

    B[i+1] = B[i] + temp_dZ[2]*deltaS
    B[i+1] = B[i+1]/np.linalg.norm(B[i+1])

    R0 = R0 + T[i]*deltaS

coords = cumtrapz(
    [
        [i[0] for i in T], [i[1] for i in T], [i[2] for i in T]
    ]
    , x=np.arange(num)*deltaS
)


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

plt.show()

Answer 3

Вот код, использующий интерфейс 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,).

Для постоянной кривизны и без скручивания, в результате получается круг: