¿Alguna forma de resolver un sistema de ecuaciones diferenciales acopladas en python?

He estado trabajando con sympy y scipy, pero no puedo encontrar o averiguar cómo resolver un sistema de ecuaciones diferenciales acopladas (no lineales, de primer orden).

Entonces, ¿hay alguna manera de resolver ecuaciones diferenciales acopladas?

Las ecuaciones son de la forma:

V11'(s) = -12*v12(s)**2 v22'(s) = 12*v12(s)**2 v12'(s) = 6*v11(s)*v12(s) - 6*v12(s)*v22(s) - 36*v12(s) 

con condiciones iniciales para v11 (s), v22 (s), v12 (s).

Para la solución numérica de EDO con scipy, vea la función scipy.integrate.odeint o la clase scipy.integrate.ode .

Algunos ejemplos se dan en el Libro de cocina de SciPy (desplácese hacia abajo hasta la sección “Ecuaciones diferenciales ordinarias”).

Además de los métodos SciPy odeint y ode que ya se mencionaron, ahora tiene solve_ivp que es más nuevo ya menudo más conveniente. Un ejemplo completo, que codifica [v11, v22, v12] como una matriz v :

 from scipy.integrate import solve_ivp def rhs(s, v): return [-12*v[2]**2, 12*v[2]**2, 6*v[0]*v[2] - 6*v[2]*v[1] - 36*v[2]] res = solve_ivp(rhs, (0, 0.1), [2, 3, 4]) 

Esto resuelve el sistema en el intervalo (0, 0.1) con el valor inicial [2, 3, 4] . El resultado tiene variables independientes (s en su notación) como res.t :

 array([ 0. , 0.01410735, 0.03114023, 0.04650042, 0.06204205, 0.07758368, 0.0931253 , 0.1 ]) 

Estos valores fueron elegidos automáticamente. Se puede proporcionar t_eval para evaluar la solución en los puntos deseados: por ejemplo, t_eval=np.linspace(0, 0.1) .

La variable dependiente (la función que estamos buscando) está en res.y :

 array([[ 2. , 0.54560138, 0.2400736 , 0.20555144, 0.2006393 , 0.19995753, 0.1998629 , 0.1998538 ], [ 3. , 4.45439862, 4.7599264 , 4.79444856, 4.7993607 , 4.80004247, 4.8001371 , 4.8001462 ], [ 4. , 1.89500744, 0.65818761, 0.24868116, 0.09268216, 0.0345318 , 0.01286543, 0.00830872]]) 

Con Matplotlib, esta solución se grafica como plt.plot(res.t, res.yT) (la plot sería más suave si proporcionara el t_eval como se mencionó).

diagrama de solución

Finalmente, si el sistema involucrara ecuaciones de orden superior a 1, uno tendría que usar la reducción a un sistema de primer orden .