在python中解决边值问题DE [英] Solving a boundary value problem DE in python

问题描述

我正在尝试解决以下DE:

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

符合条件:

 <代码> x(0)= y(0)= x(1)= 0y(1)= 0.6F(0)= 0.38a(0)= -0.5 

我尝试遵循

带有组件

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

解决方案

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

这篇关于在python中解决边值问题DE的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！