NumPy Numerical Analysis Skill
Master NumPy for efficient numerical computations, matrix operations, FFT analysis, and linear algebra in marine and offshore engineering applications.
When to Use This Skill
Use NumPy numerical analysis when you need:
- Matrix operations - 6DOF equations of motion, mass matrices, stiffness matrices
- FFT analysis - Frequency domain analysis, spectral density, response spectra
- Linear algebra - Solve linear systems, eigenvalue analysis, matrix decomposition
- Array operations - Efficient computations on large datasets
- Numerical integration - Time-stepping, ODE solvers
- Signal processing - Filtering, windowing, convolution
Avoid when:
- Symbolic mathematics needed (use SymPy)
- Sparse matrices dominate (use SciPy sparse)
- GPU acceleration required (use CuPy or JAX)
- Distributed computing needed (use Dask)
Core Capabilities
1. Array Creation and Operations
Array Creation:
import numpy as np
# Create arrays
zeros = np.zeros((3, 3))
ones = np.ones((3, 3))
identity = np.eye(3)
arange = np.arange(0, 10, 0.1) # 0 to 10 with step 0.1
linspace = np.linspace(0, 10, 100) # 100 points from 0 to 10
# From list
arr = np.array([1, 2, 3, 4, 5])
# Multi-dimensional
matrix = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
# Random arrays
random_uniform = np.random.rand(3, 3) # Uniform [0, 1)
random_normal = np.random.randn(3, 3) # Standard normal
random_int = np.random.randint(0, 100, size=(3, 3))Array Operations:
# Element-wise operations
a = np.array([1, 2, 3, 4, 5])
b = np.array([10, 20, 30, 40, 50])
c = a + b # [11, 22, 33, 44, 55]
d = a * b # [10, 40, 90, 160, 250]
e = a ** 2 # [1, 4, 9, 16, 25]
# Mathematical functions
sin_a = np.sin(a)
cos_a = np.cos(a)
exp_a = np.exp(a)
log_a = np.log(a)
sqrt_a = np.sqrt(a)
# Statistical operations
mean = np.mean(a)
std = np.std(a)
var = np.var(a)
min_val = np.min(a)
max_val = np.max(a)2. Matrix Operations
Matrix Multiplication:
def compute_force_response(
mass_matrix: np.ndarray,
stiffness_matrix: np.ndarray,
force_vector: np.ndarray
) -> np.ndarray:
"""
Compute structural response: F = K * x
Solve for displacement: x = K^-1 * F
Args:
mass_matrix: Mass matrix [M]
stiffness_matrix: Stiffness matrix [K]
force_vector: Applied force vector {F}
Returns:
Displacement vector {x}
"""
# Static response (ignoring mass for now)
displacement = np.linalg.solve(stiffness_matrix, force_vector)
return displacement
# Example: 3DOF spring-mass system
K = np.array([
[200, -100, 0],
[-100, 200, -100],
[0, -100, 100]
]) # Stiffness matrix (N/m)
F = np.array([1000, 0, 0]) # Force at first node (N)
x = compute_force_response(None, K, F)
print(f"Displacements: {x} m")Matrix Properties:
def analyze_matrix_properties(matrix: np.ndarray) -> dict:
"""
Analyze matrix properties for structural analysis.
Args:
matrix: Input matrix (mass or stiffness)
Returns:
Dictionary with matrix properties
"""
properties = {}
# Determinant
properties['determinant'] = np.linalg.det(matrix)
# Condition number (numerical stability indicator)
properties['condition_number'] = np.linalg.cond(matrix)
# Rank
properties['rank'] = np.linalg.matrix_rank(matrix)
# Eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(matrix)
properties['eigenvalues'] = eigenvalues
properties['eigenvectors'] = eigenvectors
# Is symmetric?
properties['is_symmetric'] = np.allclose(matrix, matrix.T)
# Is positive definite? (all eigenvalues > 0)
properties['is_positive_definite'] = np.all(eigenvalues > 0)
return properties
# Example: Check stiffness matrix properties
K = np.array([
[200, -100, 0],
[-100, 200, -100],
[0, -100, 100]
])
props = analyze_matrix_properties(K)
print(f"Determinant: {props['determinant']:.2f}")
print(f"Condition number: {props['condition_number']:.2f}")
print(f"Eigenvalues: {props['eigenvalues']}")
print(f"Positive definite: {props['is_positive_definite']}")3. 6DOF Equations of Motion
6DOF Dynamics:
def solve_6dof_equation_of_motion(
mass_matrix: np.ndarray,
damping_matrix: np.ndarray,
stiffness_matrix: np.ndarray,
force_vector: np.ndarray,
displacement: np.ndarray,
velocity: np.ndarray,
dt: float
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Solve 6DOF equation of motion using Newmark-Beta method.
[M]{ẍ} + [C]{ẋ} + [K]{x} = {F}
Args:
mass_matrix: 6x6 mass matrix [M]
damping_matrix: 6x6 damping matrix [C]
stiffness_matrix: 6x6 stiffness matrix [K]
force_vector: 6x1 force vector {F}
displacement: Current displacement {x_n}
velocity: Current velocity {ẋ_n}
dt: Time step
Returns:
(acceleration, velocity, displacement) at next time step
"""
# Newmark-Beta parameters
beta = 0.25 # Constant average acceleration
gamma = 0.5
# Effective stiffness
K_eff = (
mass_matrix / (beta * dt**2) +
damping_matrix * gamma / (beta * dt) +
stiffness_matrix
)
# Effective force at next time step
F_eff = (
force_vector +
mass_matrix @ (
displacement / (beta * dt**2) +
velocity / (beta * dt)
) +
damping_matrix @ (
displacement * gamma / (beta * dt) -
velocity * (1 - gamma / beta)
)
)
# Solve for displacement at next time step
displacement_next = np.linalg.solve(K_eff, F_eff)
# Calculate velocity at next time step
velocity_next = (
gamma / (beta * dt) * (displacement_next - displacement) +
(1 - gamma / beta) * velocity
)
# Calculate acceleration at next time step
acceleration_next = (
(displacement_next - displacement) / (beta * dt**2) -
velocity / (beta * dt)
)
return acceleration_next, velocity_next, displacement_next
# Example: Single time step for floating vessel
M = np.diag([100000, 100000, 100000, 5e6, 5e6, 5e6]) # Mass matrix
C = np.diag([50000, 50000, 50000, 2e5, 2e5, 2e5]) # Damping
K = np.diag([1000, 1000, 2000, 5e4, 5e4, 1e4]) # Restoring
F = np.array([1e6, 0, 0, 0, 0, 0]) # Surge force 1 MN
x = np.zeros(6) # Initial displacement
v = np.zeros(6) # Initial velocity
a, v_new, x_new = solve_6dof_equation_of_motion(M, C, K, F, x, v, dt=0.01)
print(f"Acceleration: {a}")
print(f"Velocity: {v_new}")
print(f"Displacement: {x_new}")4. FFT and Frequency Analysis
FFT for Spectral Analysis:
def compute_fft_spectrum(
time_series: np.ndarray,
dt: float,
window: str = 'hann'
) -> tuple[np.ndarray, np.ndarray]:
"""
Compute FFT spectrum of time series.
Args:
time_series: Time series data
dt: Time step
window: Window function ('hann', 'hamming', 'blackman')
Returns:
(frequencies, amplitude_spectrum)
"""
n = len(time_series)
# Apply window
if window == 'hann':
windowed = time_series * np.hanning(n)
elif window == 'hamming':
windowed = time_series * np.hamming(n)
elif window == 'blackman':
windowed = time_series * np.blackman(n)
else:
windowed = time_series
# Compute FFT
fft_result = np.fft.fft(windowed)
# Compute frequencies
frequencies = np.fft.fftfreq(n, d=dt)
# Amplitude spectrum (single-sided)
amplitude = np.abs(fft_result)[:n//2] * 2 / n
frequencies_positive = frequencies[:n//2]
return frequencies_positive, amplitude
# Example: Analyze wave elevation time series
import numpy as np
# Generate sample wave elevation (3 components)
t = np.linspace(0, 100, 10000) # 100 seconds, 10000 points
dt = t[1] - t[0]
# Wave components: 6s, 8s, 10s periods
wave = (
2.0 * np.sin(2*np.pi*t / 6) +
1.5 * np.sin(2*np.pi*t / 8) +
1.0 * np.sin(2*np.pi*t / 10)
)
# Add noise
wave += 0.2 * np.random.randn(len(t))
# Compute spectrum
freq, amplitude = compute_fft_spectrum(wave, dt, window='hann')
# Find peaks
peak_indices = np.argsort(amplitude)[-3:] # Top 3 peaks
peak_frequencies = freq[peak_indices]
peak_periods = 1 / peak_frequencies
print("Detected wave periods:")
for period in sorted(peak_periods, reverse=True):
print(f" T = {period:.2f} s")Power Spectral Density:
def compute_power_spectral_density(
time_series: np.ndarray,
dt: float,
nfft: int = None
) -> tuple[np.ndarray, np.ndarray]:
"""
Compute power spectral density using Welch's method.
Args:
time_series: Time series data
dt: Time step
nfft: FFT length (None = length of time series)
Returns:
(frequencies, PSD)
"""
from scipy import signal
# Compute PSD using Welch's method
frequencies, psd = signal.welch(
time_series,
fs=1/dt,
nperseg=nfft or len(time_series)//8,
window='hann'
)
return frequencies, psd
# Example: Wave spectral analysis
t = np.linspace(0, 3600, 36000) # 1 hour, 10 Hz sampling
dt = t[1] - t[0]
# JONSWAP spectrum simulation (simplified)
wave_elevation = np.random.randn(len(t)) * 2.0 # Simplified
freq, psd = compute_power_spectral_density(wave_elevation, dt)
# Calculate Hs from PSD
m0 = np.trapz(psd, freq) # Zero-order moment
Hs = 4 * np.sqrt(m0)
print(f"Significant wave height: {Hs:.2f} m")5. Linear Algebra Operations
Eigenvalue Analysis:
def natural_frequency_analysis(
mass_matrix: np.ndarray,
stiffness_matrix: np.ndarray
) -> dict:
"""
Perform eigenvalue analysis to find natural frequencies and mode shapes.
[K]{ϕ} = ω²[M]{ϕ}
Args:
mass_matrix: Mass matrix [M]
stiffness_matrix: Stiffness matrix [K]
Returns:
Dictionary with natural frequencies and mode shapes
"""
# Solve generalized eigenvalue problem
eigenvalues, eigenvectors = np.linalg.eig(
np.linalg.solve(mass_matrix, stiffness_matrix)
)
# Natural frequencies (rad/s)
natural_frequencies_rad = np.sqrt(eigenvalues)
# Natural frequencies (Hz)
natural_frequencies_hz = natural_frequencies_rad / (2 * np.pi)
# Sort by frequency
sort_indices = np.argsort(natural_frequencies_hz)
natural_frequencies_hz = natural_frequencies_hz[sort_indices]
eigenvectors = eigenvectors[:, sort_indices]
# Periods
periods = 1 / natural_frequencies_hz
return {
'frequencies_hz': natural_frequencies_hz,
'frequencies_rad_s': natural_frequencies_rad[sort_indices],
'periods_s': periods,
'mode_shapes': eigenvectors
}
# Example: FPSO natural frequencies
# 6DOF system
M = np.diag([150000, 150000, 150000, 1e7, 1e7, 5e6]) # Mass matrix (tonnes, tonne-m²)
K = np.diag([500, 500, 3000, 5e5, 5e5, 1e5]) # Stiffness (kN/m, kN-m/rad)
results = natural_frequency_analysis(M, K)
print("Natural Frequencies:")
for i, (freq, period) in enumerate(zip(results['frequencies_hz'], results['periods_s'])):
dof_names = ['Surge', 'Sway', 'Heave', 'Roll', 'Pitch', 'Yaw']
print(f" Mode {i+1} ({dof_names[i]}): f = {freq:.4f} Hz, T = {period:.2f} s")Matrix Decomposition:
def lu_decomposition_solve(A: np.ndarray, b: np.ndarray) -> np.ndarray:
"""
Solve linear system using LU decomposition.
Args:
A: Coefficient matrix
b: Right-hand side vector
Returns:
Solution vector x
"""
from scipy.linalg import lu
# LU decomposition
P, L, U = lu(A)
# Solve Ly = Pb
y = np.linalg.solve(L, P @ b)
# Solve Ux = y
x = np.linalg.solve(U, y)
return x
# Example
A = np.array([[4, -1, 0],
[-1, 4, -1],
[0, -1, 3]])
b = np.array([15, 10, 10])
x = lu_decomposition_solve(A, b)
print(f"Solution: {x}")6. Numerical Integration
Time-Stepping Integration:
def runge_kutta_4th_order(
derivative_func,
y0: np.ndarray,
t: np.ndarray
) -> np.ndarray:
"""
4th-order Runge-Kutta integration.
Args:
derivative_func: Function dy/dt = f(t, y)
y0: Initial conditions
t: Time array
Returns:
Solution array
"""
n = len(t)
y = np.zeros((n, len(y0)))
y[0] = y0
for i in range(n - 1):
dt = t[i+1] - t[i]
k1 = derivative_func(t[i], y[i])
k2 = derivative_func(t[i] + dt/2, y[i] + k1*dt/2)
k3 = derivative_func(t[i] + dt/2, y[i] + k2*dt/2)
k4 = derivative_func(t[i] + dt, y[i] + k3*dt)
y[i+1] = y[i] + (dt/6) * (k1 + 2*k2 + 2*k3 + k4)
return y
# Example: Simple harmonic oscillator
# m*x'' + k*x = 0
# y = [x, v], dy/dt = [v, -k/m*x]
def oscillator_derivatives(t, y):
"""Simple harmonic oscillator."""
m = 1.0 # Mass
k = 4.0 # Stiffness (omega = 2 rad/s, T = pi s)
x, v = y
dxdt = v
dvdt = -k/m * x
return np.array([dxdt, dvdt])
# Initial conditions
y0 = np.array([1.0, 0.0]) # x = 1, v = 0
# Time array
t = np.linspace(0, 10, 1000)
# Solve
solution = runge_kutta_4th_order(oscillator_derivatives, y0, t)
print(f"Final position: {solution[-1, 0]:.4f}")
print(f"Final velocity: {solution[-1, 1]:.4f}")Trapezoidal Integration:
def integrate_spectrum(
frequencies: np.ndarray,
spectral_density: np.ndarray
) -> float:
"""
Integrate spectral density to get variance.
m_0 = ∫ S(f) df
Args:
frequencies: Frequency array
spectral_density: Spectral density array
Returns:
Integral (variance)
"""
variance = np.trapz(spectral_density, frequencies)
return variance
# Example: Calculate significant wave height from spectrum
freq = np.linspace(0.05, 0.5, 100) # Hz
S = 10 * freq**(-5) # Simplified spectrum
m0 = integrate_spectrum(freq, S)
Hs = 4 * np.sqrt(m0)
print(f"Significant wave height: {Hs:.2f} m")Complete Examples
Example 1: 6DOF Time-Domain Simulation
import numpy as np
import plotly.graph_objects as go
def simulate_6dof_vessel_motion(
mass_matrix: np.ndarray,
damping_matrix: np.ndarray,
stiffness_matrix: np.ndarray,
force_time_series: np.ndarray,
time: np.ndarray
) -> dict:
"""
Complete 6DOF time-domain simulation of vessel motion.
Args:
mass_matrix: 6x6 mass matrix
damping_matrix: 6x6 damping matrix
stiffness_matrix: 6x6 stiffness matrix
force_time_series: Force time series (n_steps x 6)
time: Time array
Returns:
Dictionary with motion time series
"""
n_steps = len(time)
dt = time[1] - time[0]
# Initialize arrays
displacement = np.zeros((n_steps, 6))
velocity = np.zeros((n_steps, 6))
acceleration = np.zeros((n_steps, 6))
# Initial conditions (at rest)
displacement[0] = np.zeros(6)
velocity[0] = np.zeros(6)
# Time-stepping loop
for i in range(n_steps - 1):
a, v, d = solve_6dof_equation_of_motion(
mass_matrix,
damping_matrix,
stiffness_matrix,
force_time_series[i],
displacement[i],
velocity[i],
dt
)
acceleration[i+1] = a
velocity[i+1] = v
displacement[i+1] = d
return {
'time': time,
'displacement': displacement,
'velocity': velocity,
'acceleration': acceleration
}
# Define FPSO properties
M = np.diag([150000, 150000, 150000, 1e7, 1e7, 5e6]) # Mass
C = np.diag([50000, 50000, 100000, 5e5, 5e5, 2e5]) # Damping
K = np.diag([500, 500, 3000, 5e4, 5e4, 1e4]) # Stiffness
# Simulate wave forces (simplified)
time = np.linspace(0, 100, 10000)
dt = time[1] - time[0]
# Wave force in surge (1 MN amplitude, 10s period)
F = np.zeros((len(time), 6))
F[:, 0] = 1e6 * np.sin(2*np.pi*time / 10)
# Simulate
results = simulate_6dof_vessel_motion(M, C, K, F, time)
# Plot results
fig = go.Figure()
dof_names = ['Surge', 'Sway', 'Heave', 'Roll', 'Pitch', 'Yaw']
for i in range(3): # Plot first 3 DOFs
fig.add_trace(go.Scatter(
x=results['time'],
y=results['displacement'][:, i],
name=dof_names[i],
mode='lines'
))
fig.update_layout(
title='Vessel Motion - Translation DOFs',
xaxis_title='Time (s)',
yaxis_title='Displacement (m)',
hovermode='x unified'
)
fig.write_html('reports/vessel_motion_6dof.html')
# Print statistics
print("Motion Statistics:")
for i, name in enumerate(dof_names):
print(f"{name}:")
print(f" Max: {np.max(np.abs(results['displacement'][:, i])):.3f}")
print(f" Std: {np.std(results['displacement'][:, i]):.3f}")Example 2: RAO Calculation from FFT
def calculate_rao_from_time_series(
wave_elevation: np.ndarray,
vessel_response: np.ndarray,
dt: float
) -> tuple[np.ndarray, np.ndarray]:
"""
Calculate Response Amplitude Operator (RAO) from time series.
RAO(ω) = |Response(ω)| / |Wave(ω)|
Args:
wave_elevation: Wave elevation time series
vessel_response: Vessel response time series
dt: Time step
Returns:
(frequencies, RAO)
"""
# FFT of wave and response
freq_wave, amplitude_wave = compute_fft_spectrum(wave_elevation, dt)
freq_response, amplitude_response = compute_fft_spectrum(vessel_response, dt)
# Calculate RAO
# Avoid division by zero
rao = np.where(
amplitude_wave > 1e-6,
amplitude_response / amplitude_wave,
0
)
return freq_wave, rao
# Example: Calculate heave RAO
time = np.linspace(0, 100, 10000)
dt = time[1] - time[0]
# Regular wave (Hs = 2m, T = 8s)
wave = 1.0 * np.sin(2*np.pi*time / 8)
# Vessel heave response (with RAO ~1.5)
heave = 1.5 * np.sin(2*np.pi*time / 8 - np.pi/6) # 30° phase lag
# Calculate RAO
freq, rao = calculate_rao_from_time_series(wave, heave, dt)
# Plot
fig = go.Figure()
fig.add_trace(go.Scatter(
x=freq,
y=rao,
name='Heave RAO',
mode='lines'
))
fig.update_layout(
title='Heave Response Amplitude Operator',
xaxis_title='Frequency (Hz)',
yaxis_title='RAO (m/m)',
xaxis_range=[0, 0.5]
)
fig.write_html('reports/heave_rao.html')
# Find peak
peak_idx = np.argmax(rao[freq < 0.5])
peak_freq = freq[peak_idx]
peak_period = 1 / peak_freq
print(f"Peak RAO: {rao[peak_idx]:.2f} at T = {peak_period:.2f} s")Example 3: Mooring Stiffness Matrix
def calculate_mooring_stiffness_matrix(
num_lines: int,
pretension: float,
fairlead_radius: float,
fairlead_depth: float,
line_azimuth: np.ndarray,
weight_per_length: float
) -> np.ndarray:
"""
Calculate mooring system stiffness matrix.
Args:
num_lines: Number of mooring lines
pretension: Pretension per line (kN)
fairlead_radius: Horizontal distance from center to fairlead (m)
fairlead_depth: Depth of fairlead below waterline (m)
line_azimuth: Azimuth angle of each line (degrees)
weight_per_length: Line weight per unit length (kN/m)
Returns:
6x6 mooring stiffness matrix
"""
K_mooring = np.zeros((6, 6))
# Individual line stiffness
# Simplified: K_line = w * T / d
k_line = weight_per_length * pretension / fairlead_depth
for i in range(num_lines):
theta = np.radians(line_azimuth[i])
# Direction cosines
cx = np.cos(theta)
cy = np.sin(theta)
# Surge-surge
K_mooring[0, 0] += k_line * cx**2
# Sway-sway
K_mooring[1, 1] += k_line * cy**2
# Surge-sway coupling
K_mooring[0, 1] += k_line * cx * cy
K_mooring[1, 0] += k_line * cx * cy
# Yaw-yaw (moment = force * lever arm)
K_mooring[5, 5] += k_line * fairlead_radius**2
# Surge-yaw coupling
K_mooring[0, 5] += k_line * fairlead_radius * cy
K_mooring[5, 0] += k_line * fairlead_radius * cy
# Sway-yaw coupling
K_mooring[1, 5] -= k_line * fairlead_radius * cx
K_mooring[5, 1] -= k_line * fairlead_radius * cx
return K_mooring
# Example: 12-line spread mooring
num_lines = 12
azimuths = np.array([0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330])
K = calculate_mooring_stiffness_matrix(
num_lines=12,
pretension=2000, # kN
fairlead_radius=100, # m
fairlead_depth=25, # m
line_azimuth=azimuths,
weight_per_length=1.2 # kN/m
)
print("Mooring Stiffness Matrix:")
print(K)
# Check symmetry
print(f"\nIs symmetric: {np.allclose(K, K.T)}")Example 4: Statistical Analysis of Extremes
def extreme_value_statistics(
data: np.ndarray,
method: str = '3hr_max'
) -> dict:
"""
Perform extreme value statistical analysis.
Args:
data: Time series data
method: '3hr_max' or 'annual_max'
Returns:
Statistical parameters
"""
if method == '3hr_max':
# Extract 3-hour maxima
chunk_size = int(3 * 3600 / 0.1) # Assuming 0.1s time step
n_chunks = len(data) // chunk_size
maxima = np.array([
np.max(data[i*chunk_size:(i+1)*chunk_size])
for i in range(n_chunks)
])
elif method == 'annual_max':
maxima = data # Assume already annual maxima
# Fit Gumbel distribution
# μ = mean, σ = std
mu = np.mean(maxima)
sigma = np.std(maxima)
# Gumbel parameters
beta = sigma * np.sqrt(6) / np.pi
mu_gumbel = mu - 0.5772 * beta
# Extreme values for different return periods
return_periods = np.array([1, 10, 100, 10000]) # years
extreme_values = mu_gumbel - beta * np.log(-np.log(1 - 1/return_periods))
return {
'mean': mu,
'std': sigma,
'gumbel_location': mu_gumbel,
'gumbel_scale': beta,
'return_periods': return_periods,
'extreme_values': extreme_values
}
# Example
# Generate sample tension data
time = np.linspace(0, 86400, 864000) # 24 hours, 0.1s timestep
tension = 2000 + 500 * np.random.rayleigh(scale=1.0, size=len(time))
stats = extreme_value_statistics(tension, method='3hr_max')
print("Extreme Value Statistics:")
print(f"Mean: {stats['mean']:.1f} kN")
print(f"Std: {stats['std']:.1f} kN")
print("\nExtreme Values:")
for T, val in zip(stats['return_periods'], stats['extreme_values']):
print(f" {T:>5} year: {val:.1f} kN")Example 5: Convolution for Impulse Response
def convolve_impulse_response(
impulse_response: np.ndarray,
force_time_series: np.ndarray,
dt: float
) -> np.ndarray:
"""
Convolve impulse response with force time series.
Response(t) = ∫ h(τ) * F(t-τ) dτ
Args:
impulse_response: Impulse response function
force_time_series: Force time series
dt: Time step
Returns:
Response time series
"""
# Convolution
response = np.convolve(impulse_response, force_time_series, mode='same') * dt
return response
# Example: System response to random force
t = np.linspace(0, 10, 1000)
dt = t[1] - t[0]
# Impulse response (damped oscillator)
omega = 2 * np.pi # 1 Hz natural frequency
zeta = 0.1 # 10% damping
h = np.exp(-zeta * omega * t) * np.sin(omega * np.sqrt(1 - zeta**2) * t)
# Random force
F = np.random.randn(len(t))
# Compute response
response = convolve_impulse_response(h, F, dt)
print(f"Max response: {np.max(np.abs(response)):.3f}")
print(f"Std response: {np.std(response):.3f}")Best Practices
1. Use Vectorization
# ❌ Bad: Loop
result = np.zeros(len(x))
for i in range(len(x)):
result[i] = x[i]**2 + y[i]**2
# ✅ Good: Vectorized
result = x**2 + y**22. Avoid Unnecessary Copies
# ❌ Bad: Creates copies
a = np.array([1, 2, 3])
b = a
b[0] = 10 # Modifies original
# ✅ Good: Explicit copy when needed
a = np.array([1, 2, 3])
b = a.copy()
b[0] = 10 # Original unchanged3. Use In-Place Operations
# ❌ Bad: Creates new array
a = a + 1
# ✅ Good: In-place
a += 14. Choose Appropriate Data Types
# Use float32 for large arrays when precision allows
large_array = np.zeros((10000, 10000), dtype=np.float32) # 400 MB instead of 800 MBResources
- NumPy Documentation: https://numpy.org/doc/
- NumPy for MATLAB Users: https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
- Linear Algebra: https://numpy.org/doc/stable/reference/routines.linalg.html
- FFT Module: https://numpy.org/doc/stable/reference/routines.fft.html
- SciPy (extends NumPy): https://scipy.org/
Use this skill for all numerical computations in DigitalModel!