振動工学Pythonコード

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 定数
m = 1.0  # 質量
k = 10.0 # バネ定数
c = 0.5  # 粘性係数
F0 = 1.0 # 外力の振幅
omega = 1.0 # 外力の角周波数

# 運動方程式
def model(y, t):
    x, v = y
    dydt = [v, (F0 * np.sin(omega * t) - k * x - c * v) / m]
    return dydt

# 初期条件
y0 = [0.0, 0.0]

# 時間配列
t = np.linspace(0, 20, 500)

# 数値解法
sol = odeint(model, y0, t)

# プロット
plt.figure()
plt.plot(t, sol[:, 0], 'b-', label='Displacement')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.title('Forced Vibration Response')
plt.legend()
plt.grid()
plt.show()

omega_values = np.linspace(0.5, 2.5, 100)
displacements = []

for omega in omega_values:
    sol = odeint(model, y0, t)
    displacements.append(sol[-1, 0])  # 最終時刻の変位

# プロット
plt.figure()
plt.plot(omega_values, displacements, 'r-')
plt.xlabel('Frequency (ω)')
plt.ylabel('Displacement')
plt.title('Frequency Response')
plt.grid()
plt.show()

from scipy import signal

# 系の伝達関数
num = [F0]
den = [m, c, k]
sys = signal.TransferFunction(num, den)

# 周波数応答
w, mag, phase = signal.bode(sys)

# プロット
plt.figure()
plt.subplot(2, 1, 1)
plt.semilogx(w, mag)
plt.xlabel('Frequency [rad/s]')
plt.ylabel('Magnitude [dB]')
plt.grid()

plt.subplot(2, 1, 2)
plt.semilogx(w, phase)
plt.xlabel('Frequency [rad/s]')
plt.ylabel('Phase [degrees]')
plt.grid()

plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 定数
m1, m2 = 1.0, 1.0  # 質量
k1, k2, k3 = 10.0, 10.0, 10.0  # バネ定数
c1, c2 = 0.5, 0.5  # 粘性係数

# 運動方程式
def model(Y, t):
    x1, v1, x2, v2 = Y
    dydt = [v1, 
            (k2 * (x2 - x1) - c2 * (v2 - v1) - k1 * x1 - c1 * v1) / m1,
            v2,
            (F0 * np.sin(omega * t) - k2 * (x2 - x1) - c2 * (v2 - v1) - k3 * x2 - c2 * v2) / m2]
    return dydt

# 初期条件
Y0 = [0.0, 0.0, 0.0, 0.0]

# 時間配列
t = np.linspace(0, 20, 500)

# 数値解法
sol = odeint(model, Y0, t)

# プロット
plt.figure()
plt.plot(t, sol[:, 0], 'b-', label='Displacement of m1')
plt.plot(t, sol[:, 2], 'r-', label='Displacement of m2')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.title('Displacements of a 2-DOF System')
plt.legend()
plt.grid()
plt.show()

from scipy.linalg import eig

# 行列の定義
K = np.array([[k1 + k2, -k2],
              [-k2, k2 + k3]])

M = np.array([[m1, 0],
              [0, m2]])

# 固有値と固有ベクトルの計算
eigenvalues, eigenvectors = eig(K, M)

# 固有振動数
natural_frequencies = np.sqrt(np.abs(eigenvalues))

# プロット
plt.figure()
plt.bar(range(len(natural_frequencies)), natural_frequencies)
plt.xlabel('Mode')
plt.ylabel('Natural Frequency (rad/s)')
plt.title('Natural Frequencies of the 2-DOF System')
plt.grid()
plt.show()

# モードの表示
print("Eigenvectors (Mode Shapes):")
print(eigenvectors)

# 更新された運動方程式
def forced_model(Y, t):
    x1, v1, x2, v2 = Y
    dydt = [v1, 
            (k2 * (x2 - x1) - c2 * (v2 - v1) - k1 * x1 - c1 * v1) / m1,
            v2,
            (F0 * np.sin(omega * t) - k2 * (x2 - x1) - c2 * (v2 - v1) - k3 * x2 - c2 * v2) / m2]
    return dydt

# 初期条件
Y0 = [0.0, 0.0, 0.0, 0.0]

# 外力の定義
F0 = 1.0
omega = 1.0

# 数値解法
sol = odeint(forced_model, Y0, t)

# プロット
plt.figure()
plt.plot(t, sol[:, 0], 'b-', label='Displacement of m1')
plt.plot(t, sol[:, 2], 'r-', label='Displacement of m2')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.title('Forced Vibration Response of a 2-DOF System')
plt.legend()
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定数
L = 1.0  # 棒の長さ
n_modes = 5  # モード数

