需要帮助解决python中的二阶非线性ODE [英] Need help solving a second order non-linear ODE in python
问题描述
我真的不知道从哪里开始这个问题,因为我对此没有很多经验,但是需要使用计算机来解决项目的这一部分.
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.
我有一个二阶ODE,即:
I have a 2nd order ODE which is:
m = 1220
k = 35600
g = 17.5
a = 450000
和b在1000到10000之间,增量为500.
and b is between 1000 and 10000 with increments of 500.
x(0)= 0
x'(0)= 5
m*x''(t) + b*x'(t) + k*x(t)+a*(x(t))^3 = -m*g
我需要找到最小的b,以便解永远不会是正数.我知道图形应该是什么样子,但是我只是不知道如何使用odeint来获得微分方程的解. 这是我到目前为止的代码:
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:
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()
实际上什么都不产生.
有什么想法吗?谢谢
推荐答案
要使用scipy.integrate.odeint
求解二阶ODE,应将其编写为一阶ODE系统:
To solve a second-order ODE using scipy.integrate.odeint
, you should write it as a system of first-order ODEs:
我先定义z = [x', x]
,然后定义z' = [x'', x']
,这就是您的系统!当然,您必须插入自己的真实关系:
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:
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]
z[0]' = -1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g)
z[1]' = z[0]
或者,只需将其命名为d(z)
:
Or, just call it 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
:
Now you can feed it to the odeint
as such:
_, x = odeint(d, x0, t).T
(_
是我们制作的x'
变量的空白占位符)
(The _
is a blank placeholder for the x'
variable we made)
为了在x
的最大值始终为负的约束条件下最小化b
,可以使用scipy.optimize.minimize
.我将通过实际上最大化x
的最大值来实现它,但要遵守x
保持为负的约束,因为我无法考虑如何在不求函数反转的情况下最小化参数.
In order to minimize b
subject to the constraint that the maximum of x
is 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.
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)
这篇关于需要帮助解决python中的二阶非线性ODE的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!