Two-Mass Oscillator: Plotting Amplitudes Over Frequency in Hertz

In summary, the "Two-Mass Oscillator" explores the relationship between the amplitudes of two coupled oscillators as a function of frequency in Hertz. The study emphasizes how varying the frequency affects the amplitude response of each mass, illustrating resonance phenomena and the interaction between the oscillators. Graphical representations demonstrate the amplitude peaks at specific frequencies, highlighting the system's dynamic behavior and the importance of frequency in oscillatory systems.
  • #1
tehsportsmaen
3
0
TL;DR Summary
I try to simulate a two-mass oscillator for a speaker-system and i don't know if my approach makes sense because the physiscal and code-wise details are very complex for me right now.
Idea:
Given a system of two coupled oscillators in which 2 masses are connected to a spring in the middle. Each of the two masses is coupled to another spring on the left and right, which have fixed ends but are not connected to each other. So we have 3 springs, two masses and the springs also have their own damping parameters in addition to the spring constants. Mass m1 is driven by a constant sinusoidal force. So parameters should be: m1, m2, c1, c2, c3, k1, k2, k3, F1. I try to plot the amplitudes over the frequncies in hertz (not omega)
  • Define equations and parameters
  • Calculate Amplitudes of x1,x2, then for-loop with fast fourier tranformation-calculate in a frequency-range (np.logspace) from 10-1000Hz in 81 steps
  • Plot the amplitudes of x1,x2 over the f-range
  • Change parameters to see how the system reacts

Code::
import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from scipy.integrate import odeint
def equations(Y, t, m1, m2, c1, c2, c3, k1, k2, k3, F1, omega):

    x1, v1, x2, v2 = Y

    dx1dt = v1

    dv1dt = (-k1*x1 + k2*(x2 - x1) - c1*v1 + F1*np.sin(omega*t)) / m1

    dx2dt = v2

    dv2dt = (-k3*x2 - k2*(x2 - x1) - ((c2+c3)*v2))/ m2

    return [dx1dt, dv1dt, dx2dt, dv2dt]#SYSTEM

m1, m2 = 0.001, 0.090

c1, c2, c3 = 2.85, 0.5, 1

k1, k2, k3 = 1872.5, 1777.3, 130

F1 = 1.63#TIME AND FREQUENCIES

t = np.linspace(0, 10, 1000)

frequencies = np.logspace(1, 3, 81)def integrate_system(frequency, m1, m2, c1, c2, c3, k1, k2, k3, F1, t):

    omega = 2 * np.pi * frequency

    Y0 = [0, 0, 0, 0]  # initial conditions: [x1, v1, x2, v2]

    solution = odeint(equations, Y0, t, args=(m1, m2, c1, c2, c3, k1, k2, k3, F1, omega))

    return solution[:, 0], solution[:, 2]  #return x1 x2amplitudes_x1 = []

amplitudes_x2 = []for frequency in frequencies:

    x1, x2 = integrate_system(frequency, m1, m2, c1, c2, c3, k1, k2, k3, F1, t)

    #FFT

    fft_x1 = np.fft.fft(x1)

    fft_x2 = np.fft.fft(x2)

    #compute amplitudes

    amplitude_spectrum_x1 = np.abs(fft_x1) / len(fft_x1)

    amplitude_spectrum_x2 = np.abs(fft_x2) / len(fft_x2)

    #index of fft frequencies

    frequency_index = np.argmin(np.abs(frequency - np.fft.fftfreq(len(t), d=t[1]-t[0])))

    #amplitudes of driving frequencies

    amplitudes_x1.append(amplitude_spectrum_x1[frequency_index])

    amplitudes_x2.append(amplitude_spectrum_x2[frequency_index])pd.set_option('display.max_rows', None)

df = pd.DataFrame({'Frequency_Hz': frequencies,'Amplitude_m1': amplitudes_x1,'Amplitude_m2': amplitudes_x2})

print(df)plt.figure(figsize=(12, 6))

plt.loglog(frequencies, amplitudes_x1, label='Amplitude of m1', color='red')

plt.loglog(frequencies, amplitudes_x2, label='Amplitude of m2', color='blue')

plt.xlabel('Frequency (Hz)')

plt.ylabel('Amplitude')

plt.title('LS FREQUENCIES')

plt.legend()

plt.grid(True)

plt.show()
The result looks strange but maybe it makes sense.Would could be potentially difficult with my approach?=> (timearray, integration, FFT, calculation, initial conditions, parameter interpretation)I am happy to get advices, comments and remarks on this.
 
Last edited:
Physics news on Phys.org
  • #2
I have not used odeint in python before, so I'd have to read up on that documentation.

When you put your equations into code, it looks like you account for 2 damping coefficients for m2 but only 1 damping coefficient for m1. Try drawing free body diagrams for each mass.
 
Last edited:
  • #3
Also I think when you have the formula for

dv2dt = (-k3*x2 - k2*(x2 - x1) - ((c2+c3)*v2))/ m2

