用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

我试着跟着a similar problem，但就是没法让它工作。有没有可能，我的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

我认为你首先要解决的是一个物理问题，将物理情况转化为一个ODE系统

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解算器中出现奇异雅可比矩阵。最后我得到了x-y平面上的解

与组件