Need help solving a second order non-linear ODE in python
Solution 1
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)
Solution 2
I don't think you can solve your problem as stated: your initial conditions, with x = 0
and x' > 0
imply that the solution will be positive for some values very close to the starting point. So there is no b
for which the solution is never positive...
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:
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()
Admin
Updated on November 06, 2020Comments
-
Admin over 3 years
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
-
Admin over 10 yearsAh, I forgot that the 5 should be negative. Also I don't quite understand how to call x' = y or y' = -b/my - k/mx - a/m*x**3 - g. Wouldn't those just be the starts of strings?
-
askewchan over 10 yearsUsing the prime
'
is just notation to explain the equation, don't put them in your code. The actual code is given byfun
or (d
in my answer), where instead ofx' = y
andy' = blah
, you writefun([x, y]) = [y, blah]
-
Jaime over 10 yearsCool, didn't know that
odeint
supported broadcasting of the arguments. -
askewchan over 10 years@Jaime Actually without fixing the return format of
d(z)
it doesn't, i meant to come back and fix it but forgot. -
usethedeathstar over 10 yearswhat does the "_, x" thing do? especially the _-part
-
askewchan over 10 years@usethedeathstar,
_
is just a variable name, it's equivalent toy, x = odeint(...)
, or evenz = odeint(...)
, sinceodeint
returns an array with first dimension of2
, I can "unpack" it. I chose to use_
because I have no use for they
part of the array, since it representsx'
, which I don't need. I could alternatively have done:x = odeint(...)[1,:]
-
usethedeathstar over 10 yearsso it is basically a dummy that you dont use?
-
askewchan over 10 years@usethedeathstar exactly :) It's just a convention, not official. Unused variable naming in python, Underscore in Python,
_
as a variable name in Python.