Métodos Numéricos: Interpolación, Regresión, Fourier, EDP Elípticas, Parabólicas e Hiperbólicas

· INTERPOLACIÓN

Interpolación de Newton por diferencias divididas

Polinomio de Newton: f_n(x) = f(x₀) + (x−x₀)·f[x₁,x₀] + (x−x₀)(x−x₁)·f[x₂,x₁,x₀] + … + (x−x₀)…(x−x_n−1)·f[x_n,…,x₀].

Diferencias divididas: f[x_i,x_j] = (f(x_i) − f(x_j)) / (x_i − x_j); f[x₂,x₁,x₀] = (f[x₂,x₁] − f[x₁,x₀]) / (x₂ − x₀); de forma recursiva.

Interpolación lineal: f₁(x) = f(x₀) + ((f(x₁) − f(x₀)) / (x₁ − x₀))·(x − x₀).

Interpolación cuadrática: f₂(x) = b₀ + b₁(x − x₀) + b₂(x − x₀)(x − x₁), con b₀ = f(x₀), b₁ = f[x₁,x₀], b₂ = f[x₂,x₁,x₀].

Error: R_n ≈ f[x_n+1,…,x₀]·(x − x₀)(x − x₁)…(x − x_n).

Newton-Gregory (puntos equiespaciados)

Con Δx constante y α = (x − x₀)/Δx:

f_n = f(x₀) + Δf₀·α + Δ²f₀·α(α−1)/2! + … + Δⁿf₀·α(α−1)…(α−n+1)/n!

Interpolación de Lagrange

f_n(x) = Σᵢ Lᵢ(x)·f(xᵢ), con Lᵢ(x) = Πⱼ≠ᵢ (x − xⱼ)/(xᵢ − xⱼ).

Lineal: f₁ = ((x − x₁)/(x₀ − x₁))·f(x₀) + ((x − x₀)/(x₁ − x₀))·f(x₁).

Cuadrática: f₂ = ((x−x₁)(x−x₂))/((x₀−x₁)(x₀−x₂))·f(x₀) + ((x−x₀)(x−x₂))/((x₁−x₀)(x₁−x₂))·f(x₁) + ((x−x₀)(x−x₁))/((x₂−x₀)(x₂−x₁))·f(x₂).

Newton y Lagrange proporcionan el mismo polinomio.

Splines

Spline cuadrático

fᵢ(x) = aᵢx² + bᵢx + cᵢ en cada intervalo. Con n+1 puntos hay 3n incógnitas. Condiciones: paso por los extremos de cada intervalo (2n), primera derivada continua en los n−1 nudos interiores (n−1) y condición arbitraria a₁ = 0 (primer tramo recta). Sistema efectivo de 3n−1 ecuaciones.

Spline cúbico natural

fᵢ(x) = aᵢx³ + bᵢx² + cᵢx + dᵢ. Condición natural: f″(x₀) = f″(xₙ) = 0.

Sistema tridiagonal en los nudos interiores:

(xᵢ − xᵢ₋₁)·f″(xᵢ₋₁) + 2(xᵢ₊₁ − xᵢ₋₁)·f″(xᵢ) + (xᵢ₊₁ − xᵢ)·f″(xᵢ₊₁) = (6/(xᵢ₊₁ − xᵢ))·[f(xᵢ₊₁) − f(xᵢ)] + (6/(xᵢ − xᵢ₋₁))·[f(xᵢ₋₁) − f(xᵢ)]

Si Δx es constante: f″ᵢ₋₁ + 4f″ᵢ + f″ᵢ₊₁ = (6/Δx²)(fᵢ₋₁ − 2fᵢ + fᵢ₊₁).

Polinomio del tramo (hₓ = xᵢ − xᵢ₋₁): fᵢ(x) = (f″(xᵢ₋₁)/6)(xᵢ − x)³ + (f″(xᵢ)/6)(x − xᵢ₋₁)³ + [f(xᵢ₋₁) − f″(xᵢ₋₁)·h²/6]·(xᵢ − x)/h + [f(xᵢ) − f″(xᵢ)·h²/6]·(x − xᵢ₋₁)/h.

Interpolación bilineal 2D

f(xᵢ, yᵢ) = Σₖ wₖ·fₖ, con Σwₖ = 1.

· REGRESIÓN POR MÍNIMOS CUADRADOS

Regresión lineal

y = a₀ + a₁x, con:

a₁ = (nΣxᵢyᵢ − ΣxᵢΣyᵢ) / (nΣxᵢ² − (Σxᵢ)²); a₀ = ȳ − a₁x̄.