# 固有振動数の計算
def eigenfrequencies(n_modes, L):
    return [(n*np.pi/L)**2 for n in range(1, n_modes+1)]

# 固有振動数
frequencies = eigenfrequencies(n_modes, L)

# プロット
plt.figure()
plt.plot(range(1, n_modes+1), frequencies, 'o-')
plt.xlabel('Mode number')
plt.ylabel('Frequency squared (rad^2/s^2)')
plt.title('Eigenfrequencies of a Uniform Rod')
plt.grid()
plt.show()

from scipy.linalg import eig

# 定数
J = 1.0  # 軸の断面二次モーメント
G = 1.0  # 剛性係数

# モデルの定義
def torsional_eigenvalues(n_modes, L, J, G):
    return [(n*np.pi/L)**2 * J * G for n in range(1, n_modes+1)]

# 固有振動数
torsional_frequencies = torsional_eigenvalues(n_modes, L, J, G)

# プロット
plt.figure()
plt.plot(range(1, n_modes+1), torsional_frequencies, 'o-')
plt.xlabel('Mode number')
plt.ylabel('Torsional Frequency squared (rad^2/s^2)')
plt.title('Torsional Eigenfrequencies of a Shaft')
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定数
L = 1.0  # はりの長さ
EI = 1.0  # 曲げ剛性
n_modes = 5  # モード数

# モード形状の計算
x = np.linspace(0, L, 1000)
def mode_shape(n, x, L):
    return np.sin(n * np.pi * x / L)

# モード形状のプロット
plt.figure()
for n in range(1, n_modes+1):
    plt.plot(x, mode_shape(n, x, L), label=f'Mode {n}')
plt.xlabel('Position along the beam')
plt.ylabel('Displacement')
plt.title('Mode Shapes of a Beam')
plt.legend()
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定数
a = 1.0  # 平板の一辺の長さ
b = 1.0  # 平板の他の辺の長さ
D = 1.0  # 曲げ剛性
n_modes = 3  # モード数

# モード形状の計算
x = np.linspace(0, a, 100)
y = np.linspace(0, b, 100)
X, Y = np.meshgrid(x, y)

def mode_shape(m, n, X, Y, a, b):
    return np.sin(m * np.pi * X / a) * np.sin(n * np.pi * Y / b)

# モード形状のプロット
plt.figure()
for m in range(1, n_modes+1):
    for n in range(1, n_modes+1):
        Z = mode_shape(m, n, X, Y, a, b)
        plt.contourf(X, Y, Z, levels=20, cmap='viridis')
        plt.title(f'Mode Shape (m={m}, n={n})')
        plt.colorbar()
        plt.xlabel('x')
        plt.ylabel('y')
        plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定数
I = 1.0  # 慣性モーメント
k = 1.0  # 剛性係数
n_modes = 5  # モード数

# 固有振動数の計算
def rotational_frequencies(n_modes, I, k):
    return [(n*np.pi)**2 * I / k for n in range(1, n_modes+1)]

# 固有振動数
rotational_frequencies = rotational_frequencies(n_modes, I, k)

# プロット
plt.figure()
plt.plot(range(1, n_modes+1), rotational_frequencies, 'o-')
plt.xlabel('Mode number')
plt.ylabel('Rotational Frequency squared (rad^2/s^2)')
plt.title('Rotational Eigenfrequencies')
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定数
J = 1.0  # 慣性モーメント
G = 1.0  # 剛性係数
n_modes = 5  # モード数

# 固有振動数の計算
def torsional_frequencies(n_modes, J, G):
    return [(n*np.pi)**2 * J / G for n in range(1, n_modes+1)]

# 固有振動数
torsional_frequencies = torsional_frequencies(n_modes, J, G)

# プロット
plt.figure()
plt.plot(range(1, n_modes+1), torsional_frequencies, 'o-')
plt.xlabel('Mode number')
plt.ylabel('Torsional Frequency squared (rad^2/s^2)')
plt.title('Torsional Eigenfrequencies')
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# ダッフリング振動子のモデル
def duffing_oscillator(y, t, delta, gamma, omega):
    x, v = y
    dydt = [v, -delta*v - x**3 + gamma*np.cos(omega*t)]
    return dydt

# パラメータ設定
delta = 0.2
gamma = 0.5
omega = 1.0

# 初期条件
y0 = [1.0, 0.0]
t = np.linspace(0, 20, 1000)

