
<Callout type="success" heading="💬 Mid-course check-in - how's it going??">

I hope you’re enjoying the course so far.

<Image
	src="/images/mid-course-check-in.jpg"
	alt="How's the course going? | EngineeringSkills.com"
	width="400"
	height="438"
/>

I’d appreciate it if you would take 5 minutes to let me know how it’s going by leaving a quick mid-course review. You can update your review again at the end of the course.

[Complete your review here](/leave-a-review) by selecting your course from the dropdown and telling me how it's going.

Thanks, and enjoy the rest of the course.

Seán

</Callout>

The SciPy library contains many helpful functions and tools. Among them is a numerical ordinary [differential equation solver](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html). We won't delve deep into its inner workings here, there is plenty of information available online if you want to explore further. Our only goal is to explain how to use the tool and then to verify the response we've calculated using our implementation of the central difference method. Going forward, if you have differential equations to solve, you can either implement your own algorithm as we've done, or more practically, just use the mature tools already available.

To use the SciPy ODE solver we need to provide it with a system of first order differential equations. You'll recall that we have a system of second order differential equations, so we need to rearrange our system.

This is easier than you might imagine and is a pretty formulaic and routine process once you get familiar with it. We need to introduce new variables to represent the velocity of each mass in our system,

$$
y_1 = \dot{u}_1
\tag{1}
$$

$$
y_2 = \dot{u}_2
\tag{2}
$$

$$
y_3 = \dot{u}_3
\tag{3}
$$

This allows us to write the acceleration expression as the first derivative of $y$ rather than the second derivative of $u$. Remember this is simply a requirement of the ODE solver tool we're using. Matlab's (a mathematical modelling programme) popular ODE45 solver imposes the same requirement. After we rearrange the equations of motion for each mass to make the acceleration the subject of the equation, our system of 3 coupled second order differential equations can be represented by a system of 6 coupled first order equations. Note that I have unpacked the equations from their matrix representation.

$$
\dot{u}_1 = y_1
\tag{4}
$$

$$
\tag{5}
\dot{y}_1 = \ddot{u}_1 = [-\dot{u}_1(c_1+c_2)  + \dot{u}_2 c_2 -u_1(k_1 + k_2) + u_2 k_2]/m_1
$$

$$
\dot{u}_2 = y_2
\tag{6}
$$

$$
\dot{y}_2 = \ddot{u}_2 = [\dot{u}_1c_2-\dot{u}_2(c_2+c_3)+\dot{u}_3c_3+u_1k_2-u_2(k_2+k_3)+u_3k_3]/m_2
\tag{7}
$$

$$
\dot{u}_3 = y_3
\tag{8}
$$

$$
\dot{y}_3 = \ddot{u}_3 = [p_3 + \dot{u}_2c_3-\dot{u}_3c_3+u_2k_3-u_3k_3]/m_3
\tag{9}
$$

This format can now be solved by direct integration using SciPy. First we need to define a function that will return this system of equations. Then this function can be passed into the odeint solver along with, system initial conditions, a time vector and any additional dynamic system parameters. ODEINT will then work its magic and return a matrix where the columns contain all position or velocity values for each mass.

Remember to import SciPy into your notebook. You'll almost certainly have this installed (comes with Anaconda), if not it can be installed at the terminal just like any other python library.

```python
# Use ODEINT to solve the differential equations
from scipy.integrate import odeint
```

Now we can define the function that returns the system of equations to be solved.

```python
def systemEquations(inits,t,params,):
"""
Defines the differential equations for the coupled system

var: vector of the state variables: var = [u1,y1,u2,y2,u3,y3]
    u1 = mass 1 position
    y1 = mass 1 velociity
    u2 = mass 2 position
    y2 = mass 2 velocity
    u3 = mass 3 position
    y3 = mass 3 velocity
t: time
params:  vector of parameters: params = [m1,m2,m3,k1,k2,k3,c1,c2,c3]
"""

# Unpack position and velocity variables and the system parameters
u1, y1, u2, y2, u3, y3 = inits # Position and velocity variables
m1, m2, m3, k1, k2, k3, c1, c2, c3, P, w = params # System parameters

#Define the force vector
F = F3*np.sin(w*t)

# Define a system of first order differential equations
# (introduce dummy variables y1 and y2 for this)
f = [y1,
     ( -y1*(c1+c2) + y2*c2 -(k1+k2)*u1 + k2*u2)/m1,
     y2,
     (y1*c2 - y2*(c2+c3) + y3*c3 + u1*k2 - u2*(k2+k3) + u3*k3)/m2,
    y3,
    (F + y2*c3 - y3*c3 + u2*k3 - u3*k3)/m3]

return f
```

Next we define our initial conditions and pass those along with the other system parameters and our previously defined function into ODEINT.