Errores y coeficientes de determinación:

• S_r = Σ(yᵢ − a₀ − a₁xᵢ)²
• S_t = Σ(yᵢ − ȳ)²
• s_y/x = √(S_r/(n−2))
• r² = (S_t − S_r)/S_t
• r = (nΣxy − ΣxΣy) / √((nΣx² − (Σx)²)(nΣy² − (Σy)²))

r² = 1 indica ajuste perfecto; r² = 0 indica que el modelo no aporta información.

Modelos linealizables

• Exponencial y = α·e^βx: ln y frente a x (pendiente β, intercepto ln α).
• Potencial y = α·x^β: log y frente a log x (pendiente β, intercepto log α).
• Saturación y = α·x/(β + x): 1/y frente a 1/x (pendiente β/α, intercepto 1/α).

Regresión polinómica cuadrática

y = a₀ + a₁x + a₂x². Sistema normal 3×3:

[n Σx Σx²; Σx Σx² Σx³; Σx² Σx³ Σx⁴]·[a₀; a₁; a₂] = [Σy; Σxy; Σx²y]

Error estándar: s_y/x = √(S_r/(n − (m+1))).

Regresión lineal múltiple (2 variables)

y = a₀ + a₁x₁ + a₂x₂:

[n Σx₁ Σx₂; Σx₁ Σx₁² Σx₁x₂; Σx₂ Σx₁x₂ Σx₂²]·[a₀; a₁; a₂] = [Σy; Σx₁y; Σx₂y]

Regresión lineal general (matricial)

y = a₀z₀ + a₁z₁ + … + a_mz_m + e, con zⱼ funciones base.

Forma matricial: {Y} = [Z]{A} + {E}.

Ecuaciones normales: [Z]ᵀ[Z]{A} = [Z]ᵀ{Y}. Solución: {A} = ([Z]ᵀ[Z])⁻¹·[Z]ᵀ{Y}.

Varianzas: var(aⱼ) = z⁻¹ⱼⱼ·s²_y/x. Intervalo de confianza: aⱼ ± t_{α/2, n−2}·s(aⱼ).

Regresión no lineal · Gauss-Newton

Función no lineal en los parámetros, por ejemplo f(x) = a₀(1 − e^−a₁x).

Desarrollo de Taylor de primer orden: f(xᵢ)_j+1 ≈ f(xᵢ)_j + (∂f/∂a₀)Δa₀ + (∂f/∂a₁)Δa₁.

Definiendo el vector de residuos Dᵢ = yᵢ − f(xᵢ)_j y el jacobiano [Zⱼ]_ik = ∂f(xᵢ)/∂aₖ:

Sistema linealizado: {D} = [Zⱼ]{ΔA} + {E}. Ecuaciones normales: [Zⱼ]ᵀ[Zⱼ]{ΔA} = [Zⱼ]ᵀ{D}.

Actualización: aₖ^j+1 = aₖ^j + Δaₖ. Criterio de parada: |Δaₖ / aₖ^j+1|·100% < ε_s.

Para f = a₀(1 − e^−a₁x): ∂f/∂a₀ = 1 − e^−a₁x; ∂f/∂a₁ = a₀·x·e^−a₁x.

· SERIES Y TRANSFORMADA DE FOURIER

Sinusoide

f(t) = A₀ + C₁cos(ω₀t + θ) = A₀ + A₁cos(ω₀t) + B₁sen(ω₀t), con A₁ = C₁cosθ, B₁ = −C₁senθ.

C₁ = √(A₁² + B₁²); θ = arctan(−B₁/A₁); ω₀ = 2πf = 2π/T.

Ajuste por mínimos cuadrados de una sinusoide con N puntos equiespaciados

A₀ = Σy/N; Aⱼ = (2/N)Σy·cos(jω₀t); Bⱼ = (2/N)Σy·sen(jω₀t).

Serie continua de Fourier

Para una función T-periódica con ω₀ = 2π/T:

f(t) = a₀ + Σₖ₌₁^∞ [aₖcos(kω₀t) + bₖsen(kω₀t)]

Coeficientes mediante integrales:

• a₀ = (1/T)∫₀ᵀ f(t) dt
• aₖ = (2/T)∫₀ᵀ f(t)·cos(kω₀t) dt
• bₖ = (2/T)∫₀ᵀ f(t)·sen(kω₀t) dt

Simetrías: si f es par, f(−t) = f(t) ⇒ todos los bₖ = 0; si f es impar, f(−t) = −f(t) ⇒ todos los aₖ = 0 (incluido a₀).

