python-scientific-computing

What this skill does

# Python Scientific Computing Skill

Master Python for engineering analysis, numerical simulations, and scientific workflows using industry-standard libraries.

## When to Use This Skill

Use Python scientific computing when you need:
- **Numerical analysis** - Solving equations, optimization, integration
- **Engineering calculations** - Stress, strain, dynamics, thermodynamics
- **Matrix operations** - Linear algebra, eigenvalue problems
- **Symbolic mathematics** - Analytical solutions, equation manipulation
- **Data analysis** - Statistical analysis, curve fitting
- **Simulations** - Physical systems, finite element preprocessing

**Avoid when:**
- Real-time performance critical (use C++/Fortran)
- Simple calculations (use calculator or Excel)
- No numerical computation needed

## Core Capabilities

### 1. NumPy - Numerical Arrays and Linear Algebra

**Array Operations:**
```python
import numpy as np

# Create arrays
array_1d = np.array([1, 2, 3, 4, 5])
array_2d = np.array([[1, 2, 3], [4, 5, 6]])

# Special arrays
zeros = np.zeros((3, 3))
ones = np.ones((2, 4))
identity = np.eye(3)
linspace = np.linspace(0, 10, 100)  # 100 points from 0 to 10

# Array operations (vectorized - fast!)
x = np.linspace(0, 2*np.pi, 1000)
y = np.sin(x) * np.exp(-x/10)
```

**Linear Algebra:**
```python
# Matrix operations
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Matrix multiplication
C = A @ B  # or np.dot(A, B)

# Inverse
A_inv = np.linalg.inv(A)

# Eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)

# Solve linear system Ax = b
b = np.array([1, 2])
x = np.linalg.solve(A, b)

# Determinant
det_A = np.linalg.det(A)
```

### 2. SciPy - Scientific Computing

**Optimization:**
```python
from scipy import optimize

# Minimize function
def rosenbrock(x):
    return (1 - x[0])**2 + 100*(x[1] - x[0]**2)**2

result = optimize.minimize(rosenbrock, x0=[0, 0], method='BFGS')
print(f"Minimum at: {result.x}")

# Root finding
def equations(vars):
    x, y = vars
    eq1 = x**2 + y**2 - 4
    eq2 = x - y - 1
    return [eq1, eq2]

solution = optimize.fsolve(equations, [1, 1])
```

**Integration:**
```python
from scipy import integrate

# Numerical integration
def integrand(x):
    return x**2

result, error = integrate.quad(integrand, 0, 1)  # Integrate from 0 to 1
print(f"Result: {result}, Error: {error}")

# ODE solver
def ode_system(t, y):
    # dy/dt = -2y
    return -2 * y

solution = integrate.solve_ivp(
    ode_system,
    t_span=[0, 10],
    y0=[1],
    t_eval=np.linspace(0, 10, 100)
)
```

**Interpolation:**
```python
from scipy import interpolate

# 1D interpolation
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 0.5, 1.0, 1.5, 2.0])

f_linear = interpolate.interp1d(x, y, kind='linear')
f_cubic = interpolate.interp1d(x, y, kind='cubic')

x_new = np.linspace(0, 4, 100)
y_linear = f_linear(x_new)
y_cubic = f_cubic(x_new)

# 2D interpolation
from scipy.interpolate import griddata
points = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
values = np.array([0, 1, 1, 2])
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:100j]
grid_z = griddata(points, values, (grid_x, grid_y), method='cubic')
```

### 3. SymPy - Symbolic Mathematics

**Symbolic Expressions:**
```python
from sympy import symbols, diff, integrate, solve, simplify, expand
from sympy import sin, cos, exp, log, sqrt, pi

# Define symbols
x, y, z = symbols('x y z')
t = symbols('t', real=True, positive=True)

# Create expressions
expr = x**2 + 2*x + 1
simplified = simplify(expr)
expanded = expand((x + 1)**3)

# Differentiation
f = x**3 + 2*x**2 + x
df_dx = diff(f, x)  # 3*x**2 + 4*x + 1
d2f_dx2 = diff(f, x, 2)  # 6*x + 4

# Integration
indefinite = integrate(x**2, x)  # x**3/3
definite = integrate(x**2, (x, 0, 1))  # 1/3

# Solve equations
equation = x**2 - 4
solutions = solve(equation, x)  # [-2, 2]

# System of equations
eq1 = x + y - 5
eq2 = x - y - 1
sol = solve([eq1, eq2], [x, y])  # {x: 3, y: 2}
```

## Complete Examples