# 微分方程式を解く
solution = odeint(duffing_oscillator, y0, t, args=(delta, gamma, omega))

# 結果のプロット
plt.plot(t, solution[:, 0])
plt.title('Duffing Oscillator - Nonlinear Free Vibration')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.grid(True)
plt.show()

# 非線形強制振動のモデル
def forced_duffing_oscillator(y, t, delta, gamma, omega, F0, omega_ext):
    x, v = y
    dydt = [v, -delta*v - x**3 + F0*np.cos(omega_ext*t)]
    return dydt

# パラメータ設定
delta = 0.2
gamma = 0.5
omega = 1.0
F0 = 0.5
omega_ext = 1.2

# 初期条件
y0 = [1.0, 0.0]
t = np.linspace(0, 20, 1000)

# 微分方程式を解く
solution = odeint(forced_duffing_oscillator, y0, t, args=(delta, gamma, omega, F0, omega_ext))

# 結果のプロット
plt.plot(t, solution[:, 0])
plt.title('Forced Duffing Oscillator - Nonlinear Forced Vibration')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.grid(True)
plt.show()

# 実際の機械システムの例として、単純な非線形ダンパー付き振動子のモデル
def nonlinear_damper_oscillator(y, t, k, c, m, F0, omega_ext):
    x, v = y
    dydt = [v, (F0*np.cos(omega_ext*t) - c*v - k*x**3) / m]
    return dydt

# パラメータ設定
k = 1.0
c = 0.5
m = 1.0
F0 = 0.5
omega_ext = 1.2

# 初期条件
y0 = [0.5, 0.0]
t = np.linspace(0, 20, 1000)

# 微分方程式を解く
solution = odeint(nonlinear_damper_oscillator, y0, t, args=(k, c, m, F0, omega_ext))

# 結果のプロット
plt.plot(t, solution[:, 0])
plt.title('Nonlinear Damper Oscillator - Practical Mechanical System')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.grid(True)
plt.show()

from scipy.signal import correlate, welch

# サンプルデータの生成
np.random.seed(0)
x = np.random.randn(1000)

# 相関関数
corr = correlate(x, x, mode='full')
lags = np.arange(-len(x) + 1, len(x))

plt.plot(lags, corr)
plt.title('Autocorrelation Function')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.grid(True)
plt.show()

# スペクトル密度
frequencies, power_spectrum = welch(x)

plt.semilogy(frequencies, power_spectrum)
plt.title('Power Spectral Density')
plt.xlabel('Frequency')
plt.ylabel('Power/Frequency')
plt.grid(True)
plt.show()

from scipy.integrate import odeint

def linear_oscillator(y, t, k, c, m, F0, omega_ext):
    x, v = y
    dydt = [v, (F0*np.cos(omega_ext*t) - c*v - k*x) / m]
    return dydt

# パラメータ設定
k = 1.0
c = 0.1
m = 1.0
F0 = 1.0
omega_ext = 2.0

# 初期条件
y0 = [0.0, 0.0]
t = np.linspace(0, 10, 1000)

# 微分方程式を解く
solution = odeint(linear_oscillator, y0, t, args=(k, c, m, F0, omega_ext))

# 結果のプロット
plt.plot(t, solution[:, 0])
plt.title('Linear Oscillator with Irregular Vibration')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.grid(True)
plt.show()

def parametric_excitation(y, t, k, c, m, gamma, omega_ext):
    x, v = y
    dydt = [v, (gamma * np.cos(omega_ext * t) - k * x - c * v) / m]
    return dydt

# パラメータ設定
k = 1.0
c = 0.1
m = 1.0
gamma = 1.0
omega_ext = 2.0

# 初期条件
y0 = [0.0, 0.0]
t = np.linspace(0, 10, 1000)

# 微分方程式を解く
solution = odeint(parametric_excitation, y0, t, args=(k, c, m, gamma, omega_ext))

# 結果のプロット
plt.plot(t, solution[:, 0])
plt.title('Parametric Excitation Oscillator')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.grid(True)
plt.show()

def lorentz_attractor(state, t, sigma, beta, rho):
    x, y, z = state
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]

# パラメータ設定
sigma = 10.0
beta = 8.0 / 3.0
rho = 28.0

# 初期条件
state0 = [1.0, 0.0, 0.0]
t = np.linspace(0, 30, 10000)

# 微分方程式を解く
solution = odeint(lorentz_attractor, state0, t, args=(sigma, beta, rho))

# 結果のプロット
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot(solution[:, 0], solution[:, 1], solution[:, 2])
ax.set_title('Lorentz Attractor - Chaotic Oscillator')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

