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

Question

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

我有一个二阶 ODE，即：

m = 1220

k = 35600

g = 17.5

a = 450000

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

我需要找到最小的 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()

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

有什么想法吗？谢谢

Answer 1

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

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

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

变成：

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

或者，就叫它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]'
                   ))

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

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

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

为了在x 的最大值始终为负的约束下最小化b，您可以使用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)

Answer 2

我认为您无法解决上述问题：您的初始条件，x = 0 和 x' > 0 意味着对于非常接近起点的某些值，解决方案将是积极的。所以没有b 的解决方案永远不会是积极的......

除此之外，要求解二阶微分方程，您首先需要将其重写为两个一阶微分方程的系统。定义y = x'，我们可以将你的单个方程改写为：

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

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

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

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

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

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