Biblioteca SciPy

Tabla de contenido

Python - Este artículo es parte de una serie.

Parte 1: Por qué Python

Parte 2: Instalación de Python

Parte 3: Primeros pasos en Python

Parte 4: Bibliotecas de uso en ingeniería

Parte 5: Biblioteca NumPy

Parte 6: Biblioteca Matplotlib

Parte 7: Este artículo

Parte 8: Biblioteca Pandas

Parte 9: Biblioteca SymPy

📊 SciPy: Más Allá de los Cálculos Básicos
#

Si bien NumPy es la base para el cálculo numérico, SciPy (abreviatura de Scientific Python) construye sobre ella para resolver problemas científicos y de ingeniería más complejos. SciPy ofrece un vasto conjunto de algoritmos y funciones para optimización, estadísticas, procesamiento de señales, y otras áreas, convirtiéndola en una herramienta indispensable para el análisis avanzado.

Para usar SciPy, generalmente importarás submódulos específicos en lugar de la biblioteca completa. El módulo principal se importa como scipy.

1. Módulos Clave para Ingeniería Mecánica
#

SciPy está organizada en submódulos, cada uno dedicado a un área de la ciencia y la ingeniería. Aquí tienes una descripción detallada de los más relevantes, con sus funciones más comunes y ejemplos prácticos.

Módulo `scipy.integrate`
#

Este módulo es fundamental para la integración numérica, es decir, para calcular el área bajo una curva cuando no se conoce la función exacta o cuando se trabaja con datos discretos de un experimento.

Aplicación: Cálculo del trabajo, momento, energía, o cualquier cantidad que sea la integral de otra.
Funciones principales:
- trapezoid(y, x): Aproxima la integral usando la regla del trapecio. Es la más común para datos experimentales.
- quad(f, a, b): Para integrar funciones matemáticas de forma analítica en un intervalo.

Ejemplo Práctico: Cálculo del Trabajo

import numpy as np
from scipy.integrate import trapezoid

# Datos experimentales de fuerza vs. desplazamiento de un actuador
desplazamiento_m = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.5])
fuerza_N = np.array([0, 20, 45, 60, 80, 100])

# El trabajo es el área bajo la curva de fuerza vs. desplazamiento
trabajo_joules = trapezoid(fuerza_N, desplazamiento_m)
print(f"Trabajo total realizado: {trabajo_joules:.2f} J")

Módulo `scipy.optimize`
#

Este módulo se especializa en la optimización, encontrando los valores que minimizan o maximizan una función. Es crucial para el diseño, el ajuste de modelos y la resolución de sistemas de ecuaciones no lineales.

Aplicación: Optimización de componentes, ajuste de curvas a datos experimentales, minimización de errores en un modelo.
Funciones principales:
- curve_fit(f, xdata, ydata): Ajusta una función a datos experimentales. Es la más utilizada para el ajuste de curvas.
- minimize(fun, x0): Minimiza una función de forma general, útil para problemas de diseño donde buscas el mínimo costo, peso, etc.

Ejemplo Práctico: Ajuste de Curvas

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Datos de un ensayo de resorte
desplazamiento = np.array([0.1, 0.2, 0.3, 0.4, 0.5])
fuerza = np.array([1.1, 2.0, 3.2, 4.1, 5.0])

# Definir la función del modelo (Ley de Hooke: F = k*x)
def ley_hooke(x, k):
    return k * x

# Ajustar los datos al modelo para encontrar k
parametros_optimos, covarianza = curve_fit(ley_hooke, desplazamiento, fuerza)
constante_k = parametros_optimos[0]

# Graficar los resultados
x_modelo = np.linspace(0, 0.6, 100)
y_modelo = ley_hooke(x_modelo, constante_k)

plt.scatter(desplazamiento, fuerza, label='Datos Experimentales')
plt.plot(x_modelo, y_modelo, 'r-', label='Curva Ajustada')
plt.title("Ajuste de Curva para la Ley de Hooke")
plt.xlabel("Desplazamiento (m)")
plt.ylabel("Fuerza (N)")
plt.legend()
plt.grid(True)
plt.show()

Módulo `scipy.interpolate`
#

Este módulo permite la interpolación, es decir, la estimación de valores que se encuentran entre los puntos de datos conocidos. Es una tarea fundamental al trabajar con tablas de datos discretos, como las propiedades de los materiales.

Aplicación: Interpolación de tablas de propiedades de fluidos, curvas de esfuerzo-deformación, o cualquier dato tabulado.
Funciones principales:
- interp1d: Crea una función de interpolación a partir de un conjunto de datos, permitiendo interpolación lineal o cúbica.