def critically_damped_oscillator(y, t, c, m, k):
    x, v = y
    dydt = [v, (-c*v - k*x) / m]
    return dydt

# パラメータ設定
c = 2.0  # 臨界減衰係数
m = 1.0
k = 1.0

# 初期条件
y0 = [1.0, 0.0]
t = np.linspace(0, 10, 1000)

# 微分方程式を解く
solution = odeint(critically_damped_oscillator, y0, t, args=(c, m, k))

# 結果のプロット
plt.plot(t, solution[:, 0])
plt.title('Critically Damped Oscillator')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.grid(True)
plt.show()

import numpy as np

# 減衰係数と固有振動数の設定
c = 0.5
m = 1.0
k = 1.0

# 対数減衰率の計算
wn = np.sqrt(k / m)  # 固有振動数
zeta = c / (2 * np.sqrt(k * m))  # 減衰比
log_decrement = (1 / zeta) * np.log(1 / (1 - zeta**2)**0.5)

print(f'Logarithmic Decrement: {log_decrement:.2f}')

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

# パラメータ設定
N = 4  # 自由度の数

# 質量行列と剛性行列の定義（例として対角行列を使用）
M = np.diag([1.0] * N)  # 質量行列
K = np.diag([2.0] * N)  # 剛性行列
for i in range(N-1):
    K[i, i+1] = -1.0
    K[i+1, i] = -1.0

# 固有値問題を解く
eigenvalues, eigenvectors = eigh(K, M)

# 固有振動数
natural_frequencies = np.sqrt(eigenvalues)

# 結果の表示
print("Natural Frequencies:", natural_frequencies)
print("Mode Shapes:\n", eigenvectors)

# モードプロット
fig, axs = plt.subplots(N, 1, figsize=(10, 2*N))

for i in range(N):
    axs[i].plot(eigenvectors[:, i], marker='o', linestyle='-', label=f'Mode {i+1}')
    axs[i].set_title(f'Mode {i+1} Shape')
    axs[i].set_xlabel('DOF')
    axs[i].set_ylabel('Displacement')
    axs[i].grid(True)
    axs[i].legend()

plt.tight_layout()
plt.show()

def calculate_beam_frequencies(n_modes, L):
    frequencies = []
    for n in range(1, n_modes + 1):
        frequency = n * np.pi / L
        frequencies.append(frequency)
    return np.array(frequencies)

# 固有振動数の計算
n_modes = 5
frequencies = calculate_beam_frequencies(n_modes, L)
print("Natural Frequencies of the Beam:", frequencies)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# はりのパラメータ設定
L = 1.0  # はりの長さ
E = 1.0  # ヤング率
I = 1.0  # 断面2次モーメント
rho = 1.0  # 密度
A = 1.0  # 断面積

# 運動方程式の定義
def bending_vibration(w, t, n):
    x = np.linspace(0, L, len(w))
    dw_dt = np.zeros_like(w)
    for i in range(1, len(w)-1):
        dw_dt[i] = (E * I / (rho * A)) * (w[i-1] - 2*w[i] + w[i+1])
    return dw_dt

# 初期条件と時間設定
x = np.linspace(0, L, 100)
w0 = np.sin(np.pi * x / L)  # 初期の変位（例: 最初のモード）
t = np.linspace(0, 10, 200)

# 微分方程式を解く
solution = odeint(bending_vibration, w0, t, args=(1,))

# 結果のプロット
plt.imshow(solution.T, extent=[t[0], t[-1], x[0], x[-1]], origin='lower', aspect='auto')
plt.title('Bending Vibration of the Beam')
plt.xlabel('Time')
plt.ylabel('Position along the beam')
plt.colorbar(label='Displacement')
plt.show()

import numpy as np

def calculate_beam_frequencies(n_modes, L):
    frequencies = []
    for n in range(1, n_modes + 1):
        frequency = (n * np.pi)**2 * np.sqrt(E * I / (rho * A)) / L**2
        frequencies.append(frequency)
    return np.array(frequencies)

# 固有振動数の計算
n_modes = 5
frequencies = calculate_beam_frequencies(n_modes, L)
print("Natural Frequencies of the Beam:", frequencies)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 平板のパラメータ設定
Lx = 1.0  # x方向の長さ
Ly = 1.0  # y方向の長さ
E = 1.0   # ヤング率
h = 0.1   # 板の厚さ
nu = 0.3  # ポアソン比
D = E * h**3 / (12 * (1 - nu**2))
rho = 1.0 # 密度

