Ecuaciones diferenciales ordinarias

J3.6 · Nivel J3 · Computación científica

Objetivos

  • Programar a mano Euler explícito, Euler implícito y Runge–Kutta, y medir su orden, su estabilidad y su coste.
  • Escribir una ecuación de movimiento como sistema de primer orden y resolverla con ODEProblem y solve.
  • Fijar reltol y abstol sabiendo que no son el error de la solución, y comprobarlo.
  • Reconocer un sistema rígido y elegir método con un argumento.
  • Simular un oscilador forzado, uno no lineal (Duffing) y la respuesta de un mástil al viento.
  • Detectar sucesos con un callback.

Euler explícito

\[\mathbf u_{n+1}=\mathbf u_n+h\,\mathbf f(\mathbf u_n)=(\mathbf I+h\mathbf A)\,\mathbf u_n,\qquad \mathbf A=\begin{bmatrix}0&1\\-\omega_0^2&-2\zeta\omega_0\end{bmatrix}\]

paso_ee(A, u, h) = u + h * (A * u)
  • \(\zeta=0\), \(h=T/40\), cinco periodos: la energía se multiplica por 131.
  • Cada paso: \(E\to(1+(\omega_0h)^2)\,E\), siempre \(>1\).
  • Con \(\zeta=2\ \%\) solo es seguro si \(h<2\zeta/\omega_0=2{,}17\) ms (157 pasos por periodo).

Euler implícito

\[(\mathbf I-h\mathbf A)\,\mathbf u_{n+1}=\mathbf u_n\]

paso_ei(A, u, h) = (I - h * A) \ u
  • Cada paso: \(E\to E/(1+(\omega_0h)^2)\): queda el 0,76 % a los cinco periodos.
  • Es amortiguamiento numérico: confundiría un diagnóstico de daño.
  • Con \(h=0{,}1\) s no explota, pero la oscilación desaparece: estable no es preciso.

Runge–Kutta: orden y coste

método orden medido evaluaciones por paso
Euler explícito 1,0 1
Euler implícito 1,0 (un sistema lineal)
Heun 2,0 2
RK4 4,0 4

Error ≤ 10⁻³ mm en 10 s: Euler explícito 2,3 millones de evaluaciones, Heun 25 600, RK4 3 200, Tsit5 1 407.

Estabilidad por paso

Oscilador sin amortiguar, \(z=\omega_0h\); \(\lvert g\rvert=1\) conserva la energía.

z Euler expl. Euler impl. Heun RK4 trapecio
0,5 1,118 0,894 1,008 1,000 1
1 1,414 0,707 1,118 0,994 1
2 2,236 0,447 2,236 0,745 1
3 3,162 0,316 4,610 1,505 1

RK4 es estable hasta \(z=2\sqrt2\approx2{,}83\) (0,153 s).

Newmark = trapecio

  • Con \(\beta=1/4\), \(\gamma=1/2\) y un sistema lineal coincide con la regla del trapecio (diferencia de \(10^{-14}\) mm).
  • Implícito, de orden 2, conserva la energía con cualquier \(h\).
  • No es exacto: con \(h=0{,}1\) s la frecuencia numérica es el 81 % de la real.

De la ecuación al problema

\[\ddot x+2\zeta\omega_0\dot x+\omega_0^2x=0\;\Longrightarrow\;\begin{cases}\dot x=v\\ \dot v=-2\zeta\omega_0 v-\omega_0^2x\end{cases}\]

function libre!(du, u, p, t)
    du[1] = u[2]
    du[2] = -2p.ζ * p.ω * u[2] - p.ω^2 * u[1]
end
prob = ODEProblem(libre!, [1.0, 0.0], (0.0, 10.0), (ω=ω0, ζ=ζ))
sol  = solve(prob, Tsit5(); reltol=1e-9, abstol=1e-9)

du se modifica en su sitio; sol(t) interpola.

Flujo de trabajo

Cuatro cajas unidas por flechas: ecuación, problema, solución y uso, con el algoritmo y las tolerancias sobre la flecha central