Forma exponencial compleja: f(t) = Σₖ₌₋∞^∞ c̃ₖ·e^ikω₀t, con c̃ₖ = (1/T)∫_−T/2^T/2 f(t)·e^−ikω₀t dt; c̃₀ = a₀; c̃ₖ = (aₖ − ibₖ)/2.

Transformada de Fourier (aperiódica, frecuencia continua)

F(iω) = ∫_−∞^∞ f(t)·e^−iωt dt; inversa: f(t) = (1/2π)∫_−∞^∞ F(iω)·e^iωt dω.

Serie ↔ señales periódicas y espectro discreto; Transformada ↔ señales aperiódicas y espectro continuo.

DFT (N puntos discretos, ω₀ = 2π/N)

Fₖ = Σₙ₌₀^{N−1} fₙ·e^−ikω₀n, para k = 0,…,N−1.

Inversa: fₙ = (1/N)Σₖ₌₀^{N−1} Fₖ·e^ikω₀n. Coste ∼ N².

FFT (N = 2^M)

Coste ∼ N·log₂N.

Sande-Tukey (decimación en frecuencia), con W = e^−i·2π/N:

gₙ = fₙ + f_n+N/2; hₙ = (fₙ − f_n+N/2)·Wⁿ ⇒ F_2k = DFT{gₙ} y F_2k+1 = DFT{hₙ}.

Se aplica recursivamente hasta DFT de 2 puntos (mariposas o butterflies). El resultado sale en orden bit-reversal.

Cooley-Tukey (decimación en tiempo): mismo coste.

Teorema de Parseval: (1/T)∫_−T/2^T/2 f²(t) dt = Σₖ |Fₖ|².

· EDP ELÍPTICAS

Clasificación de EDP de 2.º orden

Para A·u_xx + B·u_xy + C·u_yy + D = 0: B² − 4AC < 0 elíptica (Laplace); = 0 parabólica (calor); > 0 hiperbólica (onda).

Laplace y Poisson

Laplace: ∂²T/∂x² + ∂²T/∂y² = 0. Poisson: ∂²T/∂x² + ∂²T/∂y² = f(x,y).

Discretización con diferencias centrales (Δx = Δy): T_i+1,j + T_i−1,j + T_i,j+1 + T_i,j−1 − 4T_i,j = 0.

Métodos iterativos

• Jacobi (explícito iterativo): T_i,j^k+1 = (T_i+1,j^k + T_i−1,j^k + T_i,j+1^k + T_i,j−1^k)/4. Usa solo valores antiguos (iteración k). Lento.
• Gauss-Seidel / Liebmann (implícito iterativo): T_i,j = (T_i+1,j + T_i−1,j + T_i,j+1 + T_i,j−1)/4 usando valores nuevos conforme se calculan. Converge siempre (matriz diagonalmente dominante) y es más rápido que Jacobi.
• SOR (sobrerrelajación): T_i,j^new = λ·T_i,j^GS + (1 − λ)·T_i,j^old, con 1 < λ < 2. Acelera la convergencia. λ = 1 equivale a Gauss-Seidel; el λ óptimo depende del problema.

Criterio de parada: |(T^new − T^old)/T^new|·100% < ε_s.

