# Untitled

unknown
python
a month ago
5.1 kB
2
Indexable
Never
```# -*- coding: utf-8 -*-
"""
Created on Thu Aug 15 18:11:39 2024

@author: Danny
"""

import numpy as np
from numpy.polynomial import Polynomial
from numpy.polynomial.polynomial import polyval
import matplotlib.pyplot as plt

def zero_cross(ar):
H1t=np.sign(ar)
H1s=np.abs(H1t[:-1]-H1t[1:])
return np.nonzero(H1s)[0]

def Power_Spectrum(t,x):
# t is the time array x is the signal array
T=t[-1] # T is the last element of the time array
dt=t[1]-t[0]
df=1/T # 1/T is the freq resolution
X=np.fft.fft(x) # calculate the fast fourier transform of the signal
P=np.abs(X) # find the power spectrum
f=np.arange(0,1/dt,df )
halfN=t.size//2 # only the first half of the power spectrum is of interest
return f[:halfN],P[:halfN]
#input
f_in=float(input('Enter frequncy (Hz): '))

#(1)
# Function to simulate an ideal spring
def ideal(x0, v0, k):
w0 = np.sqrt(k/m)
A, f = [], []

for x1 in x0:
x = np.zeros_like(t)
v = np.zeros_like(t)
x[0], v[0] = x1, v0

for i in range(1, len(t)):
v[i] = v[i-1] - (w0**2) * x[i-1] * dt1
x[i] = x[i-1] + v[i] * dt1

plt.plot(t, x, label=f"x0 = {x1} m")
A.append(np.max(x))

zero_crossings = zero_cross(x)
if len(zero_crossings) > 1:
f.append(1 / (t[zero_crossings[1]] - t[zero_crossings[0]]))

plt.title("Ideal Spring Motion")
plt.xlabel("Time (s)")
plt.ylabel("Position (m)")
plt.legend()
plt.show()

plt.plot(A, f, '-o')
plt.title("Amplitude vs Frequency (Ideal Spring)")
plt.xlabel("Amplitude (m)")
plt.ylabel("Frequency (Hz)")
plt.show()

f_spectrum, P_spectrum = Power_Spectrum(t, x)
plt.plot(f_spectrum, P_spectrum)
plt.title('Power Spectrum (Ideal Spring)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.grid(True)
plt.show()
#(2)
# Function to simulate a non-ideal spring
def not_ideal(x0, v0, k):
global poly
w0 = np.sqrt(k/m)
A, f = [], []

for x1 in x0:
x = np.zeros_like(t)
v = np.zeros_like(t)
x[0], v[0] = x1, v0

for i in range(1, len(t)):
v[i] = v[i-1] - (w0**2) * np.power(x[i-1], 3) * dt1
x[i] = x[i-1] + v[i] * dt1

plt.plot(t, x, label=f"x0 = {x1} m")
A.append(np.max(x))

zero_crossings = zero_cross(x)
if len(zero_crossings) > 1:
f.append(1 / (t[zero_crossings[2]] - t[zero_crossings[0]]))

plt.title("Non-Ideal Spring Motion")
plt.xlabel("Time (s)")
plt.ylabel("Position (m)")
plt.legend()
plt.show()

poly = Polynomial.fit(f, A, deg=3)
plt.plot(A, f, '-o')
plt.title("Amplitude vs Frequency (Non-Ideal Spring)")
plt.xlabel("Amplitude (m)")
plt.ylabel("Frequency (Hz)")
plt.show()

f_spectrum, P_spectrum = Power_Spectrum(t, x)
plt.plot(f_spectrum, P_spectrum)
plt.title('Power Spectrum (Non-Ideal Spring)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.grid(True)
plt.show()

#parameters
dt1 = 1e-2
t = np.arange(0, 10, dt1)
K_ideal = 20
K_not_ideal = 5
m = 1
x0 = [1, 1.2, 1.4, 1.6, 1.8, 2]

# Run ideal spring and not_ideal
ideal(x0, 0, K_ideal)
not_ideal(x0, 0, K_not_ideal)

# Calculate output for given input frequency
f_in = float(input('Enter frequency (Hz): '))
poly_coef = poly.convert().coef
A_of_f = polyval(f_in, poly_coef)
print(f'{A_of_f:.4g}')

#(3)
#parameters
x0 = float(input("Initial position: "))
v0 = float(input("Initial velocity: "))
t_input = float(input("Time: "))
dt = 1e-3
t_array = np.arange(0, 10, dt)
m = 2
k = 40
b = 5

gamma = b / (2 * m)
w0 = np.sqrt(k / m)
w = np.sqrt(w0**2 - gamma**2)

# Initialize arrays for position and velocity
v = np.zeros_like(t_array)
x = np.zeros_like(t_array)

# Initial conditions
v[0] = v0
x[0] = x0

# Numerical integration using Euler's method
for i in range(len(t_array) - 1):
v[i + 1] = v[i] + (-2 * gamma * v[i] - w0**2 * x[i]) * dt
x[i + 1] = x[i] + v[i + 1] * dt

# Find velocity at the input time
t_index = np.argmin(np.abs(t_array - t_input))
print(f"Velocity at t={t_input}s: {v[t_index]:.5f} m/s")

# Analytical solution
alpha = -gamma
A = x0 / np.cos(np.arctan((-x0 * gamma - v0) / (x0 * w0)))
phi = np.arctan((-x0 * gamma - v0) / (x0 * w0))
x_theory = A * np.exp(alpha * t_array) * np.cos(w0 * t_array + phi)

# Plot position as a function of time
plt.figure(figsize=(10, 6))
plt.plot(t_array, x, label='Numerical')
plt.plot(t_array, x_theory, '--', label='Analytical')
plt.xlabel('Time (s)')
plt.ylabel('Position (m)')
plt.title('Position vs Time')
plt.legend()
plt.grid(True)
plt.show()

# Plot Power Spectrum
f, P = Power_Spectrum(t_array, x)
plt.plot(f, P)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.title('Power Spectrum')
plt.grid(True)
plt.show()```