Ejemplo Práctico: Estimación de un Valor de Esfuerzo

import numpy as np
from scipy.interpolate import interp1d

# Datos de una curva de esfuerzo-deformación
deformacion_conocida = np.array([0.0, 0.01, 0.02, 0.03, 0.04])
esfuerzo_conocido = np.array([0, 200, 380, 550, 680])

# Crear una función de interpolación lineal
interpolacion = interp1d(deformacion_conocida, esfuerzo_conocido)

# Estimar el esfuerzo para una deformación de 0.025
deformacion_a_estimar = 0.025
esfuerzo_estimado = interpolacion(deformacion_a_estimar)
print(f"El esfuerzo estimado para la deformación {deformacion_a_estimar} es: {esfuerzo_estimado:.2f} MPa")

Módulo `scipy.stats`
#

Este módulo proporciona funciones para el análisis estadístico. Esencial para el control de calidad, el análisis de datos experimentales y la modelización de fenómenos aleatorios.

Aplicación: Pruebas de hipótesis, análisis de distribuciones de datos de fabricación, cálculo de probabilidades.
Funciones principales:
- norm: Proporciona las funciones de la distribución normal (gaussiana).
- t.test: Para pruebas t, que comparan las medias de dos grupos de datos.

Ejemplo Práctico: Análisis de Distribución de Datos

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# Datos de mediciones de la longitud de 500 piezas
mediciones = np.random.normal(loc=150.0, scale=0.5, size=500)

# Obtener el histograma para visualizar la distribución
plt.hist(mediciones, bins=25, density=True, alpha=0.6)

# Ajustar la curva de distribución normal
media, std_dev = norm.fit(mediciones)
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, media, std_dev)
plt.plot(x, p, linewidth=2, label=f'Curva normal ajustada\n media={media:.2f}, std={std_dev:.2f}')

plt.title("Distribución de Medidas de Piezas")
plt.xlabel("Longitud (mm)")
plt.ylabel("Densidad de probabilidad")
plt.legend()
plt.show()

Módulo `scipy.linalg`
#

Este módulo ofrece funciones más avanzadas que el submódulo de NumPy para el álgebra lineal. Es ideal para resolver problemas complejos que involucren matrices de gran tamaño, como en el análisis de elementos finitos o la dinámica estructural.

Aplicación: Descomposición de matrices, cálculo de autovalores y autovectores, resolución de sistemas de ecuaciones a gran escala.
Funciones principales:
- inv(a): Calcula la inversa de una matriz.
- eig(a): Calcula los autovalores y autovectores de una matriz.

Ejemplo Práctico: Cálculo de Autovalores de una Matriz En dinámica de sistemas, los autovalores y autovectores de la matriz de rigidez y masa de un sistema de múltiples grados de libertad se relacionan con las frecuencias naturales y las formas de modo.

import numpy as np
from scipy.linalg import eig

# Matriz de rigidez de un sistema simple de dos masas y dos resortes
K = np.array([[200, -100],
              [-100, 150]])

# Matriz de masa
M = np.array([[10, 0],
              [0, 5]])

# Resolver el problema de autovalores generalizado K*x = lambda*M*x
autovalores, autovectores = eig(K, M)

print("Autovalores (frecuencias al cuadrado):", autovalores.real)
print("Autovectores (formas de modo):\n", autovectores)

Dominar estas funciones de SciPy te permitirá realizar análisis más profundos y sofisticados de tus datos, llevándote más allá de los cálculos básicos de NumPy.

Ejercicios Prácticos de SciPy: ¡Lleva tus Cálculos al Siguiente Nivel!
#

Pon a prueba tu dominio de SciPy con estos ejercicios, diseñados para reforzar los conceptos clave de cada módulo.

1. `scipy.integrate`
#

Cálculo de Energía: Tienes los datos de la fuerza aplicada a un pistón a lo largo de un desplazamiento. Usa trapezoid() para calcular el trabajo total (energía) realizado.

desplazamiento = np.array([0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])

fuerza = np.array([0, 150, 280, 420, 500, 550, 580])

Haz clic para ver la solución

import numpy as np
from scipy.integrate import trapezoid