# 運動方程式の定義
def bending_vibration_2D(w, t, n):
    x = np.linspace(0, Lx, len(w))
    y = np.linspace(0, Ly, len(w[0]))
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    
    dw_dt = np.zeros_like(w)
    for i in range(1, len(x)-1):
        for j in range(1, len(y)-1):
            dw_dt[i, j] = (D / (rho * h)) * (
                (w[i-1, j] - 2*w[i, j] + w[i+1, j]) / dx**2 +
                (w[i, j-1] - 2*w[i, j] + w[i, j+1]) / dy**2
            )
    return dw_dt.flatten()

# 初期条件と時間設定
x = np.linspace(0, Lx, 50)
y = np.linspace(0, Ly, 50)
X, Y = np.meshgrid(x, y)
w0 = np.sin(np.pi * X) * np.sin(np.pi * Y)
w0 = w0.flatten()
t = np.linspace(0, 10, 200)

# 微分方程式を解く
solution = odeint(bending_vibration_2D, w0, t)
solution = solution.reshape(-1, len(x), len(y))

# 結果のプロット
plt.imshow(solution[-1], extent=[x[0], x[-1], y[0], y[-1]], origin='lower', aspect='auto')
plt.title('Bending Vibration of the Plate')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar(label='Displacement')
plt.show()

def calculate_plate_frequencies(n_modes_x, n_modes_y, Lx, Ly):
    frequencies = []
    for nx in range(1, n_modes_x + 1):
        for ny in range(1, n_modes_y + 1):
            frequency = np.sqrt(
                (nx**2 * np.pi**2 + ny**2 * np.pi**2) * D / (rho * h)
            )
            frequencies.append(frequency)
    return np.array(frequencies)

# 固有振動数の計算
n_modes_x = 3
n_modes_y = 3
frequencies = calculate_plate_frequencies(n_modes_x, n_modes_y, Lx, Ly)
print("Natural Frequencies of the Plate:", frequencies)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 系のパラメータ設定
m = 1.0   # 質量
k = 1.0   # バネ定数
c = 0.2   # 線形ダンパー定数
F0 = 1.0  # 強制力の振幅
omega = 2.0  # 強制力の周波数（固有振動数の2倍）

# 非線形ダンパーの定義
def nonlinear_damper(v):
    return 0.2 * v**3

# 運動方程式の定義
def equation(y, t, omega, F0):
    x, v = y
    dxdt = v
    dvdt = (F0 * np.sin(omega * t) - k * x - c * v - nonlinear_damper(v)) / m
    return [dxdt, dvdt]

# 初期条件と時間設定
y0 = [0.0, 0.0]
t = np.linspace(0, 50, 1000)

# 微分方程式を解く
solution = odeint(equation, y0, t, args=(omega, F0))
x = solution[:, 0]

# 結果のプロット
plt.plot(t, x)
plt.title('Secondary Resonance in Nonlinear System')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Parameters
L = 1.0  # Length of the rod (m)
E = 2.1e11  # Young's modulus (Pa)
rho = 7800  # Density (kg/m^3)
A = 1e-4  # Cross-sectional area (m^2)

# Calculate natural frequencies
n = np.arange(1, 6)  # Mode numbers
omega_n = n * np.pi / L * np.sqrt(E / (rho * A))

# Plot
plt.figure()
plt.plot(n, omega_n, 'o-')
plt.xlabel('Mode Number (n)')
plt.ylabel('Natural Frequency (rad/s)')
plt.title('Natural Frequencies of Vertical Vibrations of a Rod')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Parameters
L = 1.0  # Length of the beam (m)
E = 2.1e11  # Young's modulus (Pa)
I = 1e-6  # Second moment of area (m^4)
rho = 7800  # Density (kg/m^3)
A = 1e-4  # Cross-sectional area (m^2)

# Calculate natural frequencies
n = np.arange(1, 6)  # Mode numbers
omega_n = (n * np.pi / L) ** 2 * np.sqrt(E * I / (rho * A))

# Plot
plt.figure()
plt.plot(n, omega_n, 'o-')
plt.xlabel('Mode Number (n)')
plt.ylabel('Natural Frequency (rad/s)')
plt.title('Natural Frequencies of Bending Vibrations of a Beam')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Parameters
Lx = 1.0  # Length of the plate (m)
Ly = 1.0  # Width of the plate (m)
E = 2.1e11  # Young's modulus (Pa)
nu = 0.3  # Poisson's ratio
rho = 7800  # Density (kg/m^3)
h = 0.01  # Thickness of the plate (m)

