Necesita ayuda para resolver una EDO no lineal de segundo orden en Python

Realmente no sé por dónde empezar con este problema, ya que no he tenido mucha experiencia con esto, pero es necesario resolver esta parte del proyecto con una computadora.

Tengo una EDO de 2º orden que es:

m = 1220 k = 35600 g = 17.5 a = 450000 

y b es entre 1000 y 10000 con incrementos de 500.

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

Necesito encontrar la b más pequeña para que la solución nunca sea positiva. Sé cómo debería verse la gráfica, pero no sé cómo usar odeint para obtener una solución a la ecuación diferencial. Este es el código que tengo hasta ahora:

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

Lo que realmente no produce mucho de nada.

¿Alguna idea? Gracias

Para resolver una EDO de segundo orden utilizando scipy.integrate.odeint , debe escribirla como un sistema de EDO de primer orden:

Definiré z = [x', x] , luego z' = [x'', x'] , ¡y ese es tu sistema! Por supuesto, tienes que conectar tus relaciones reales:

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

se convierte en:

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

O, simplemente llámalo 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]' )) 

Ahora puedes alimentarlo a la odeint como tal:

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

(El _ es un marcador de posición en blanco para la variable x' que hicimos)

Para minimizar b sujeto a la restricción de que el máximo de x es siempre negativo, puede usar scipy.optimize.minimize . Lo implementaré al maximizar el máximo de x , sujeto a la restricción de que sigue siendo negativo, porque no puedo pensar en cómo minimizar un parámetro sin poder invertir la función.

 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) 

No creo que pueda resolver su problema como se indica: sus condiciones iniciales, con x = 0 y x' > 0 implican que la solución será positiva para algunos valores muy cercanos al punto de partida. Entonces no hay b para la cual la solución nunca sea positiva …

Dejando eso de lado, para resolver una ecuación diferencial de segundo orden, primero debe volver a escribirla como un sistema de dos ecuaciones diferenciales de primer orden. Definiendo y = x' podemos reescribir su única ecuación como:

 x' = y y' = -b/m*y - k/m*x - a/m*x**3 - g x[0] = 0, y[0] = 5 

Entonces tu función debería verse algo como esto:

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

Y puedes resolver y trazar tus ecuaciones haciendo:

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

introduzca la descripción de la imagen aquí