numpy-numerical-analysis

What this skill does

# 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:**
```python
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:**
```python
# 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:**
```python
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:**
```python
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:**
```python
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:**
```python
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] -