# Stiffness coefficient
D = E * h**3 / (12 * (1 - nu**2))

# Mode numbers
n = np.arange(1, 4)
m = np.arange(1, 4)

# Calculate natural frequencies
omega_n = np.zeros((len(n), len(m)))
for i, ni in enumerate(n):
    for j, mj in enumerate(m):
        omega_n[i, j] = np.pi * np.sqrt(D * (ni**2 / Lx**2 + mj**2 / Ly**2) / (rho * h))

# Plot
X, Y = np.meshgrid(n, m)
plt.figure()
contour = plt.contourf(X, Y, omega_n.T, cmap='viridis')
plt.colorbar(contour, label='Natural Frequency (rad/s)')
plt.xlabel('Mode Number (n)')
plt.ylabel('Mode Number (m)')
plt.title('Natural Frequencies of Bending Vibrations of a Plate')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Parameters
L = 1.0  # Length of the beam (m)
E = 2.1e11  # Young's modulus (Pa)
I = 1e-6  # Moment of inertia (m^4)
rho = 7800  # Density (kg/m^3)
A = 1e-4  # Cross-sectional area (m^2)

# Transfer Matrix for a Beam Element
def transfer_matrix(length, E, I, rho, A):
    k = E * I / length**3  # Stiffness per unit length
    m = rho * A * length  # Mass per unit length

    # Element transfer matrix
    T = np.array([
        [1, length, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 1, length],
        [0, 0, 0, 1]
    ])
    
    return T

# Formulate transfer matrix for the entire beam
T_total = transfer_matrix(L, E, I, rho, A)

# Calculate natural frequencies
eigenvalues, _ = np.linalg.eig(T_total)
natural_frequencies = np.sqrt(np.abs(eigenvalues))

# Plot natural frequencies
plt.figure()
plt.plot(natural_frequencies, 'o-')
plt.xlabel('Mode Number')
plt.ylabel('Natural Frequency (rad/s)')
plt.title('Natural Frequencies of Simply Supported Beam')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Mass and stiffness matrices
M = np.array([[2, 0], [0, 1]])  # Mass matrix
K = np.array([[4, -2], [-2, 2]])  # Stiffness matrix

# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(M) @ K)

# Natural frequencies (rad/s)
natural_frequencies = np.sqrt(np.abs(eigenvalues))

# Mode shapes
mode_shapes = eigenvectors

# Print natural frequencies and mode shapes
print("Natural Frequencies (rad/s):", natural_frequencies)
print("Mode Shapes:\n", mode_shapes)

# Plot mode shapes
plt.figure()
for i in range(mode_shapes.shape[1]):
    plt.plot(mode_shapes[:, i], label=f'Mode {i+1}')

plt.xlabel('DOF')
plt.ylabel('Mode Shape')
plt.title('Mode Shapes of 2-DOF System')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
from scipy.linalg import eig

# Function to compute the mass matrix
def eval_Mmatrix(m_vec):
    return np.diag(m_vec)

# Function to compute the stiffness matrix
def eval_Kmatrix(k_vec):
    n = len(k_vec)
    K = np.zeros((n, n))
    for i in range(n):
        if i == 0:
            K[i, i] = k_vec[i]
        else:
            K[i-1:i+1, i-1:i+1] += np.array([[1, -1], [-1, 1]]) * k_vec[i]
    return K

# Define mass and stiffness vectors
m_vec = np.ones(100)
k_vec = np.arange(1, 101) * 10**3

# Compute mass and stiffness matrices
M = eval_Mmatrix(m_vec)
K = eval_Kmatrix(k_vec)

# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = eig(K, M)

# Compute natural frequencies
wn = np.sqrt(np.abs(eigenvalues))  # Angular frequencies
fn = wn / (2 * np.pi)  # Natural frequencies in Hz

print("Natural Frequencies (Hz):", fn)

# Function to compute the Modal Assurance Criterion (MAC) matrix
def compute_MAC(V):
    num_modes = V.shape[1]
    MAC = np.zeros((num_modes * num_modes, 1))
    index = 0
    for i in range(num_modes):
        for j in range(num_modes):
            MAC[index] = np.dot(V[:, i].T, V[:, j])
            index += 1
    return MAC

# Compute MAC values and reshape them
MAC = compute_MAC(eigenvectors)
MAC_matrix = MAC.reshape((eigenvectors.shape[1], eigenvectors.shape[1]))

print("MAC Matrix:")
print(MAC_matrix)

import numpy as np

# Example mass and stiffness matrices
M = np.array([[2, 0], [0, 3]])  # 2x2 mass matrix
K = np.array([[6, -2], [-2, 4]])  # 2x2 stiffness matrix