# Datos de fuerza vs. desplazamiento
desplazamiento = np.array([0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
fuerza = np.array([0, 150, 280, 420, 500, 550, 580])

trabajo = trapezoid(y=fuerza, x=desplazamiento)
print(f"Trabajo total realizado: {trabajo:.2f} J")

Integración de una Curva Teórica: La potencia de un motor puede modelarse con la función \(P(t) = 100 \cdot sin(t) \cdot e^{-0.1t}\). Usa quad() para calcular la energía total generada por el motor en los primeros 10 segundos.
Área bajo una Curva de Esfuerzo-Deformación: Un ensayo de tracción proporciona los siguientes datos. Calcula el módulo de resiliencia (área bajo la curva hasta la fluencia, asumiendo que es el punto final).
- deformacion = np.array([0.00, 0.02, 0.04, 0.06, 0.08])
- esfuerzo = np.array([0, 400, 750, 1000, 1100])
Integración de Datos de Presión: En un sistema hidráulico, la fuerza sobre un émbolo está dada por el área y la presión. Si tienes datos de presión variable con el tiempo, usa trapezoid() para encontrar el impulso total ejercido por la presión.
Flujo a través de una Tubería: El flujo volumétrico a través de una sección transversal se puede calcular integrando el perfil de velocidad. Si el perfil de velocidad es v(r) = 1 - r^2 (donde r es la posición radial), usa quad() para encontrar el flujo total en una tubería de radio 1.

2. `scipy.optimize`
#

Ajuste de la Envolvente de Amplitud: Tienes los datos experimentales de la amplitud máxima (A) de la vibración libre amortiguada de un sistema masa-resorte a lo largo del tiempo (t). Esta envolvente de amplitud sigue la forma \(A(t) = A_0 \cdot e^{-\zeta \omega_n t}\).
- \(\zeta\) (zeta): Constante de amortiguamiento (adimensional).
- \(\omega_n\) (omega-n): Frecuencia natural no amortiguada (en rad/s).
- \(A_0\): Amplitud inicial.
Usa curve_fit() para encontrar los parámetros \(\zeta\) y \(A_0\) que mejor se ajustan a tus datos, asumiendo una frecuencia natural de \(\omega_n = 10\) rad/s.

Datos del experimento:
- tiempo = np.array([0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0])
- amplitud = np.array([10.0, 6.0, 3.5, 2.0, 1.0, 0.5, 0.3, 0.1, 0.05])

Haz clic para ver la solución

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Datos de amplitud vs. tiempo
tiempo = np.array([0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0])
amplitud = np.array([10.0, 6.0, 3.5, 2.0, 1.0, 0.5, 0.3, 0.1, 0.05])

def envolvente_amortiguada(t, zeta, A0):
    wn = 10  # Frecuencia natural no amortiguada en rad/s
    return A0 * np.exp(-zeta * wn * t)

parametros, _ = curve_fit(envolvente_amortiguada, tiempo, amplitud, p0=[0.1, 10])
zeta, A0 = parametros
print(f"Zeta (zeta): {zeta:.4f}")
print(f"Amplitud inicial (A0): {A0:.4f}")

plt.plot(tiempo, amplitud, 'o', label='Datos experimentales')
plt.plot(tiempo, envolvente_amortiguada(tiempo, zeta, A0), 'r-', label='Ajuste')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud')
plt.legend()
plt.show()

Minimización de Costos: Una compañía fabrica un contenedor cilíndrico. El costo del material depende del área de la superficie. Encuentra el radio y la altura que minimizan el área de la superficie para un volumen fijo de 100 m³. Utiliza minimize().
Ajuste de Curva Potencial: El esfuerzo-deformación en la región plástica de un material a menudo se ajusta a la forma \(σ = K \epsilon^n\). Usa curve_fit() para encontrar los parámetros \(K\) y \(n\) a partir de los siguientes datos experimentales.
- deformacion = np.array([0.1, 0.2, 0.3, 0.4])
- esfuerzo = np.array([200, 250, 285, 310])
Optimización del Diseño de una Viga: Una viga debe soportar una carga dada. La deflexión de la viga es una función de su inercia (I). El costo (C) es inversamente proporcional a I. Usa minimize() para encontrar el valor de I que minimiza el costo mientras se mantiene la deflexión por debajo de un umbral máximo.
Ajuste de un Modelo de Carga-Deformación No Lineal: Ajusta los siguientes datos de carga-deformación a una función de la forma \(y = a \cdot x + b \cdot x^2\).
- x = np.array([1, 2, 3, 4, 5])
- y = np.array([3.1, 7.8, 15.2, 25.5, 38.0])

3. `scipy.interpolate`
#

Interpolación de la Densidad de un Fluido: Tienes una tabla de densidad vs. temperatura de un fluido. Usa interp1d() para crear una función que estime la densidad a una temperatura de 150 °C.

temperatura = np.array([100, 200, 300])

densidad = np.array([960, 920, 880])

Haz clic para ver la solución

import numpy as np
from scipy.interpolate import interp1d

# Datos de densidad vs. temperatura
temperatura = np.array([100, 200, 300])
densidad = np.array([960, 920, 880])

# Crear una función de interpolación
interpolacion = interp1d(temperatura, densidad, kind='linear')

# Estimar la densidad a 150 °C
temperatura_a_estimar = 150
densidad_estimada = interpolacion(temperatura_a_estimar)
print(f"La densidad estimada para una temperatura de {temperatura_a_estimar} °C es: {densidad_estimada:.2f} kg/m³")

Estimación de la Viscosidad: Un experimento te proporciona la viscosidad de un lubricante a diferentes temperaturas. Estima la viscosidad a 75 °C utilizando una interpolación lineal.
- temperatura = np.array([50, 100, 150])
- viscosidad = np.array([0.08, 0.05, 0.02])
Interpolación de una Curva de Potencia: Un motor de combustión tiene los siguientes puntos de potencia y velocidad. Estima la potencia a 3500 rpm.
- rpm = np.array([2000, 3000, 4000, 5000])
- potencia = np.array([120, 150, 180, 160])
Cálculo de Coeficientes de Fricción: Utiliza interp1d() con una interpolación cúbica para estimar el coeficiente de fricción a una velocidad de 1.5 m/s, dada una tabla de datos experimentales.
Interpolación de Datos de Sensores: Tienes lecturas de un sensor de presión cada 5 minutos. Crea una función de interpolación para estimar la presión en un momento intermedio (por ejemplo, a los 7.5 minutos).

4. `scipy.stats`
#

Análisis de Control de Calidad: Tienes 200 mediciones de la longitud de un eje. Usa hist() de Matplotlib para visualizar la distribución y norm.fit() para encontrar la media y la desviación estándar de la muestra.

Haz clic para ver la solución

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# Datos de mediciones de la longitud de 200 piezas
mediciones = np.random.normal(loc=150.0, scale=0.5, size=200)

plt.hist(mediciones, bins=25, density=True, alpha=0.6)

media, std_dev = norm.fit(mediciones)
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, media, std_dev)
plt.plot(x, p, linewidth=2, label='Curva normal ajustada\n media={media:.2f}, std={std_dev:.2f}')

plt.title("Distribución de Medidas de Piezas")
plt.xlabel("Longitud (mm)")
plt.ylabel("Densidad de probabilidad")
plt.legend()
plt.show()

Prueba t de Muestras Pares: Un fabricante de bombas prueba un nuevo diseño. Las mediciones de flujo (L/s) antes y después del cambio son las siguientes. Usa ttest_rel() para determinar si el nuevo diseño es significativamente mejor.
- flujo_antes = np.array([10.5, 11.0, 10.8, 10.9])
- flujo_despues = np.array([11.2, 11.5, 11.0, 11.1])
Distribución de Esfuerzos: Los esfuerzos en un componente se distribuyen normalmente con una media de 500 MPa y una desviación estándar de 15 MPa. Usa norm.cdf() para encontrar la probabilidad de que el esfuerzo en una pieza supere los 530 MPa.
Comparación de Tolerancias: Un control de calidad compara las tolerancias de dos lotes de producción. Realiza una prueba t de muestras independientes para ver si hay una diferencia significativa en las medias.
Análisis de la Vida Útil de un Componente: La vida útil de una pieza sigue una distribución de Weibull. Usa el módulo scipy.stats para calcular la vida útil media esperada basándote en un conjunto de datos.

5. `scipy.linalg`
#

Resolución de Sistemas Lineales a Gran Escala: Resuelve un sistema de 10x10 de ecuaciones lineales utilizando una función de scipy.linalg.

Haz clic para ver la solución

import numpy as np
from scipy.linalg import solve
matriz = np.random.rand(10, 10)
vector = np.random.rand(10)

solucion = solve(matriz, vector)
print("Solución del sistema:")
print(solucion)

Autovalores y Frecuencias Naturales: Tienes las matrices de rigidez (K) y masa (M) de un sistema mecánico. Usa eig() para encontrar las frecuencias naturales de vibración.
- K = np.array([[300, -100], [-100, 250]])
- M = np.array([[20, 0], [0, 15]])
Inversa de una Matriz: Calcula la inversa de una matriz de rigidez (K_inv) y utilízala para encontrar el desplazamiento (x) para un vector de fuerza dado (F), donde \(x = K^{-1}F\).
Descomposición Cholesky: La descomposición de Cholesky se usa en la simulación de Monte Carlo. Descompón una matriz de covarianza y explica brevemente su utilidad.
Análisis de Vibración Forzada: Un sistema de dos grados de libertad se somete a una fuerza sinusoidal. La respuesta del sistema se puede encontrar resolviendo un sistema de ecuaciones lineales. Usa scipy.linalg para encontrar la amplitud de la respuesta.