Решение краевой задачи DE в python

Question

Я пытаюсь решить следующий набор DE:

dx' = cos(a)
dy' = sin(a)
dF' = - b * x * cos(a) + sin(a)
da' = (b * x * sin(a) + cos(a)) / F

с условиями:

x(0) = y(0) = x(1) = 0
y(1) = 0.6
F(0) = 0.38
a(0) = -0.5

Я пытался после аналогичная проблема , но Я просто не могу заставить его работать. Возможно ли, что мои F (0) и a (0) полностью отключены, я даже не уверен в них.

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

beta = 5

def fun(x, y):
    x, dx, y, dy,  F, dF, a, da, = y;
    dxds=np.cos(a)
    dyds=np.sin(a)
    dFds=-beta * x * np.cos(a) + np.sin(a)
    dads=(beta * x * np.sin(a) + np.cos(a) ) / F

    return dx, dxds, dy, dyds, dF, dFds, da, dads

def bc(ya, yb):

    return ya[0], yb[0], ya[2], yb[2] + 0.6, ya[4] + 1, yb[4] + 1, ya[6], yb[6]

x = np.linspace(0, 0.5, 10)
y = np.zeros((8, x.size))
y[4] = 0.38
y[6] = 2.5

res = solve_bvp(fun, bc, x, y)
print(res.message)
x_plot = np.linspace(0, 0.5, 200)
plt.plot(x_plot, res.sol(x_plot)[0])

Answer 1

Я думаю, что у вас есть прежде всего физическая проблема, переводящая физическую ситуацию в систему ОДУ.

x(s) и y(s) - это координаты веревки, где s - это длина вдоль веревка. Следовательно, (x'(s),y'(s)) является единичным вектором, который уникально характеризуется своим углом a(s), что дает

x'(s) = cos(a(s))
y'(s) = sin(a(s))

. Чтобы получить форму, теперь нужно рассмотреть механику. Предполагается, что веревка вращается без спирали вокруг оси вращения, оставаясь в одной плоскости. Кроме того, из равновесия сил вы также получаете, что два других уравнения действительно являются уравнениями первого порядка, а не уравнениями второго порядка. Таким образом, ваше состояние имеет только 4 компонента, и поэтому системная функция ODE должна быть

def fun(s, u):
    x, y, F, a = u;
    dxds=np.cos(a)
    dyds=np.sin(a)
    dFds=-beta * x * np.cos(a) + np.sin(a)
    dads=(beta * x * np.sin(a) + np.cos(a) ) / F

    return dxds, dyds, dFds, dads

Теперь доступно только 4 слота граничных условий, которые являются координатами начала и конца каната.

def bc(ua, ub):
    return ua[0], ub[0], ua[1], ub[1] - 0.6

Кроме того, длина интервала для s также является длиной веревки, поэтому значение 0,5 невозможно для заданных координат на полюсе, попробуйте 1,0. Существует некоторый эксперимент, необходимый для получения первоначального предположения, которое не приводит к единственному якобиану в решателе BVP. В конце я получаю решение в плоскости xy

с компонентами