El coste se mide en evaluaciones de \(\mathbf f\) (sol.stats), no en segundos.

La tolerancia no es el error

tolerancia pasos error máximo (mm)
10⁻³ 186 2,8·10⁻³
10⁻⁶ 745 7,2·10⁻⁷
10⁻⁹ 3 001 5,2·10⁻¹⁰
  • Controla el error local; el global se acumula.
  • Con 10⁻³ el error supera la tolerancia.

Conservar la energía

Oscilador sin amortiguamiento, 100 s: \(E(100)/E(0)\) debería valer 1.

tolerancia E(100 s) / E(0)
10⁻³ 1,3131
10⁻⁶ 0,99994
10⁻⁹ 0,9999999

Una simulación sin comprobación no es un resultado.

Oscilador forzado

  • Transitorio \(\propto e^{-\zeta\omega_0t}\) y régimen permanente \(F_0\lvert H(f_e)\rvert\).
  • En resonancia: 3,66 mm \(=F_0/(2\zeta k)\).
  • 95 % a los 8,1 s (\(\tau=2{,}71\) s): medir después del transitorio.
  • La simulación reproduce frf_1gdl con cuatro cifras.

Sistemas rígidos

Edificio de tres plantas + sensor de 5 kHz (razón de rigidez 227 777), tolerancia 10⁻⁸:

método pasos evaluaciones de f
Tsit5 174 021 1 044 189
Rodas5P 2 685 46 254
FBDF 9 481 15 675

Mismo error (≈ 2·10⁻¹⁰ m): cambia el trabajo. El paso lo manda la estabilidad.

Oscilador de Duffing

\[\ddot x+2\zeta\dot x+x+\varepsilon x^3=F\cos\phi(t)\]

  • Frecuencia libre con ContinuousCallback en los cruces por cero.
  • \(\varepsilon=0{,}5\), \(A=2\): simulación 1,57; perturbaciones \(1+\tfrac38\varepsilon A^2=1{,}75\).
  • Una aproximación no se usa fuera de su dominio.

Histéresis

  • Barrido hacia arriba y hacia abajo de \(\Omega\) (\(\varepsilon=1\), \(\zeta=0{,}02\), \(F=0{,}02\)).
  • A \(\Omega=1{,}1\): amplitud 0,44 al subir y 0,09 al bajar.
  • Un sistema lineal da la misma curva en ambos sentidos.

Mástil y vórtices (datos inventados)

  • \(f_s=\mathrm{St}\,U/D\); resonancia en \(U_\mathrm{cr}=f_0D/\mathrm{St}=6\) m/s.
  • A 5 m/s, 0,19 mm; a 6 m/s, 8,2 mm: ×40 con un 20 % más de viento.
  • Simulación y teoría \(F(U)\lvert H(f_s)\rvert\) coinciden a tres cifras.
  • \(\tau=15{,}9\) s: simular 300 s y medir en los últimos 50.

Detectar un umbral

cb = ContinuousCallback((u, t, int) -> abs(u[1]) - umbral,
                        int -> (t1[] = int.t; terminate!(int)))
  • Supera 5 mm a los 15,0 s; 7,5 mm, a los 39,25 s.
  • No se recorre una malla de tiempos, así que no se puede saltar el cruce.

Viento de banda ancha

  • Suma de 176 armónicos (0,5–4 Hz), fases aleatorias con semilla fija, 40 N eficaces.
  • \(\sigma_x^2=\sum_i\tfrac12A_i^2\lvert H(f_i)\rvert^2\) (rms_respuesta_multiseno).
  • Simulado 0,4993 mm; teórico 0,4993 mm.
  • Duplicar \(\zeta\) divide por dos la amplitud en resonancia: 8,2 a 4,1 mm.

La doy por superada si…

  • Resuelvo una ecuación de movimiento con ODEProblem y solve.
  • Explico por qué reltol no es el error y lo compruebo.
  • Reconozco un sistema rígido por su razón de rigidez y sol.stats.
  • Mido el régimen permanente después del transitorio.
  • Detecto un suceso con un callback y estudio una respuesta no lineal.