Is the bolded portion solely dependent on v2, or is it the relative speed? The free body diagrams should help with this.
 
  • #4
It might be easier to first try just doing oscillators with no damping and no driving. That can be solved analytically and you can compare to numerical results. This way you can check the barebones of your code is working properly.
 
  • Like
Likes scottdave
  • #5
scottdave said:
Also I think when you have the formula for

dv2dt = (-k3*x2 - k2*(x2 - x1) - ((c2+c3)*v2))/ m2

Is the bolded portion solely dependent on v2, or is it the relative speed? The free body diagrams should help with this.
Edited Equations:
def equations(Y, t, m1, m2, c1, c2, c3, k1, k2, k3, F1, omega):
    x1, v1, x2, v2 = Y
    dx1dt = v1
    dx1dt = v1
    dv1dt = (-k1*x1 + k2*(x2 - x1) - c1*v1 - c2*(v1 - v2) + F1*np.sin(omega*t)) / m1

    dx2dt = v2
    dv2dt = (-k3*x2 - k2*(x2 - x1) - c3*v2 - c2*(v2 - v1)) / m2
    return [dx1dt, dv1dt, dx2dt, dv2dt]

Is this better?
 
  • #6
1704718581438.png

This is the Output now
 
  • #7
Does this look like your problem?

1719043959669.png


I am still not done with it but I'm trying to calculate the equations of motion for an "n" body system like the one shown. You didn't ask for that generalization but it allows you to do some checks in the process although I'm still struggling with n>2.

For that I'm setting up the equations using Lagrange with the added dissipative term and then finding the inertia, stiffness, and dissipative matrixes.
Now it's where the fun starts and hopefully, I will be able to spend more time on it this weekend because I think it's an interesting problem and good for practice.

For now, for the two-body problem with the initial conditions you described in the OP, this is what I'm getting.

Natural angular velocities (rad/s):
[107.449357049886, 1912.96807550251]

Natural frecuencies (Hz):
[17.1010963065354, 304.458325193215]

Normalized natural modes (m) (Column = vector):
Matrix([[0.438930605481940, 0.999985269730640], [0.898520964458402, -0.00542773633656283]])


If you make all the ##c## smaller and adjust the simulation to calculate from 0→500 Hz which already covers both natural frequencies it might be easier to see the resonance in the frequency response you provided in post #6.
 
  • #8
Juanda said:
You didn't ask for that generalization but it allows you to do some checks in the process although I'm still struggling with n>2.
I'm still struggling with the generalization...

After getting the equations of motion, I wrote the information in matrix form obtaining the inertia matrix ##\bar{\bar{M}}## and the stiffness matrix ##\bar{\bar{K}}##. From that, I am able to calculate the system's natural frequencies and natural modes.
$$\bar{\bar{M}}\ddot{\bar{x}}-\bar{\bar{K}}\bar{x}=0$$
Assuming (knowing) the solution of the differential equation takes the form:
$$\bar{x}(t)=e^{\omega t}\bar{C}$$
Where the vector ##\bar{C}## is just a vector of still unknown constant.
Inputting that in the first matrix differential equation gives that non-trivial solutions must fulfill the following.
$$\det(\bar{\bar{M}}\bar{x}\omega^2-\bar{\bar{K}}\bar{x})=0$$
From that, the natural angular velocities ##\omega## are found. Then, once they are known, they can be plugged into the first system of equations and, by removing one of the equations (the first one or last one because they don't contain all the variables) and assuming the value of one of the constants, the ratio between amplitudes at the natural frequencies is found.

Using this link as a reference I can at least check my results for up to four bodies. My code seems to work fine for up to 3 bodies. Once the fourth one is added, the frequencies are still calculated correctly but the 4th natural mode is not correct.

Here is the output.
The number of bodies is: 4
The mass of each body from 0 to n is: [1, 1, 1, 1] [kg]
The stiffness of each spring from 0 to n+1 is: [1, 1, 1, 1, 1] [N/m]
The inertia matrix is: Matrix([
[1, 0, 0, 0]
[0, 1, 0, 0]
[0, 0, 1, 0]
[0, 0, 0, 1]]) [kg]

The stiffness matrix is: Matrix([
[-2, 1, 0, 0]
[1, -2, 1, 0]
[0, 1, -2, 1]
[0, 0, 1, -2]]) [N/m]

The natural angular velocities are: [0.618, 1.175, 1.618, 1.902] [rad/s]
Normalized natural modes (m) (Column = vector): Matrix([
[0.371, 0.601, 0.601, 0.371]
[0.601, 0.3710, -0.371, 0.601]
[0.601, -0.371, -0.371, -0.371]
[0.371, -0.601, 0.601, -0.601]]) [m]

I checked the natural angular velocities with the link and it's correct.
Due to symmetries in the system, the rows and columns in the last matrix are the same and by comparing it with the link previously mentioned, I'd expect the red values to be -0.601, 0.601, and -0.371. At least that's what I believe it should happen and it is what happens when there are 2 or 3 bodies.
There clearly is something wrong with my code but I believe it might be I'm not indexing the values correctly because the underlying math seems solid.
 
  • #9