# Displacement and force vectors
u = np.array([1, 2])
F = np.array([3, 4])

# Compute acceleration
M_inv = np.linalg.inv(M)
acceleration = M_inv @ (F - K @ u)

print("Acceleration Vector:", acceleration)

import numpy as np
from scipy.linalg import eigh

# Define the mass and stiffness matrices
M = np.array([[2, 0], [0, 3]])  # Mass matrix
K = np.array([[6, -2], [-2, 4]])  # Stiffness matrix

# Define the displacement and force vectors
u = np.array([1, 2])
F = np.array([3, 4])

# Compute acceleration using the equation of motion
M_inv = np.linalg.inv(M)
acceleration = M_inv @ (F - K @ u)
print("Acceleration Vector:", acceleration)

# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = eigh(K, M)

# Normalize eigenvectors
eigenvectors = eigenvectors / np.linalg.norm(eigenvectors, axis=0)

# Verify orthogonality
orthogonality_matrix = eigenvectors.T @ eigenvectors
print("Orthogonality Matrix:\n", orthogonality_matrix)

# Transform the force vector into modal coordinates
F_modal = np.linalg.inv(eigenvectors.T @ M @ eigenvectors) @ (eigenvectors.T @ F)

# Assume the solution in modal coordinates (for simplicity, use zero initial conditions)
q_modal = np.zeros_like(F_modal)

# Transform the solution back to physical coordinates
u_response = eigenvectors @ q_modal
print("Response Vector:", u_response)

# Compute modal coordinates
q = np.linalg.inv(eigenvectors) @ u
print("Modal Coordinates:", q)

import numpy as np
import matplotlib.pyplot as plt

# Define system parameters
m = 1.0  # Mass (kg)
k = 100.0  # Stiffness (N/m)
c = 0.0  # Damping coefficient (N·s/m) - for undamped system

# Define frequency range
omega = np.linspace(0.1, 10.0, 500)  # Frequency range (rad/s)

# Calculate frequency response for an undamped system
def frequency_response_undamped(omega, m, k):
    return 1 / np.sqrt((k/m - omega**2)**2)

# Calculate frequency response for a damped system
def frequency_response_damped(omega, m, k, c):
    return 1 / np.sqrt(((k/m - omega**2)**2 + (c*omega/m)**2))

# Calculate the response for undamped and damped cases
response_undamped = frequency_response_undamped(omega, m, k)
response_damped = frequency_response_damped(omega, m, k, c)

# Plot frequency response curves
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(omega, response_undamped, label='Undamped System')
plt.title('Frequency Response Curve (Undamped)')
plt.xlabel('Frequency (rad/s)')
plt.ylabel('Amplitude Response')
plt.grid(True)
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(omega, response_damped, label='Damped System', color='orange')
plt.title('Frequency Response Curve (Damped)')
plt.xlabel('Frequency (rad/s)')
plt.ylabel('Amplitude Response')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# System parameters
m = 1.0  # Mass (kg)
k = 100.0  # Spring constant (N/m)
m_e = 0.1  # Eccentric mass (kg)
e = 0.2  # Eccentric distance (m)
omega = 2 * np.pi * 1.0  # Angular velocity (rad/s)

# Define the system of differential equations
def system_dynamics(t, y):
    # y[0] = displacement, y[1] = velocity
    dydt = [y[1], (-k * y[0] + m_e * e * omega**2 * np.sin(omega * t)) / m]
    return dydt

# Initial conditions
y0 = [0.0, 0.0]  # Initial displacement and velocity

# Time span
t_span = (0, 10)  # Start and end time
t_eval = np.linspace(t_span[0], t_span[1], 500)  # Evaluation times

# Solve the differential equations
sol = solve_ivp(system_dynamics, t_span, y0, t_eval=t_eval, method='RK45')

# Plot the results
plt.figure(figsize=(10, 5))
plt.plot(sol.t, sol.y[0], label='Displacement')
plt.title('Response of a System with Eccentric Mass')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import scipy.linalg

# System parameters
E = 210e9  # Young's modulus (Pa)
I = 1e-4   # Moment of inertia (m^4)
L = 1.0    # Length of the beam (m)
n = 10     # Number of elements

# Discretize the beam into elements
dx = L / n
k_local = np.array([[12, 6*dx, -12, 6*dx],
                    [6*dx, 4*dx**2, -6*dx, 2*dx**2],
                    [-12, -6*dx, 12, -6*dx],
                    [6*dx, 2*dx**2, -6*dx, 4*dx**2]])

