python-scientific-computing

name	python-scientific-computing
version	1.0.0
description	Python for engineering analysis, numerical computing, and scientific workflows using NumPy, SciPy, SymPy
author	workspace-hub
category	programming
tags	["python","numpy","scipy","sympy","numerical-computing","engineering","scientific-computing"]
platforms	["python"]

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:

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:

# 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:

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:

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:

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:

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

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** + (w * water_depth)**)  

    
    x_profile = np.linspace(, horizontal_distance, )
    z_profile = a * (np.cosh(x_profile/a) - )

     {
        : s_chain,
        : H,
        : T_top,
        : T_bottom,
        : a,
        : x_profile,
        : z_profile
    }


result = catenary_mooring_analysis(
    water_depth=,  
    horizontal_distance=,  
    chain_weight=,  
    required_tension=  
)

()
()
()

Example 2: Structural Dynamics - Natural Frequency

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, , -],
    [, -, ]
])  

result = calculate_natural_frequencies(M, K, num_modes=)

 i, (f, T)  ((result[], result[])):
    ()

Example 3: Hydrodynamic Analysis - Wave Spectrum

import numpy as np
from scipy.integrate import trapz

def jonswap_spectrum(
    frequencies: np.ndarray,
    Hs: float,
    Tp: float,
    gamma: float = 3.3
) -> np.ndarray:
    """
    Calculate JONSWAP wave spectrum.

    Parameters:
        frequencies: Frequency array (Hz)
        Hs: Significant wave height (m)
        Tp: Peak period (s)
        gamma: Peak enhancement factor (default 3.3)

    Returns:
        Spectral density S(f) in m²/Hz
    """
    fp = 1 / Tp  # Peak frequency
    omega_p = 2 * np.pi * fp
    omega = 2 * np.pi * frequencies

    # Pierson-Moskowitz spectrum
    alpha = 0.0081
    beta = 0.74
    S_PM = (alpha * 9.81**2 / omega**5) * np.exp(-beta * (omega_p / omega)**4)

    # Peak enhancement
    sigma = np.where(omega <= omega_p, 0.07, 0.09)
    r = np.exp(-(omega - omega_p)**2 / (2 * sigma**2 * omega_p**2))

    S_JONSWAP = S_PM * gamma**r

    return S_JONSWAP

def wave_statistics(spectrum: np.ndarray, frequencies: np.ndarray) -> dict:
    """
    Calculate wave statistics from spectrum.

    Parameters:
        spectrum: Spectral density S(f)
        frequencies: Frequency array (Hz)

    Returns:
        Wave statistics
    """
    df = frequencies[1] - frequencies[0]

    # Spectral moments
    m0 = trapz(spectrum, frequencies)
    m2 = trapz(spectrum * frequencies**, frequencies)
    m4 = trapz(spectrum * frequencies**, frequencies)

    
    Hm0 =  * np.sqrt(m0)  
    Tz = np.sqrt(m0 / m2)   
    Te = np.sqrt(m0 / m4)  m4 >     

    
    n_waves =  *  / Tz  
    H_max = Hm0 /  * np.sqrt( * np.log(n_waves))

     {
        : Hm0,
        : Tz,
        : Te,
        : H_max,
        : m0,
        : m2
    }


frequencies = np.linspace(, , )
spectrum = jonswap_spectrum(frequencies, Hs=, Tp=, gamma=)
stats = wave_statistics(spectrum, frequencies)

()
()
()

Example 4: Numerical Integration - Velocity to Displacement

import numpy as np
from scipy.integrate import cumtrapz

def integrate_motion_time_history(
    time: np.ndarray,
    acceleration: np.ndarray
) -> dict:
    """
    Integrate acceleration to get velocity and displacement.

    Parameters:
        time: Time array (s)
        acceleration: Acceleration time history (m/s²)

    Returns:
        Dictionary with velocity and displacement
    """
    # Integrate acceleration to get velocity
    velocity = cumtrapz(acceleration, time, initial=0)

    # Integrate velocity to get displacement
    displacement = cumtrapz(velocity, time, initial=0)

    # Remove drift (if needed)
    # Fit linear trend and subtract
    from numpy.polynomial import polynomial as P
    coef_vel = P.polyfit(time, velocity, 1)
    velocity_detrended = velocity - P.polyval(time, coef_vel)

    coef_disp = P.polyfit(time, displacement, 1)
    displacement_detrended = displacement - P.polyval(time, coef_disp)

    return {
        'velocity': velocity,
        'displacement': displacement,
        'velocity_detrended': velocity_detrended,
        'displacement_detrended': displacement_detrended
    }