```python
# Initial conditions
# x1, x2 and x3 are the initial displacements; y1, y2 and y3 are the initial velocities
x1 = 0.0 #initial displacement of mass 1
y1 = 0.0 #initial velocity of mass 1
x2 = 0.0 #initial displacement of mass 2
y2 = 0.0 #initial velocity of mass 2
x3 = 0.0 #initial displacement of mass 3
y3 = 0.0 #initial velocity of mass 3

# Pack up the parameters and initial conditions:
params = [m1, m2, m3, k1, k2, k3, c1, c2, c3, F3, w]
inits = [x1, y1, x2, y2, x3, y3]

# Call the ODE solver.
sol = odeint(systemEquations, inits, t, args=(params,))

#Assign displacements to new variables
U1_direct = sol[:,0] #Displacement of mass 1
U2_direct = sol[:,2] #Displacement of mass 2
U3_direct = sol[:,4] #Displacement of mass 3

#Assign velocities to new variables
V1_direct = sol[:,1] #Velocity of mass 1
V2_direct = sol[:,3] #Velocity of mass 2
V3_direct = sol[:,5] #Velocity of mass 3
```

Next we can plot the output.

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1.5])
axes.plot(t,U1_direct,label='Mass/floor 1')
axes.plot(t,U2_direct,label='Mass/floor 2')
axes.plot(t,U3_direct,label='Mass/floor 3')
axes.set_xlim([0,tmax])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Structural Displacement (m)')
axes.set_title('Forced vibration response (ODEINT)')
axes.grid()
axes.legend(loc='lower right')
plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/using-scipys-odeint-solver/img1.jpg"
	alt="Forced vibration response of each mass/floor level. | EngineeringSkills.com"
	width="1000"
	height="521"
	caption="Fig 1. Forced vibration response of each mass/floor level."
/>

### Compare solutions

Now we can compare the solutions for our system response obtained using our implementation of the central difference method with the ODEINT solution.

```python
fig, axes = plt.subplots(figsize=(12,14),nrows=3,ncols=1)
axes[0].plot(t,U1,'-',label='Central diff')
axes[0].plot(t,U1_direct,'--',label='SciPy ODEINT')
axes[0].set_xlim([0,tmax])
axes[0].set_ylabel('Structural Displacement (m)')
axes[0].set_title('Mass 1')
axes[0].grid()
axes[0].legend(loc='lower right')

axes[1].plot(t,U2,'-',label='Central diff')
axes[1].plot(t,U2_direct,'--',label='SciPy ODEINT')
axes[1].set_xlim([0,tmax])
axes[1].set_ylabel('Structural Displacement (m)')
axes[1].set_title('Mass 2')
axes[1].grid()
axes[1].legend(loc='lower right')

axes[2].plot(t,U3,'-',label='Central diff')
axes[2].plot(t,U3_direct,'--',label='SciPy ODEINT')
axes[2].set_xlim([0,tmax])
axes[2].set_xlabel('Time (sec)')
axes[2].set_ylabel('Structural Displacement (m)')
axes[2].set_title('Mass 3')
axes[2].grid()
axes[2].legend(loc='lower right')

plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/using-scipys-odeint-solver/img2.jpg"
	alt="Comparison of solutions obtained using our central difference algorithm and the ODEINT solver toolbox for mass 1, 2 and 3 (top to bottom). | EngineeringSkills.com"
	width="1000"
	height="1083"
	caption="Fig 2. Comparison of solutions obtained using our central difference algorithm and the ODEINT solver toolbox for mass 1, 2 and 3 (top to bottom)."
/>

As we can see there is pretty much perfect agreement between both numerical solutions. This obviously gives us good confidence in the validity of the code we've written.

### Free vibration

As a final demonstration of ODEINT, we can simulate the free vibration response of the system simply by setting the force amplitude to zero and specifying some non-zero initial conditions.

```python
#Extract max steady-state positions
u1_init = max(U1[2000::])
u2_init = max(U2[2000::])
u3_init = max(U3[2000::])

x1 = u1_init #initial displacement of mass 1
y1 = 0.0 #initial velocity of mass 1
x2 = u2_init #initial displacement of mass 2
y2 = 0.0 #initial velocity of mass 2
x3 = u3_init #initial displacement of mass 2
y3 = 0.0 #initial velocity of mass 2

F3 = 0 #Set forcing term to zero

# Pack up the parameters and initial conditions:
params = [m1, m2, m3, k1, k2, k3, c1, c2, c3, F3, w]
inits = [x1, y1, x2, y2, x3, y3]

# Call the ODE solver.
sol = odeint(systemEquations, inits, t, args=(params,))

U1_directFree = sol[:,0]
U2_directFree = sol[:,2]
U3_directFree = sol[:,4]
```

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1.5])
axes.plot(t,U1_directFree,label='Mass/floor 1')
axes.plot(t,U2_directFree,label='Mass/floor 2')
axes.plot(t,U3_directFree,label='Mass/floor 3')
axes.set_xlim([0,tmax])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Structural Displacement (m)')
axes.set_title('Free vibration response')
axes.grid()
axes.legend(loc='lower right')
plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/using-scipys-odeint-solver/img3.jpg"
	alt="Free vibration response of each mass/floor level. | EngineeringSkills.com"
	width="1000"
	height="517"
	caption="Fig 3. Free vibration response of each mass/floor level."
/>


Please log in or enroll to continue

Please log in or enroll to continue