Juanda said:
The natural angular velocities are: [0.618, 1.175, 1.618, 1.902] [rad/s]
Normalized natural modes (m) (Column = vector): Matrix([
[0.371, 0.601, 0.601, 0.371]
[0.601, 0.3710, -0.371, 0.601]
[0.601, -0.371, -0.371, -0.371]
[0.371, -0.601, 0.601, -0.601]]) [m]

I checked the natural angular velocities with the link and it's correct.
Due to symmetries in the system, the rows and columns in the last matrix are the same and by comparing it with the link previously mentioned, I'd expect the red values to be -0.601, 0.601, and -0.371. At least that's what I believe it should happen and it is what happens when there are 2 or 3 bodies.
There clearly is something wrong with my code but I believe it might be I'm not indexing the values correctly because the underlying math seems solid.
Indeed, indexation was the culprit.
As I was printing the middle steps to check the result I realized the following.

sol
Out[109]: {NM_10: 1.61803398874994, NM_11: -1.00000000000003, NM_9: -1.61803398874993}

The solution vector is not in order. I'll fix that and then I'll be able to attack the actual problem OP was talking about next time I got some time to put on this.


EDIT: Fixed. ChatGPT is so nice.

Now I'm getting the natural modes correctly. This is the output for n=4.

The number of bodies is: 4

The mass of each body from 0 to n is:
[1, 1, 1, 1] [kg]

The stiffness of each spring from 0 to n+1 is:
[1, 1, 1, 1, 1] [N/m]

The inertia matrix is: Matrix([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]]) [kg]

The stiffness matrix is: Matrix([
[-2, 1, 0, 0],
[1, -2, 1, 0],
[0, 1, -2, 1],
[0, 0, 1, -2]]) [N/m]

The natural angular velocities are:
[0.618, 1.175, 1.618, 1.902] [rad/s]

Normalized natural modes (m) (Column = vector): Matrix([
[0.371, 0.601, 0.601, 0.371],
[0.601, 0.371, -0.371, -0.601],
[0.601, -0.371, -0.371, 0.601],
[0.37174, -0.601, 0.601, -0.371]]) [m]
 
Last edited:
  • #10
I wrote a very sketchy code that works fine when you're testing a single frequency. The generalization to calculate multiple frequencies takes forever and also seems to be imprecise because some of the output data doesn't make much sense to me.

I have simulated 100 s and then applied FFT to the signal ignoring about the initial half of the time so the analysis is performed on what would be the stationary response.

1719844979222.png


That's for a system with the conditions OP mentioned in the beginning:

n = 2
m = [0.001, 0.09] # example masses in kg
k = [1872.5, 1777.3, 130] # example spring constants in N/m
c = [2.85, 0.5, 1] # example damper constants in Ns/m
The input force has an amplitude of 1.63 N as OP shared.

Such a system has natural frequencies at [17.101, 304.458] Hz.

Resonance at 17 Hz can be seen in the graph. The other one at 304 Hz cannot be very well appreciated. I believe this is because the input force is very small compared to the stiffness and dissipative terms.

I would keep refining the code but it takes so long to calculate... In case anyone wants to give it a try, here it is:
*Note: I once tried copying a code I wrote here and paste it in a different computer. It didn't work correctly. If you find yourself in the same situation, you can paste it in ChatGPT so it writes it and copy paste it again from its output.


N bodies system connected by springs:
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 18 20:16:23 2024