# Example: Process vessel motion
dt = 0.05  # Time step (s)
duration = 100  # Duration (s)
time = np.arange(0, duration, dt)

# Simulated acceleration (heave motion)
omega = 2 * np.pi / 10  # 10 second period
acceleration = 2 * np.sin(omega * time)

result = integrate_motion_time_history(time, acceleration)

()
()

Example 5: Optimization - Mooring Pretension

from scipy.optimize import minimize
import numpy as np

def optimize_mooring_pretension(
    num_lines: int,
    water_depth: float,
    target_offset: float,
    max_tension: float
) -> dict:
    """
    Optimize mooring line pretensions to achieve target offset.

    Parameters:
        num_lines: Number of mooring lines
        water_depth: Water depth (m)
        target_offset: Target vessel offset (m)
        max_tension: Maximum allowable tension (kN)

    Returns:
        Optimized pretensions
    """
    # Objective: Minimize difference from target offset
    def objective(pretensions):
        # Simplified model: offset = f(pretension)
        # In reality, would use catenary equations
        total_restoring = np.sum(pretensions) / water_depth
        predicted_offset = target_offset / (1 + total_restoring * 0.001)
        return (predicted_offset - target_offset)**2

    # Constraints: All pretensions > 0 and < max_tension
    bounds = [(100, max_tension) for _ in range(num_lines)]

    # Constraint: Symmetric pretensions (pairs equal)
    def constraint_symmetry(pretensions):
        if num_lines % 2 == 0:
            diffs = []
            for i in range(num_lines // 2):
                diffs.append(pretensions[i] - pretensions[i + num_lines//])
             np.array(diffs)
         np.array([])

    constraints = {: , : constraint_symmetry}

    
    x0 = np.ones(num_lines) *   

    
    result = minimize(
        objective,
        x0,
        method=,
        bounds=bounds,
        constraints=constraints
    )

     {
        : result.x,
        : result.success,
        : target_offset,
        : result.message
    }


result = optimize_mooring_pretension(
    num_lines=,
    water_depth=,
    target_offset=,
    max_tension=
)

()
 i, tension  (result[]):
    ()

Example 6: Symbolic Mathematics - Beam Deflection

from sympy import symbols, diff, integrate, simplify, lambdify
from sympy import Function, Eq, dsolve
import numpy as np
import matplotlib.pyplot as plt

def beam_deflection_symbolic():
    """
    Solve beam deflection equation symbolically.

    Beam equation: EI * d⁴y/dx⁴ = q(x)
    """
    # Define symbols
    x, E, I, L, q0 = symbols('x E I L q0', real=True, positive=True)

    # Simply supported beam with uniform load
    # Boundary conditions: y(0) = 0, y(L) = 0, y''(0) = 0, y''(L) = 0

    # Load function
    q = q0  # Uniform load

    # Integrate beam equation four times
    # EI * d⁴y/dx⁴ = q
    # EI * d³y/dx³ = qx + C1 (shear force)
    # EI * d²y/dx² = qx²/2 + C1*x + C2 (bending moment)
    # EI * dy/dx = qx³/6 + C1*x²/2 + C2*x + C3 (slope)
    # EI * y = qx⁴/24 + C1*x³/6 + C2*x²/2 + C3*x + C4 (deflection)

    # For simply supported: C2 = 0, C4 = 0
    # Apply boundary conditions to find C1, C3
    # y(L) = 0: q0*L⁴/24 + C1*L³/6 + C3*L = 0
    # y'(0) = 0 is not required for simply supported

    # Deflection equation
    y = q0*x**2*(L**2 - 2*L*x + x**2)/(24*E*I)

    # Maximum deflection (at center x = L/2)
    y_max = y.subs(x, L/2)
    y_max_simplified = simplify(y_max)

    # Slope
    slope = diff(y, x)

    # Bending moment
    M = -E*I*diff(y, x, )
    M_simplified = simplify(M)

     {
        : y,
        : y_max_simplified,
        : slope,
        : M_simplified
    }


solution = beam_deflection_symbolic()

()
()
()
()


 sympy  symbols
x_sym, E_sym, I_sym, L_sym, q0_sym = symbols()


y_numeric = lambdify(
    (x_sym, E_sym, I_sym, L_sym, q0_sym),
    solution[],
    
)


E_val =   
I_val =    
L_val =      
q0_val =  


x_vals = np.linspace(, L_val, )
y_vals = y_numeric(x_vals, E_val, I_val, L_val, q0_val)

()
()

Best Practices

1. Use Vectorization

# ❌ Slow: Loop
result = []
for x in x_array:
    result.append(np.sin(x) * np.exp(-x))

# ✅ Fast: Vectorized
result = np.sin(x_array) * np.exp(-x_array)

2. Choose Right Data Type

# Use appropriate precision
float32_array = np.array([1, 2, 3], dtype=np.float32)  # Less memory
float64_array = np.array([1, 2, 3], dtype=np.float64)  # More precision

# Use integer when possible
int_array = np.array([1, 2, 3], dtype=np.int32)

3. Avoid Matrix Inverse When Possible

# ❌ Slower and less stable
x = np.linalg.inv(A) @ b

# ✅ Faster and more stable
x = np.linalg.solve(A, b)

4. Use Broadcasting

# Broadcasting allows operations on arrays of different shapes
A = np.array([[1, 2, 3],
              [4, 5, 6]])  # Shape (2, 3)

b = np.array([10, 20, 30])  # Shape (3,)

# Broadcast adds b to each row of A
C = A + b  # Shape (2, 3)

5. Check Numerical Stability

# Check condition number
cond = np.linalg.cond(A)
if cond > 1e10:
    print("Warning: Matrix is ill-conditioned")

# Use appropriate solver for symmetric positive definite
if np.allclose(A, A.T) and np.all(np.linalg.eigvals(A) > 0):
    x = np.linalg.solve(A, b)  # Can use Cholesky internally

Common Patterns

Pattern 1: Load and Process Engineering Data

import numpy as np

# Load CSV data
data = np.loadtxt('../data/measurements.csv', delimiter=',', skiprows=1)

# Extract columns
time = data[:, 0]
temperature = data[:, 1]
pressure = data[:, 2]

# Process
mean_temp = np.mean(temperature)
std_temp = np.std(temperature)
max_pressure = np.max(pressure)

Pattern 2: Solve System of Equations

from scipy.optimize import fsolve

def system(vars):
    x, y, z = vars
    eq1 = x + y + z - 6
    eq2 = 2*x - y + z - 1
    eq3 = x + 2*y - z - 3
    return [eq1, eq2, eq3]

solution = fsolve(system, [1, 1, 1])

Pattern 3: Curve Fitting

from scipy.optimize import curve_fit

def model(x, a, b, c):
    return a * np.exp(-b * x) + c

# Fit data
params, covariance = curve_fit(model, x_data, y_data)
a_fit, b_fit, c_fit = params

Installation

# Using pip
pip install numpy scipy sympy matplotlib

# Using UV (faster)
uv pip install numpy scipy sympy matplotlib

# Specific versions
pip install numpy==1.26.0 scipy==1.11.0 sympy==1.12

Integration with DigitalModel

CSV Data Processing

import numpy as np

# Load OrcaFlex results
data = np.loadtxt('../data/processed/orcaflex_results.csv',
                   delimiter=',', skiprows=1)

# Process time series
time = data[:, 0]
tension = data[:, 1]

# Statistics
max_tension = np.max(tension)
mean_tension = np.mean(tension)
std_tension = np.std(tension)

YAML Configuration Integration

import yaml
import numpy as np

def run_analysis_from_config(config_file: str):
    with open(config_file) as f:
        config = yaml.safe_load(f)

    # Extract parameters
    L = config['geometry']['length']
    E = config['material']['youngs_modulus']

    # Run numerical analysis
    result = calculate_natural_frequency(L, E)

    return result

Resources

• NumPy Documentation: https://numpy.org/doc/
• SciPy Documentation: https://docs.scipy.org/doc/scipy/
• SymPy Documentation: https://docs.sympy.org/
• NumPy for MATLAB Users: https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
• SciPy Lecture Notes: https://scipy-lectures.org/

Performance Tips

1. Use NumPy's built-in functions - They're optimized in C
2. Avoid Python loops - Use vectorization
3. Use views instead of copies when possible
4. Choose appropriate algorithms - O(n) vs O(n²)
5. Profile your code - Find bottlenecks with cProfile

---

Use this skill for all numerical engineering calculations in DigitalModel!