Flujos de calor (variables secundarias, k' = k·ρ·C)

• q_x = −k'(T_i+1,j − T_i−1,j)/(2Δx)
• q_y = −k'(T_i,j+1 − T_i,j−1)/(2Δy)
• Módulo: q_n = √(q_x² + q_y²)
• Dirección: θ = arctan(q_y/q_x) (sumar π si q_x < 0)

Frontera de Neumann

Con ∂T/∂x dada, nodo ficticio T_−1,j = T_1,j − 2Δx·(∂T/∂x). Sustituyendo: 2T_1,j − 2Δx·(∂T/∂x) + T_0,j+1 + T_0,j−1 − 4T_0,j = 0. Frontera aislada: ∂T/∂x = 0.

· EDP PARABÓLICAS

k·∂²T/∂x² = ∂T/∂t; λ = k·Δt/Δx².

Explícito 1D

T_i^l+1 = T_i^l + λ(T_i+1^l − 2T_i^l + T_i−1^l).

Estabilidad: λ ≤ 1/2 (estable); λ ≤ 1/4 (no oscila); λ = 1/6 (mínimo error de truncamiento). Si λ > 1/2 diverge.

Explícito 2D

T_i,j^l+1 = T_i,j^l + λ[(T_i+1,j + T_i−1,j − 2T_i,j) + (T_i,j+1 + T_i,j−1 − 2T_i,j)]^l. Estabilidad: λ ≤ 1/4 (con Δx = Δy).

Implícito simple 1D

Resolver un sistema tridiagonal en cada paso: −λT_i−1^l+1 + (1 + 2λ)T_i^l+1 − λT_i+1^l+1 = T_i^l.

Con bordes f₀(t) y f_m+1(t): (1 + 2λ)T₁^l+1 − λT₂^l+1 = T₁^l + λ·f₀(t_l+1); −λT_m−1^l+1 + (1 + 2λ)T_m^l+1 = T_m^l + λ·f_m+1(t_l+1).

Incondicionalmente estable. Precisión O(Δt), O(Δx²).

Implícito 2D

−λT_i−1,j^l+1 − λT_i+1,j^l+1 + (1 + 4λ)T_i,j^l+1 − λT_i,j−1^l+1 − λT_i,j+1^l+1 = T_i,j^l.

Matriz pentadiagonal (5 diagonales no nulas), muy costosa; por ello se usa ADI.

Crank-Nicolson 1D

Promedio temporal de la segunda derivada espacial:

−λT_i−1^l+1 + 2(1 + λ)T_i^l+1 − λT_i+1^l+1 = λT_i−1^l + 2(1 − λ)T_i^l + λT_i+1^l.

Sistema tridiagonal. Incondicionalmente estable, con precisión O(Δt²) y O(Δx²). Es el estándar para parabólicas 1D.

Crank-Nicolson 2D

−(λ/2)·[T_i−1,j + T_i+1,j + T_i,j−1 + T_i,j+1]^l+1 + (1 + 2λ)T_i,j^l+1 = (λ/2)·[T_i−1,j + T_i+1,j + T_i,j−1 + T_i,j+1]^l + (1 − 2λ)T_i,j^l.

Sigue siendo pentadiagonal, motivo por el que se emplea ADI.

ADI (Alternating-Direction Implicit, 2D)

Cada paso temporal se divide en dos medios pasos alternando la dirección implícita, manteniendo siempre matrices tridiagonales.

Primer medio paso t_l → t_l+1/2 (x explícita, y implícita):

−λT_i,j−1^l+1/2 + 2(1 + λ)T_i,j^l+1/2 − λT_i,j+1^l+1/2 = λT_i−1,j^l + 2(1 − λ)T_i,j^l + λT_i+1,j^l.

Segundo medio paso t_l+1/2 → t_l+1 (y explícita, x implícita):

−λT_i−1,j^l+1 + 2(1 + λ)T_i,j^l+1 − λT_i+1,j^l+1 = λT_i,j−1^l+1/2 + 2(1 − λ)T_i,j^l+1/2 + λT_i,j+1^l+1/2.

Incondicionalmente estable, precisión O(Δt²) y O(Δx²), y solo requiere sistemas tridiagonales (eficientes).

· EDP HIPERBÓLICAS

∂²u/∂t² = c²·∂²u/∂x²; número de Courant r = c·Δt/Δx.

Explícito con diferencias centrales

u_i^l+1 = 2u_i^l − u_i^l−1 + r²·(u_i+1^l − 2u_i^l + u_i−1^l).

Estabilidad (condición CFL): r ≤ 1. Con r = 1 proporciona la solución exacta; con r > 1 diverge exponencialmente.

Requiere dos condiciones iniciales: posición u(x,0) = f(x) y velocidad ∂u/∂t(x,0) = g(x).

Primer paso (la fórmula general requiere u⁻¹):

u_i¹ = u_i⁰ + Δt·g(x_i) + (r²/2)·(u_i+1⁰ − 2u_i⁰ + u_i−1⁰).

A partir del segundo paso se usa la fórmula general.

Implícito tipo Crank-Nicolson

(u_i^l+1 − 2u_i^l + u_i^l−1)/Δt² = (c²/2)·[(δ_x²u^l+1 + δ_x²u^l−1)/Δx²].

Genera un sistema tridiagonal en cada paso. Incondicionalmente estable (sin CFL), pero más costoso.

Lax-Wendroff (advección u_t + a·u_x = 0, 2.º orden)

u_i^l+1 = u_i^l − (aΔt/(2Δx))(u_i+1^l − u_i−1^l) + ((aΔt)²/(2Δx²))(u_i+1^l − 2u_i^l + u_i−1^l).

Upwind (a > 0)

u_i^l+1 = u_i^l − (aΔt/Δx)(u_i^l − u_i−1^l).

Ambos métodos exigen la condición CFL |aΔt/Δx| ≤ 1.