@author: Juanda
"""
import sympy as sp
import numpy as np
from IPython.display import display, Math
import matplotlib.pyplot as plt
from scipy.integrate import odeint



# # Number of bodies
# n = 2

# Define symbolic positions
# x = sp.symbols(f'x:{n}')
# v = sp.symbols(f'v:{n}')
# latex_position = sp.latex(x)
# display(Math(latex_position))

# Define numeric or symbolic masses and spring constants
# m = sp.symbols(f'm:{n}')
# k = sp.symbols(f'k:{n+1}')
# c = sp.symbols(f'c:{n+1}')

# n = 1
# m = [2]  # example masses in kg
# k = [6,4]  # example spring constants in N/m
# c = [1,3]  # example damper constants in Ns/m

n = 2
m = [0.001, 0.09]  # example masses in kg
k = [1872.5, 1777.3, 130]  # example spring constants in N/m
c = [2.85, 0.5, 1]  # example damper constants in Ns/m



# n = 3
# m = [1, 1, 1]  # example masses in kg
# k = [1, 1, 1, 1]  # example spring constants in N/m
# c = [1, 1, 1, 1]  # example damper constants in Ns/m

# n = 4
# m = [1, 1, 1, 1]  # example masses in kg
# k = [1, 1, 1, 1, 1]  # example spring constants in N/m
# c = [1, 1, 1, 1, 1]  # example damper constants in Ns/m


# n = 5
# m = [1, 1, 1, 1, 1]  # example masses in kg
# k = [1, 1, 1, 1, 1, 1]  # example spring constants in N/m
# c = [1, 1, 1, 1, 1, 1]  # example damper constants in Ns/m

# Check if size of m is equal to n
if len(m) != n:
    print(f"Error: Size of vector m ({len(m)}) does not match n ({n}). Stopping execution.")
    raise SystemExit

# Check if size of k is equal to n+1
if len(k) != n+1:
    print(f"Error: Size of vector k ({len(k)}) does not match n+1 ({n+1}). Stopping execution.")
    raise SystemExit

# Check if size of c is equal to n+1
if len(c) != n+1:
    print(f"Error: Size of vector c ({len(c)}) does not match n+1 ({n+1}). Stopping execution.")
    raise SystemExit

# Define time as a symbolic variable
t = sp.symbols('t')

# Define positions as functions of time
x_t = [sp.Function(f'x{i}')(t) for i in range(n)]
# display(Math(sp.latex(x_t)))
# print()

# Define velocities as derivatives of positions with respect to time
v_t = [sp.diff(x_t[i], t) for i in range(n)]
# display(Math(sp.latex(v_t)))
# print()

# Define accelerations as derivatives of positions with respect to time
a_t = [sp.diff(x_t[i], t, t) for i in range(n)]
# display(Math(sp.latex(a_t)))
# print()

# print("The number of bodies is:", n)
# print("The mass of each body from 0 to n is:", m, "[kg]")
# print("The stiffness of each spring from 0 to n+1 is:", k, "[N/m]")



# Kinetic energy
T = sum( 1/2*m[i] * sp.diff(x_t[i], t)**2 for i in range(n))
# print("Kinetic energy T:")
# display(Math(sp.latex(T)))
# print()



# Potential energy
# Spring 1: connected to wall and body 1
U = 0.5 * k[0] * (x_t[0] - 0)**2

# Spring 2 to n-1: connected between two bodies
U += sum(0.5 * k[i+1] * (x_t[i+1] - x_t[i])**2 for i in range(n-1))

# Last spring: connected to last body and wall
U += 0.5 * k[-1] * (x_t[n-1] - 0)**2

# print("Potential energy U:")
# display(Math(sp.latex(U)))
# print()




# Disipative energy
# Damper 1: connected to wall and body 1
D = (0.5 * c[0] * ((sp.diff(x_t[0], t)**2 - 0)))

# Damper 2 to n-1: connected between two bodies
D += sum(0.5 * c[i+1] * (sp.diff(x_t[i+1], t) - sp.diff(x_t[i], t))**2 for i in range(n-1))

# Last damper: connected to last body and wall
D += 0.5 * c[-1] * (sp.diff(x_t[n-1], t) - 0)**2

D = -D

# print("Dissipative term D:")
# display(Math(sp.latex(D)))
# print()

# Define the Lagrangian
L = T - U
# print("Lagrangian:")
# display(Math(sp.latex(L)))
# print()

# d(dL/dv)/dt

# First dL/dv is calculated
# Compute partial derivatives of L with respect to velocities v_i (v_0, v_1, ...)
dLdv = [sp.diff(L, v_t[i]) for i in range(n)]

# Display partial derivatives
# for i in range(n):
#     print(f"Partial derivative of L with respect to v_{i}:")
#     display(Math(sp.latex(dLdv[i])))
#     print()

# Secondly, derive with respect to time
ddLdvdt = [sp.diff(dLdv[i], t) for i in range(n)]
# Display the resulting d(dL/dv)/dt
# for i in range(n):
#     print(f"Time derivative of the Partial derivative of L with respect to v_{i}:")
#     display(Math(sp.latex(ddLdvdt[i])))
#     print()

# Compute partial derivatives of L with respect to positions x_i (x_0, x_1, ...)
dLdx = [sp.diff(L, x_t[i]) for i in range(n)]

# Display partial derivatives of L with respect to position
# for i in range(n):
#     print(f"Partial derivative of L with respect to x_{i}:")
#     display(Math(sp.latex(dLdx[i])))
#     print()


# Compute partial derivatives of D with respect to velocities v_i (v_0, v_1, ...)
dDdv = [sp.diff(D, v_t[i]) for i in range(n)]

# Display partial derivatives of D with respect to velocity
# for i in range(n):
#     print(f"Partial derivative of D with respect to v_{i}:")
#     display(Math(sp.latex(dDdv[i])))
#     print()

# Write the system of differential equations defined through the Lagrangian.
# Sys_eq = sp.Matrix.zeros(n,1)
# for i in range(n):
#     Sys_eq[i,0] = ddLdvdt[i] - dLdx[i] + dDdv[i]
    
# print("System of equations:")
# display(Math(sp.latex(Sys_eq)))
# print()

# I get the coefficients to express it in terms of matrixes
# Create the inertia matrix
M = sp.zeros(n)

# Extract coefficients multiplying a_i (second derivative of positions)
for i in range(n):
    # Find the coefficient of a_i in ddLdvdt[i]
    for j in range(n):
        coeff_ai = sp.expand(ddLdvdt[i]).coeff(a_t[j])
        M[i,j]=np.round(float(coeff_ai), 6)

# print("The inertia matrix is:", M, "[kg]")
# print("Inertia matrix M:")
# display(Math(sp.latex(M)))
# print()



# Extract coefficients multiplying x_i (positions) to find the stiffness matrix
K = sp.zeros(n)
for i in range(n):
    # Find the coefficient of a_i in dLdx[i]
    for j in range(n):
        coeff_xi = sp.expand(dLdx[i]).coeff(x_t[j])
        K[i,j]=coeff_xi

# print("Stiffness matrix K:")
# display(Math(sp.latex(K)))
# print()
# print("The stiffness matrix is:", K, "[N/m]")

# Extract coefficients multiplying v_i (velocities) to find the dissipative matrix
D = sp.zeros(n)
for i in range(n):
    # Find the coefficient of a_i in dLdx[i]
    for j in range(n):
        coeff_vi = sp.expand(dDdv[i]).coeff(v_t[j])
        D[i,j]=coeff_vi

# print("Dissipative matrix D:")
# display(Math(sp.latex(D)))
# print()

# Define the natural frequencies w (rad/s)
r  = sp.symbols('r')

DET = (M*r**2-K).det()
# print("Determinant:")
# display(Math(sp.latex(DET)))
# print()


def process_vector(vector):
    I = sp.I

    # Helper function to treat small real parts as zero
    def zero_small_real_parts(z, tolerance=0.001):
        real_part = sp.re(z)
        imaginary_part = sp.im(z)
        if abs(real_part) < tolerance:
            real_part = 0
        return real_part + I*imaginary_part
    # Adjust the vector for small real parts
    adjusted_vector = [zero_small_real_parts(z) for z in vector]
    # print("adjusted vector")
    # print(adjusted_vector)

    # Check if all values in the vector are imaginary
    if not all (sp.im(z) != 0 and sp.re(z) == 0 for z in adjusted_vector):
        print("Error: Not all values in the vector are purely imaginary.")
        raise SystemExit
        return

    # Extract the imaginary parts of the complex numbers
    imaginary_parts = [sp.im(z) for z in vector]

    # Check for coupled values
    positive_imaginary_parts = set(im for im in imaginary_parts if im > 0)
    negative_imaginary_parts = set(-im for im in imaginary_parts if im < 0)

    if not positive_imaginary_parts <= negative_imaginary_parts:
        print("Error: Coupled values (e.g., -2 and +2) are missing.")
        raise SystemExit
        return

    # Filter out the negative values
    filtered_imaginary_parts = [im for im in imaginary_parts if im > 0]

    # print("Original vector:", vector)
    # print("Imaginary parts:", imaginary_parts)
    # print("Filtered positive imaginary parts:", filtered_imaginary_parts)
    return filtered_imaginary_parts


nat_ang_vel = sp.solve(DET, r)
nat_ang_vel = process_vector(nat_ang_vel)
nat_ang_vel = sorted(nat_ang_vel)
# print("Natural angular velocities (rad/s):")
# display(Math(sp.latex(nat_ang_vel)))
# # print(nat_ang_vel)
# print()
# # print("The natural angular velocities are:", nat_ang_vel, "[rad/s]")


nat_frec = [nat_ang_vel[i]/(2*np.pi) for i in range(n)]
# print("Natural frecuencies (Hz):")
# display(Math(sp.latex(nat_frec)))
# # print(nat_frec)
# print()

# Natural modes
# Declare variales to later solve the system of equations for each natural frecuency
NM = sp.symbols(f'NM_:{n**2-n}')
nat_mod = sp.Matrix.zeros(n)

def create_matrix(vector, n):
    if len(vector) < n * (n - 1):
        raise ValueError("The vector does not contain enough elements to fill the matrix.")

    # Initialize an empty matrix of size (n x n)
    matrix = sp.zeros(n, n)

    # Fill the matrix according to the specified pattern
    vector_index = 0
    for i in range(n):
        matrix[0, i] = 1  # First element in each column is 1
        for j in range(1, n):
            matrix[j, i] = vector[vector_index]
            vector_index += 1

    return matrix

nat_mod = create_matrix(NM, n)
# display(Math(sp.latex(nat_mod)))


def extract_subscript(var):
    # Convert the symbol to a string, split by underscore, and convert the second part to an integer
    return int(str(var).split('_')[1])

eq = sp.Matrix.zeros(n,1)
# nat_ang_vel_vect = sp.Matrix(nat_ang_vel).T
for i in range(n):
    eq = (M*(-1)*nat_ang_vel[i]**2-K)*nat_mod.col(i) #(-1) because the ang_vel is complex and square it should give a neg number.
    eq.row_del(0) #This modifies the original eq without creating a new variable somehow
    sol = sp.solve(eq)
    # Sort the keys of the dictionary based on the subscript
    sorted_keys = sorted(sol.keys(), key=extract_subscript)
    # Create a new ordered dictionary
    sorted_sol = {key: sol[key] for key in sorted_keys}
    sol = list(sorted_sol.values())
    sol.insert(0,1)
    vector = np.array(sol, dtype=float)
    magnitude = np.linalg.norm(vector)
    nat_mod_n = vector / magnitude
    # Natural modes written in columns
    for j in range(n):
        nat_mod[j,i] = np.round(float(nat_mod_n[j]), 6)
        # nat_mod[1,i] = nat_mod_n[1]
    
# print("Normalized natural modes (m) (Column = vector):")
# display(Math(sp.latex(nat_mod)))
# # print(nat_mod)
# print()
# # print("Normalized natural modes (m) (Column = vector):", nat_mod, "[m]")

# Display the movement in the natural modes.
# Create the numerical vector t_num1 to plot the results
def generate_time_vector(w_min, min_points_per_period=20):
    # Calculate the time period for the minimum angular velocity
    T_max = 2 * np.pi / w_min
    
    # Define the time interval
    t_max = 1 * T_max
    
    # Calculate the required number of steps to meet the minimum points per period requirement
    required_steps = int(20 * min_points_per_period)
    t = np.linspace(0, t_max, required_steps)
    
    return t


def calculate_sine_wave(A, w, t_num1):
    # Calculate the sine wave using the given time vector
    A = float(A)
    w = float(w)
    x = A * np.sin(w * t_num1)
    
    return x

# This plot is not as clear because graphs are often on top of another and it's hard to tell what's going on
# def plot_multiple_sine_waves(amplitudes, angular_velocities, t_num1):
#     plt.figure(figsize=(10, 6))
    
#     # Iterate over the amplitudes and angular velocities
#     for i, (A, w) in enumerate(zip(amplitudes, angular_velocities)):
#         x = calculate_sine_wave(A, w, t_num1)
#         plt.plot(t_num1, x, label=f'$x_{i+1}$')
    
#     plt.xlabel('Time (t)')
#     plt.ylabel('Amplitude (x)')
#     plt.title('Superposed Sine Waves')
#     plt.legend()
#     plt.grid(True)
#     plt.show()


def plot_multiple_sine_waves(amplitudes, angular_velocities, t_num1):
    plt.figure(figsize=(10, 6))
    
    # Iterate over the amplitudes and angular velocities
    for i, (A, w) in enumerate(zip(amplitudes, angular_velocities)):
        x = calculate_sine_wave(A, w, t_num1)
        plt.subplot(len(amplitudes), 1, i+1)
        plt.plot(t_num1, x, label=f'$x_{i+1}$')
        plt.xlabel('Time (t)')
        plt.ylabel('Amplitude (x)')
        plt.title(f'Natural Mode {i+1}')
        plt.legend()
        plt.grid(True)
    
    plt.tight_layout()
    plt.show()

# Example usage
# Assuming nat_ang_vel is a list of angular velocities obtained earlier
# nat_ang_vel = [1, 2 * np.pi, 3 * np.pi]  # Example angular velocities
# amplitudes = [1, 2, 1.5]  # Corresponding amplitudes

# Generate the time vector based on the maximum angular velocity
w_min =  np.round(float(nat_ang_vel[0]),6)
t_num1 = generate_time_vector(w_min)

# Plot the multiple sine waves
for i in range(n):
    nat_ang_vel_mode_n = [nat_ang_vel[i]]*n
    plot_multiple_sine_waves(nat_mod.row(i), nat_ang_vel_mode_n, t_num1)

# Quiero poder ver los modos naturales en movimiento. Simular una respuesta con esas frecuencias naturales y sin amortiguamiento.
# He podido hacer la gráfica. Me gustaría verlos en movimiento pero dejo eso para más adelante.

#Create a vector for the input forces on the bodies.
A_f = 1.63 #[N]
w_f = 1 #[Hz]
f = A_f*sp.sin(w_f*(2*np.pi)*t)
f_vec = np.zeros((n,1), dtype=object)
f_vec[0, 0] = f

# Transform the variables so I can use matrix multiplication
x_t_array = np.empty((n, 1), dtype=object)
for i in range(n):
    x_t_array[i, 0] = x_t[i]
v_t_array = np.empty((n, 1), dtype=object)
for i in range(n):
    v_t_array[i, 0] = v_t[i]
a_t_array = np.empty((n, 1), dtype=object)
for i in range(n):
    a_t_array[i, 0] = a_t[i]

Sys_eq = M*x_t_array-D*v_t_array-K*x_t_array-f_vec

# Initial conditions
x0 = np.zeros(n)
v0 = np.zeros(n)
initial_conditions = np.concatenate((x0, v0))

# Define the acceleration to solve the system
a_t = M.inv()*(D*v_t_array + K*x_t_array + f_vec)
a_t_list = [0]*n
for i in range(n):
    a_t_list[i] = a_t[i]
    
# Set up initial conditions
x0 = np.zeros(n)  # Initial positions
v0 = np.zeros(n)  # Initial velocities
initial_conditions = np.concatenate((x0, v0))

# Convert the system of differential equations into first-order form
vars = x_t + v_t

# Lambdify the accelerations (right-hand side of the differential equations)
a_func = [sp.lambdify(vars + [t], a_t_list[i]) for i in range(n)]

# Define a function for odeint that matches its expected format
def equations_of_motion(y, t):
    # Unpack y into positions and velocities
    positions = y[:n]
    velocities = y[n:]
    
    # Evaluate accelerations at current positions and velocities
    accelerations = [a_func[i](*positions, *velocities, t) for i in range(n)]
    
    # Return velocities and accelerations as a flat array
    return np.concatenate((velocities, accelerations))

# Define time points for integration
T = 100
q = 100
time_points = np.linspace(0, T, q)  # Example: from 0 to 10 seconds with 100 points

# Define the sampling frequency
if nat_frec[-1]>=w_f*2*np.pi:
    min_sampling_freq = nat_frec[-1] * 2
else:
    min_sampling_freq = w_f*2*np.pi * 2
        
# Check initial condition
def check_condition(time_points, min_sampling_freq, w_f):
    sampling_interval = time_points[1] - time_points[0]
    lim = 200
    return min_sampling_freq / sampling_interval < lim or (w_f / (2 * np.pi)) ** (-1) / sampling_interval < lim

if check_condition(time_points, min_sampling_freq, w_f):
    print("Error: Initial time interval needs more precision for FFT")
    p = 1
    while p == 1:
        q = 2*q
        time_points = np.linspace(0, T, int(q))
        if check_condition(time_points, min_sampling_freq, w_f):
            print(f"Current q={q}, Error: Time interval needs more precision for FFT")
        else:
            p = 0  # Exit the loop once condition is satisfied

time_points = np.linspace(0, T, 5000)
# Solve the system numerically using odeint
result = odeint(equations_of_motion, initial_conditions, time_points)

# Extract positions and velocities from the result
positions = result[:, :n]
velocities = result[:, n:]


# Plotting positions of all bodies
plt.figure(figsize=(10, 6))

for i in range(n):
    plt.plot(time_points, positions[:, i], label=f'x{i+1}(t)')

plt.title('Positions of Bodies')
plt.xlabel('Time')
plt.ylabel('Position')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# Plotting velocities of all bodies
plt.figure(figsize=(10, 6))

for i in range(n):
    plt.plot(time_points, velocities[:, i], label=f'v{i+1}(t)')

plt.title('Velocities of Bodies')
plt.xlabel('Time')
plt.ylabel('Velocity')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()



# Compute FFT for each body's position
fft_results = np.fft.fft(positions, axis=0)/(len(time_points)/2)
fft_freq = np.fft.fftfreq(len(time_points), d=(time_points[1] - time_points[0]))

# Plotting FFT results
plt.figure(figsize=(12, 8))
for i in range(n):
    plt.subplot(n, 1, i + 1)
    plt.plot(fft_freq[:len(fft_freq)//2], np.abs(fft_results[:len(fft_freq)//2, i]))  # Plotting half of spectrum due to symmetry
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Amplitude')
    plt.title(f'FFT of x{i + 1}(t)')
    plt.grid(True)

plt.tight_layout()
plt.show()


# Define the duration in seconds to ignore at the beginning
ignore_initial_seconds = T/1.5

# Calculate the number of initial points to ignore based on the sampling frequency
ignore_initial_points = int(ignore_initial_seconds / (time_points[1] - time_points[0]))

# Adjust positions and time_points arrays to ignore initial points
positions_adjusted = positions[ignore_initial_points:, :]
time_points_adjusted = time_points[ignore_initial_points:]

# Compute FFT for each body's adjusted position
fft_results_adjusted = np.fft.fft(positions_adjusted, axis=0)/(len(time_points_adjusted)/2)
fft_freq_adjusted = np.fft.fftfreq(len(time_points_adjusted), d=(time_points_adjusted[1] - time_points_adjusted[0]))



# Plotting adjusted FFT results
plt.figure(figsize=(12, 8))
for i in range(n):
    plt.subplot(n, 1, i + 1)
    plt.plot(fft_freq_adjusted[:len(fft_freq_adjusted)//2], np.abs(fft_results_adjusted[:len(fft_freq_adjusted)//2, i]))  # Plotting half of spectrum due to symmetry
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Amplitude')
    plt.title(f'FFT of x{i + 1}(t) (Ignoring initial {ignore_initial_seconds} seconds)')
    plt.grid(True)

plt.tight_layout()
plt.show()






def calculate_fft_amplitudes(w_f, ignore_initial_seconds_ratio=1.2):
    # Create a vector for the input forces on the bodies
    A_f = 1.63  # [N]
    f = A_f * sp.sin(w_f * (2 * np.pi) * t)
    f_vec = np.zeros((n, 1), dtype=object)
    f_vec[0, 0] = f

    # Sys_eq = M * x_t_array - D * v_t_array - K * x_t_array - f_vec

    # Initial conditions
    x0 = np.zeros(n)
    v0 = np.zeros(n)
    initial_conditions = np.concatenate((x0, v0))

    # Define the acceleration to solve the system
    a_t = M.inv() * (D * v_t_array + K * x_t_array + f_vec)
    a_t_list = [0] * n
    for i in range(n):
        a_t_list[i] = a_t[i]

    # Lambdify the accelerations (right-hand side of the differential equations)
    a_func = [sp.lambdify(vars + [t], a_t_list[i]) for i in range(n)]

    # Define a function for odeint that matches its expected format
    def equations_of_motion(y, t):
        # Unpack y into positions and velocities
        positions = y[:n]
        velocities = y[n:]

        # Evaluate accelerations at current positions and velocities
        accelerations = [a_func[i](*positions, *velocities, t) for i in range(n)]

        # Return velocities and accelerations as a flat array
        return np.concatenate((velocities, accelerations))

    # Solve the system numerically using odeint
    result = odeint(equations_of_motion, initial_conditions, time_points)

    # Extract positions and velocities from the result
    positions = result[:, :n]
    # velocities = result[:, n:]

    # Define the duration in seconds to ignore at the beginning
    ignore_initial_seconds = T / ignore_initial_seconds_ratio

    # Calculate the number of initial points to ignore based on the sampling frequency
    ignore_initial_points = int(ignore_initial_seconds / (time_points[1] - time_points[0]))

    # Adjust positions and time_points arrays to ignore initial points
    positions_adjusted = positions[ignore_initial_points:, :]
    time_points_adjusted = time_points[ignore_initial_points:]

    # Compute FFT for each body's adjusted position
    fft_results_adjusted = np.fft.fft(positions_adjusted, axis=0) / (len(time_points_adjusted) / 2)
    fft_freq_adjusted = np.fft.fftfreq(len(time_points_adjusted), d=(time_points_adjusted[1] - time_points_adjusted[0]))

    # Extract the maximum amplitude in the FFT results for each body
    max_amplitudes = np.max(np.abs(fft_results_adjusted[:len(fft_freq_adjusted) // 2, :]), axis=0)

    return max_amplitudes


# Define the range of frequencies
frequencies = np.linspace(1, 350, 350)  # From 1 to 1000 (excluding 0)

# Initialize arrays to store the maximum amplitudes for x1 and x2
max_amplitudes_x1 = []
max_amplitudes_x2 = []

for w_f in frequencies:
    max_amplitudes = calculate_fft_amplitudes(w_f)
    max_amplitudes_x1.append(max_amplitudes[0])
    max_amplitudes_x2.append(max_amplitudes[1])

# Convert lists to arrays
max_amplitudes_x1 = np.array(max_amplitudes_x1)
max_amplitudes_x2 = np.array(max_amplitudes_x2)


plt.figure(figsize=(10, 6))

# Plot for x1
plt.subplot(2, 1, 1)
plt.loglog(frequencies, max_amplitudes_x1, label='$x_1$', marker='o')
plt.xlabel('Input Frequency (Hz)')
plt.ylabel('Maximum Amplitude')
plt.title('Maximum Amplitude vs Input Frequency for $x_1$')
plt.grid(True)
plt.legend()

# Plot for x2
plt.subplot(2, 1, 2)
plt.loglog(frequencies, max_amplitudes_x2, label='$x_2$', marker='o')
plt.xlabel('Input Frequency (Hz)')
plt.ylabel('Maximum Amplitude')
plt.title('Maximum Amplitude vs Input Frequency for $x_2$')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()
 
  • #11
I wasn't sure the code was right so I tested it with different initial conditions to check the hypothesis I mentioned.

Hypothesis:
Juanda said:
Resonance at 17 Hz can be seen in the graph. The other one at 304 Hz cannot be very well appreciated. I believe this is because the input force is very small compared to the stiffness and dissipative terms.

So the new initial conditions for this test are:
m = [0.001, 0.05] # example masses in kg
k = [10,50, 250] # example spring constants in N/m
c = [0.01, 0.01, 0.01] # example damper constants in Ns/m

Such a system has the natural frequencies (Hz):
[11.357, 39.27]

The input force is ##f(t) = A_f \sin(2\pi ft)## so peaks should appear within the range from 0 to 60 Hz which is what I calculated and this is what I obtained.
Where, ##A_f = 50 [N]##, ##f## will range from 0 to 60 Hz and frequencies will be calculated Hertz by Hertz.

1719952364632.png



So the results make physical sense. It's no proof that the code is OK but it gives it some confidence. I'll admit though that I had to manually set some parameters so the calculation could be done in a reasonable amount of time but the point stands. The results from it should be trust worth it.
By the way, the code reveals the expected result of the natural frequencies being dragged to the left due to the damping. For example, the peak amplitude happens at 38 Hz instead of 39 Hz which is the closest value to the highest natural frequency in the simulation. The results are close enough for this to be a coincidence, but I felt it was interesting to point out.

The outlier points from the original plot in post #10 are still a little bit of a mystery. Maybe the computer was running out of memory? Who knows... The code definitely isn't optimal.
I wish I could spend an indefinite amount of time on it but I think I am done with it for now.
 
Back
Top