# Need help solving a second-order nonlinear 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.

I have a 2nd order ODE which is:

``````m = 1220

k = 35600

g = 17.5

a = 450000
```
```

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
```
```

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

Which doesn't really produce much of anything.

Any thoughts? Thanks

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

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
```
```

becomes:

`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)`:

``````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 `odeint` as such:

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

(The `_` is a blank placeholder for the `x'` variable we made)

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