m_local = np.array([[156, 22*dx, 54, -13*dx],
                    [22*dx, 4*dx**2, 13*dx, -3*dx**2],
                    [54, 13*dx, 156, -22*dx],
                    [-13*dx, -3*dx**2, -22*dx, 4*dx**2]])

K_global = np.zeros((n+1, n+1))
M_global = np.zeros((n+1, n+1))

for i in range(n):
    K_global[i:i+2, i:i+2] += k_local[:2, :2]
    M_global[i:i+2, i:i+2] += m_local[:2, :2]

# Apply boundary conditions (fixed at both ends)
K_global = K_global[1:-1, 1:-1]
M_global = M_global[1:-1, 1:-1]

# Solve the generalized eigenvalue problem
eigvals, eigvecs = scipy.linalg.eigh(K_global, M_global)

# Extract natural frequencies
natural_frequencies = np.sqrt(eigvals) / (2 * np.pi)

print("Natural Frequencies (Hz):", natural_frequencies)

import numpy as np
from scipy.linalg import eigh

# Define parameters
m1, m2 = 1.0, 1.0  # Masses
k1, k2 = 10.0, 10.0  # Stiffnesses
alpha, beta = 0.1, 0.2  # Damping coefficients

# Define matrices
M = np.array([[m1, 0],
              [0, m2]])

K = np.array([[k1 + k2, -k2],
              [-k2, k2]])

C = alpha * M + beta * K

# Eigenvalue problem
A = np.linalg.inv(M) @ K
B = np.linalg.inv(M) @ C

# Solve generalized eigenvalue problem
eigenvalues, eigenvectors = eigh(A, B)

# Convert to natural frequencies
natural_frequencies = np.sqrt(eigenvalues)

print("Natural Frequencies:", natural_frequencies)
print("Mode Shapes:\n", eigenvectors)

import math

# Constants (modify these as needed)
k = 100.0         # Spring constant (N/m)
m = 5.0           # Mass (kg)
g = 9.81          # Acceleration due to gravity (m/s^2)
l = 2.0           # Length (m)
T = 50.0          # Tension (N)
rho = 0.01        # Linear density (kg/m)
E = 200e9         # Young's modulus (Pa)
I = 1e-6          # Area moment of inertia (m^4)
A = 0.01          # Cross-sectional area (m^2)
n = 1             # Mode number
G = 80e9          # Shear modulus (Pa)
J = 1e-6          # Polar moment of inertia (m^4)
lambda_value = math.pi  # Example value for λ (π, 2π, etc.)

def calc_frequency():
    print("Calculating frequencies using predefined constants...")

    # Formula 1: ω = sqrt(k / m)
    omega_1 = math.sqrt(k / m)
    print(f"ω (Formula 1) = {omega_1:.4f} rad/s")

    # Formula 2: ω = sqrt(g / l)
    omega_2 = math.sqrt(g / l)
    print(f"ω (Formula 2) = {omega_2:.4f} rad/s")

    # Formula 3: ω_n = (nπ / l) * sqrt(T / ρ)
    omega_n = (n * math.pi / l) * math.sqrt(T / rho)
    print(f"ω_n (Formula 3) = {omega_n:.4f} rad/s")

    # Formula 4: f_n = (1 / 2π) * (nπ / l)^2 * sqrt(EI / ρA)
    f_n_4 = (1 / (2 * math.pi)) * ((n * math.pi / l) ** 2) * math.sqrt((E * I) / (rho * A))
    print(f"f_n (Formula 4) = {f_n_4:.4f} Hz")

    # Formula 5: f_n = (1 / 2π) * (nπ / l)^2 * sqrt(E / ρ) (Rectangular Cross-Section)
    f_n_5 = (1 / (2 * math.pi)) * ((n * math.pi / l) ** 2) * math.sqrt(E / rho)
    print(f"f_n (Formula 5, Rectangular Cross-Section) = {f_n_5:.4f} Hz")

    # Formula 6: f_n = (λ / 2πl) * sqrt(GJ / ρA)
    f_n_6 = (lambda_value / (2 * math.pi * l)) * math.sqrt((G * J) / (rho * A))
    print(f"f_n (Formula 6) = {f_n_6:.4f} Hz")

    # Formula 7: f_n = (λ / 2πl) * sqrt(E / ρ)
    f_n_7 = (lambda_value / (2 * math.pi * l)) * math.sqrt(E / rho)
    print(f"f_n (Formula 7) = {f_n_7:.4f} Hz")

calc_frequency()