### Example 1: Marine Engineering - Catenary Mooring Line

```python
import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

def catenary_mooring_analysis(
    water_depth: float,
    horizontal_distance: float,
    chain_weight: float,  # kg/m
    required_tension: float  # kN
) -> dict:
    """
    Analyze catenary mooring line configuration.

    Parameters:
        water_depth: Water depth (m)
        horizontal_distance: Horizontal distance to anchor (m)
        chain_weight: Chain weight per unit length in water (kg/m)
        required_tension: Required horizontal tension (kN)

    Returns:
        Dictionary with mooring line parameters
    """
    g = 9.81  # m/s²
    w = chain_weight * g / 1000  # Weight per unit length (kN/m)

    # Horizontal tension
    H = required_tension

    # Catenary parameter
    a = H / w

    # Catenary equation: z = a(cosh(x/a) - 1)
    # Need to find length of chain s such that:
    # - Horizontal distance = a*sinh(s/a)
    # - Vertical distance = a(cosh(s/a) - 1) = water_depth

    def equations(s):
        horizontal_eq = a * np.sinh(s/a) - horizontal_distance
        vertical_eq = a * (np.cosh(s/a) - 1) - water_depth
        return [horizontal_eq, vertical_eq]

    # Solve for chain length
    s_initial = np.sqrt(horizontal_distance**2 + water_depth**2)
    s_chain = fsolve(equations, s_initial)[0]

    # Calculate tensions
    T_bottom = H  # Horizontal tension at bottom
    T_top = np.sqrt(H**2 + (w * water_depth)**2)  # Tension at vessel

    # Chain profile
    x_profile = np.linspace(0, horizontal_distance, 100)
    z_profile = a * (np.cosh(x_profile/a) - 1)

    return {
        'chain_length': s_chain,
        'horizontal_tension': H,
        'top_tension': T_top,
        'bottom_tension': T_bottom,
        'catenary_parameter': a,
        'profile_x': x_profile,
        'profile_z': z_profile
    }

# Example usage
result = catenary_mooring_analysis(
    water_depth=1500,  # m
    horizontal_distance=2000,  # m
    chain_weight=400,  # kg/m
    required_tension=2000  # kN
)

print(f"Chain length required: {result['chain_length']:.1f} m")
print(f"Top tension: {result['top_tension']:.1f} kN")
print(f"Bottom tension: {result['bottom_tension']:.1f} kN")
```

### Example 2: Structural Dynamics - Natural Frequency

```python
import numpy as np
from scipy.linalg import eig

def calculate_natural_frequencies(
    mass_matrix: np.ndarray,
    stiffness_matrix: np.ndarray,
    num_modes: int = 5
) -> dict:
    """
    Calculate natural frequencies and mode shapes.

    Solves eigenvalue problem: [K - ω²M]φ = 0

    Parameters:
        mass_matrix: Mass matrix [n×n]
        stiffness_matrix: Stiffness matrix [n×n]
        num_modes: Number of modes to return

    Returns:
        Dictionary with frequencies and mode shapes
    """
    # Solve generalized eigenvalue problem
    eigenvalues, eigenvectors = eig(stiffness_matrix, mass_matrix)

    # Calculate natural frequencies (rad/s)
    omega = np.sqrt(eigenvalues.real)

    # Sort by frequency
    idx = np.argsort(omega)
    omega_sorted = omega[idx]
    modes_sorted = eigenvectors[:, idx]

    # Convert to Hz
    frequencies_hz = omega_sorted / (2 * np.pi)

    # Periods
    periods = 1 / frequencies_hz

    return {
        'natural_frequencies_rad_s': omega_sorted[:num_modes],
        'natural_frequencies_hz': frequencies_hz[:num_modes],
        'periods_s': periods[:num_modes],
        'mode_shapes': modes_sorted[:, :num_modes]
    }

# Example: 3-DOF system
M = np.array([
    [100, 0, 0],
    [0, 100, 0],
    [0, 0, 100]
])  # Mass matrix (kg)

K = np.array([
    [200, -100, 0],
    [-100, 200, -100],
    [0, -100, 100]
])  # Stiffness matrix (kN/m)

result = calculate_natural_frequencies(M, K, num_modes=3)

for i, (f, T) in enumerate(zip(result['natural_frequencies_hz'], result['periods_s'])):
    print(f"Mode {i+1}: f = {f:.3f} Hz, T = {T:.3f} s")
```

### Example 3: Hydrodynamic Analysis - Wave Spectrum

```python
import numpy as np
from