0

I have a set of second order differential equations:

enter image description here

and I would like to solve for enter image description here using odeint in python. Is this possible?

I have attempted to do it, but my code does not give me the expected result.

import numpy as np, matplotlib.pyplot as plt, matplotlib.font_manager as fm,os
from scipy.integrate import odeint
from mpl_toolkits.mplot3d.axes3d import Axes3D

initial_state = [0.1, 0, 0]
a = 4
time_points = np.linspace(0, 100, 10000)

def my_system(current_state, t):
    theta1, theta2, theta3 = current_state

    d2theta1_dt2 = -a*(2*theta1-theta2-theta3)
    d2theta2_dt2 = -a*(2*theta2-theta1-theta3)
    d2theta3_dt2 = -a*(2*theta3-theta1-theta2)

    return [d2theta1_dt2, d2theta2_dt2, d2theta3_dt2]

xyz = odeint(my_system, initial_state, time_points)

theta1 = xyz[:, 0]
theta2 = xyz[:, 1]
theta3 = xyz[:, 2]

fig = plt.figure(figsize=(12, 9))
ax = fig.gca(projection='3d')
ax.xaxis.set_pane_color((1,1,1,1))
ax.yaxis.set_pane_color((1,1,1,1))
ax.zaxis.set_pane_color((1,1,1,1))
ax.plot(theta1, theta2, theta3, color='g', alpha=0.7, linewidth=0.6)
ax.set_title('plot', fontproperties=title_font)
plt.show()
PiccolMan
  • 4,854
  • 12
  • 35
  • 53
  • 3
    `odeint` expects a 1st order system. Rewrite the system as a system of 6 first order equations. [Youtube example](https://www.youtube.com/watch?v=z3iu1qiJue4) –  Mar 29 '18 at 03:41
  • 2
    The SciPy Cookbook has [an example of rewriting two coupled second order equations as a system of four first order equations](https://scipy-cookbook.readthedocs.io/items/CoupledSpringMassSystem.html). – Warren Weckesser Mar 29 '18 at 04:04

1 Answers1

4

Actually, the fix to insert the first order derivatives is rather quick,

def my_system(current_state, t):
    theta1, theta2, theta3, omega1, omega2, omega3 = current_state

    d2theta1_dt2 = -a*(2*theta1-theta2-theta3)
    d2theta2_dt2 = -a*(2*theta2-theta1-theta3)
    d2theta3_dt2 = -a*(2*theta3-theta1-theta2)

    return [ omega1, omega2, omega3, d2theta1_dt2, d2theta2_dt2, d2theta3_dt2]

where the first three entries of the state are again the angles, and the second three are their velocities.

You will also need to add initial values for the angular velocities to the initial state vector.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51