Solving a boundary value problem DE in python

Question

I am trying to solve the following set of DE's:

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

with the conditions:

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

I tried following a similar problem, but I just can't get it to work. Is it possible, that my F(0) and a(0) are completely off, I am not even sure about them.

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

I think that you have foremost a physics problem, translating the physical situation into an ODE system.

x(s) and y(s) are the coordinates of the rope where s is the length along the rope. Consequently, (x'(s),y'(s)) is a unit vector that is uniquely characterized by its angle a(s), giving

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

To get the shape, one now has to consider the mechanics. The assumption seems to be that the rope rotates without spiraling around the rotation axis, staying in one plane. Additionally, from the equilibrium of forces you also get that the other two equations are indeed first order, not second order equations. So your state only has 4 components and the ODE system function thus has to be

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

Now there are only 4 boundary condition slots available, which are the coordinates of the start and end of the rope.

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

Additionally, the interval length for s is also the rope length, so a value of 0.5 is impossible for the given coordinates on the pole, try 1.0. There is some experimentation needed to get an initial guess that does not lead to a singular Jacobian in the BVP solver. In the end I get the solution in the x-y plane

with the components