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])

What exactly are the boundary conditions and what are the guesses of the solution? Your system has 4 equations and one parameter, in the implementation the DE are second order, so you need 8 or 9 boundary conditions, you gave only 6. Also, the BC in the text are different from the BC in the code, please make that consistent. — Lutz Lehmann, Mar 23 '20 at 08:13
Your code works perfectly well, where exactly do expectations and results differ? — Lutz Lehmann, Mar 23 '20 at 09:20
@MikeO'Connor : It is unpacking the state vector into the individual position and velocity variables, tuple assignment. — Lutz Lehmann, Mar 23 '20 at 09:22
@LutzLehmann the problem are the conditions F and a. As I understand the problem F(0) and a(0) are the starting parameters and don't have boundary values. I should get something like this: [link](https://i.ibb.co/x121SdH/Capture.png). It is supposed to be a rope, tied at two points to a rotating post. F(0) is the force at the top knot and a(0) the angle of the rope, again at the top knot. — JerihoC, Mar 23 '20 at 16:20

Lutz Lehmann · Accepted Answer · 2021-03-12T15:14:20.027

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

edited Mar 12 '21 at 15:14

answered Mar 23 '20 at 18:26

Lutz Lehmann

25,219
2
22
51

Thank you very much. I however still can't get the correct result. Initial guesses about F(0) and a(0) you mean? Desperate I tried looping through several hundred F and a around 0.4 and -0.5 with no luck. In every case I got a singular Jacobian. Does the number of breakpoints play a vital role? You used 22, correct?. I am missing something else. I will work on the problem, but any additional help would be great. – JerihoC Mar 24 '20 at 18:28
This might be a case where single shooting works better than the integrated multiple shooting of the standard bvp solvers. I'll look in the possibility of making the length variable. // The initialization is with an ODE solution with 10 points, where a suitable one was determined by trying several graphs. – Lutz Lehmann Mar 24 '20 at 18:44
Yes, single shooting using `scipy.optimize.fsolve` as non-linear solver works a little bit more robustly. However, non-physical initial points still lead to strange failed solutions. – Lutz Lehmann Mar 24 '20 at 19:24
Hello, @LutzLehmann I have learned a lot from your answer. Allow me to make question here. How do you proceed if having `fun()` with extra arguments? In my case I have `fun(u,t,A,B,D,y)` being A,B,D scalars and y an array. How can I use your code for my case? Thanks – Andres Valdez Mar 12 '21 at 15:05
@AndresValdez : As little as there is of general code. Yes, using parameters is possible. The newer scipy procedures are trending to add back the `args=` optional argument, `solve_ivp` has it, `solve_bvp` not yet. So you have to use a `lambda` construct as wrapper, `lambda u,t: fun(u,t,A,B,D,y)`. This changes if the parameters are actually also to be fitted, then use the `p=` optional parameter, but take care to define the functions accordingly. – Lutz Lehmann Mar 12 '21 at 15:13
@LutzLehmann thanks will give a try on that – Andres Valdez Mar 12 '21 at 15:16
@LutzLehmann I tried to combine your advice with this answer [https://stackoverflow.com/a/51734783/9185675]. In my case I have `solve_bvp(lambda y , t: adj_sys(y,t,voc_sol,data,A,B,D),lambda y0, y1, y2, y3: bc(y0, y1, y2, y3),t, voc_sol)`, Notice that `voc_sol` has 4 arrays of Nsize, but `y` has one array of Nsize... How can I do it properly? – Andres Valdez Mar 12 '21 at 15:47
Hopefully you used the correct argument order, `t,y`. I was apparently stuck on the order for `odeint`. `bc` should only have 2 argument arrays, the method does not support multi-point conditions (like matlab's bvp4c does). About the array question see https://stackoverflow.com/questions/61046965/scipy-odeint-different-values-of-arguments-in-each-timestep, https://stackoverflow.com/questions/43830215/defining-a-fuction-to-be-used-in-solving-equation-from-a-data-file – Lutz Lehmann Mar 12 '21 at 16:13
@LutzLehmann the fact is that my `adj_sys()` must return 4 first order derivatives `return [x_adj_der,u_adj_der,y_adj_der,v_adj_der]` – Andres Valdez Mar 12 '21 at 16:22
That does not make sense, the ODE function should have as variable input `t` and the vector `[x,u,y,v]` of length 4, this is the basis of being an ODE function. One can not predict where this ODE function will be evaluated, so any data series has to be enhanced to some kind of interpolation function. Essentially, you can not use `data[n]`, as `n` is not given in the function arguments. You have to translate `t` into `n` and then find an associated intermediate value, which is the nature of interpolation. – Lutz Lehmann Mar 12 '21 at 16:33
Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/229833/discussion-between-andres-valdez-and-lutz-lehmann). – Andres Valdez Mar 12 '21 at 17:08
@LutzLehmann, after our talk I posted a new question. I will appreciate if you can guide me to "join" in one big `solve_bvp()` statement... as you suggested me. Thanks. https://stackoverflow.com/q/66606600/9185675 – Andres Valdez Mar 14 '21 at 23:31

Solving a boundary value problem DE in python

1 Answers1

Linked