需要帮助解决 Python 中的二阶非线性 ODE

Need help solving a second order non-linear ODE in python

提问by Max Rackoff

I don't really know where to start with this problem, as I haven't had much experience with this but it is required to solve this part of the project using a computer.

我真的不知道从哪里开始解决这个问题，因为我对此没有太多经验，但是需要使用计算机解决项目的这一部分。

I have a 2nd order ODE which is:

我有一个二阶 ODE，它是：

m = 1220

k = 35600

g = 17.5

a = 450000

and b is between 1000 and 10000 with increments of 500.

b 介于 1000 和 10000 之间，增量为 500。

x(0)= 0 

x'(0)= 5


m*x''(t) + b*x'(t) + k*x(t)+a*(x(t))^3 = -m*g

I need to find the smallest b such that the solution is never positive. I know what the graph should look like, but I just don't know how to use odeint to get a solution to the differential equation. This is the code I have so far:

我需要找到最小的 b 使得解永远不会是正的。我知道图形应该是什么样子，但我只是不知道如何使用 odeint 来求解微分方程。这是我到目前为止的代码：

from    numpy    import    *    
from    matplotlib.pylab    import    *    
from    scipy.integrate    import    odeint

m = 1220.0
k = 35600.0
g  = 17.5
a = 450000.0
x0= [0.0,5.0]

b = 1000

tmax = 10
dt = 0.01

def fun(x, t):
    return (b*x[1]-k*x[0]-a*(x[0]**3)-m*g)*(1.0/m)
t_rk = arange(0,tmax,dt)   
sol = odeint(fun, x0, t_rk)
plot(t_rk,sol)
show()

Which doesn't really produce much of anything.

这并没有真正产生任何东西。

Any thoughts? Thanks

有什么想法吗？谢谢

回答by Jaime

I don't think you can solve your problem as stated: your initial conditions, with x = 0and x' > 0imply that the solution will be positive for some values very close to the starting point. So there is no bfor which the solution is never positive...

我不认为您可以按照说明解决您的问题：您的初始条件，x = 0并x' > 0暗示对于非常接近起点的某些值，该解决方案将是正的。所以没有b哪个解决方案永远不会是积极的......

Leaving that aside, to solve a second order differential equation, you first need to rewrite it as a system of two first order differential equations. Defining y = x'we can rewrite your single equation as:

撇开这一点，要求解二阶微分方程，首先需要将其重写为两个一阶微分方程的系统。定义y = x'我们可以将您的单个方程重写为：

x' = y
y' = -b/m*y - k/m*x - a/m*x**3 - g

x[0] = 0, y[0] = 5

So your function should look something like this:

所以你的函数应该是这样的：

def fun(z, t, m, k, g, a, b):
    x, y = z
    return np.array([y, -(b*y + (k + a*x*x)*x) / m - g])

And you can solve and plot your equations doing:

你可以解决和绘制你的方程：

m, k, g, a = 1220, 35600, 17.5, 450000
tmax, dt = 10, 0.01
t = np.linspace(0, tmax, num=np.round(tmax/dt)+1)
for b in xrange(1000, 10500, 500):
    print 'Solving for b = {}'.format(b)
    sol = odeint(fun, [0, 5], t, args=(m, k, g, a, b))[..., 0]
    plt.plot(t, sol, label='b = {}'.format(b))
plt.legend()

回答by askewchan

To solve a second-order ODE using scipy.integrate.odeint, you should write it as a system of first-order ODEs:

要使用 求解二阶 ODE scipy.integrate.odeint，您应该将其编写为一阶 ODE 系统：

I'll define z = [x', x], then z' = [x'', x'], and that's your system! Of course, you have to plug in your real relations:

我会定义z = [x', x]，然后z' = [x'', x']，这就是你的系统！当然，你必须插入你的真实关系：

x'' = -(b*x'(t) + k*x(t) + a*(x(t))^3 + m*g) / m

becomes:

变成：

z[0]' = -1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g)
z[1]' = z[0]

z[0]' = -1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g)
z[1]' = z[0]

Or, just call it d(z):

或者，只需调用它d(z)：

def d(z, t):
    return np.array((
                     -1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g),  # this is z[0]'
                     z[0]                                         # this is z[1]'
                   ))

Now you can feed it to the odeintas such:

现在您可以将其提供给以下内容odeint：

_, x = odeint(d, x0, t).T

(The _is a blank placeholder for the x'variable we made)

（这_是x'我们创建的变量的空白占位符）

In order to minimize bsubject to the constraint that the maximum of xis always negative, you can use scipy.optimize.minimize. I'll implement it by actually maximizing the maximum of x, subject to the constraint that it remains negative, because I can't think of how to minimize a parameter without being able to invert the function.

为了b在最大值x始终为负的约束下最小化，您可以使用scipy.optimize.minimize. 我将通过实际最大化 的最大值来实现它x，受到它仍然为负的约束，因为我想不出如何在不能够反转函数的情况下最小化参数。

from scipy.optimize import minimize
from scipy.integrate import odeint

m = 1220
k = 35600
g = 17.5
a = 450000
z0 = np.array([-.5, 0])

def d(z, t, m, k, g, a, b):
    return np.array([-1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g), z[0]])

def func(b, z0, *args):
    _, x = odeint(d, z0, t, args=args+(b,)).T
    return -x.max()  # minimize negative max

cons = [{'type': 'ineq', 'fun': lambda b: b - 1000, 'jac': lambda b: 1},   # b > 1000
        {'type': 'ineq', 'fun': lambda b: 10000 - b, 'jac': lambda b: -1}, # b < 10000
        {'type': 'ineq', 'fun': lambda b: func(b, z0, m, k, g, a)}] # func(b) > 0 means x < 0

b0 = 10000
b_min = minimize(func, b0, args=(z0, m, k, g, a), constraints=cons)