Métodos numéricos de EDO: multipaso, disparo y diferencias finitas

Un resultado numérico de una EDO no queda validado porque la curva sea suave. El algoritmo debe ser consistente con la ecuación, estable en el régimen usado y convergente bajo refinamiento. En problemas de frontera se añade una exigencia: satisfacer simultáneamente la ecuación interior y las condiciones de ambos extremos. Esta unidad desarrolla métodos multipaso para problemas iniciales y dos estrategias complementarias —disparo y diferencias finitas— para problemas de dos puntos.

Recuperación de Euler y Runge–Kutta

Antes de reutilizar valores de pasos anteriores conviene fijar los dos esquemas de un paso que servirán de referencia y de arranque. Para $y'=f(t,y)$, Euler explícito aproxima la pendiente por su valor al inicio del intervalo:

$$y_{n+1}=y_n+h f(t_n,y_n).$$

Su error local es de orden $h^2$ y, bajo las hipótesis usuales de regularidad y estabilidad, el error global es de orden $h$. El método clásico de Runge–Kutta de orden cuatro combina cuatro pendientes:

$$\begin{aligned} k_1&=f(t_n,y_n),\\ k_2&=f(t_n+h/2,y_n+hk_1/2),\\ k_3&=f(t_n+h/2,y_n+hk_2/2),\\ k_4&=f(t_n+h,y_n+hk_3),\\ y_{n+1}&=y_n+\frac h6(k_1+2k_2+2k_3+k_4). \end{aligned}$$

RK4 tiene error global de orden $h^4$ cuando la solución es suficientemente suave. En esta unidad no se toma ese orden como una garantía automática: se estima repitiendo el cálculo con $h$, $h/2$ y $h/4$, manteniendo iguales el problema, el intervalo y el criterio de comparación. Euler ofrece una línea base transparente; RK4 proporciona valores de arranque adecuados para fórmulas multipaso de orden moderado.

Identidad integral y fórmulas multipaso

El problema inicial

$$y'=f(t,y), \qquad y(t_n)=y_n$$

satisface

$$y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}}f(t,y(t))\,dt.$$

Un método de Runge–Kutta aproxima la integral evaluando el campo dentro del paso actual. Un método multipaso reutiliza valores $f_n,f_{n-1},\ldots$ ya calculados. Esto reduce evaluaciones nuevas, pero introduce memoria, necesita valores de arranque y posee una estabilidad gobernada por una ecuación en diferencias.

Para una malla uniforme $t_n=t_0+nh$, se interpola $f$ con los valores en $t_n$ y $t_{n-1}$. Integrar el polinomio lineal sobre el paso produce Adams–Bashforth de dos pasos:

$$y_{n+1}=y_n+\frac h2(3f_n-f_{n-1}).$$

Es explícito y de orden global dos. Adams–Bashforth de tres pasos es

$$y_{n+1}=y_n+\frac h{12} (23f_n-16f_{n-1}+5f_{n-2}),$$

de orden tres. Los coeficientes no son arbitrarios: integran exactamente el interpolante de $f$ del grado correspondiente.

Zill y Cullen sitúan estos métodos después de los esquemas de un paso porque requieren distinguir error local, acumulación global y datos de arranque Zill, 2009 Dennis G. Zill y Michael R. Cullen (2009) Ecuaciones diferenciales con problemas con valores en la frontera 7.ª ed. · Cengage Learning Ubicación consultada: cap. 9, métodos numéricos para problemas de valor inicial .

Adams–Moulton y predictor–corrector

Si el interpolante incluye $f_{n+1}=f(t_{n+1},y_{n+1})$, el método es implícito. El método trapezoidal o Adams–Moulton de un paso es

$$y_{n+1}=y_n+\frac h2(f_n+f_{n+1}).$$

Hay que resolver una ecuación —o un sistema— para $y_{n+1}$. En un predictor–corrector, Adams–Bashforth genera $y_{n+1}^{(0)}$ y después Adams–Moulton lo corrige:

$$y_{n+1}^{(k+1)}=y_n+\frac h2 \left[f_n+f(t_{n+1},y_{n+1}^{(k)})\right].$$

Una sola corrección produce un método práctico, pero no equivale siempre a resolver la ecuación implícita. Si el punto fijo no contrae o si la tolerancia es laxa, el resultado puede perder estabilidad y orden. Newton es más costoso por paso, pero suele ser preferible en sistemas rígidos.

Ejemplo trabajado: AB2 y corrección implícita

Para

$$y'=-y, \qquad y(0)=1,$$

tomemos $h=0.1$ y, para aislar el error multipaso, iniciemos con el valor exacto

$$y_1=e^{-0.1}\approx0.904837418.$$

Adams–Bashforth 2 da

$$y_2=y_1+\frac{0.1}{2} [3(-y_1)-(-y_0)],$$ $$y_2\approx0.819111806.$$

El valor exacto es $e^{-0.2}\approx0.818730753$, por lo que el error es aproximadamente $3.81\times10^{-4}$. Para el trapezoidal implícito,

$$y_2=y_1+0.05(-y_1-y_2),$$ $$y_2=\frac{0.95}{1.05}y_1\approx0.818662426.$$

El error baja a unos $6.83\times10^{-5}$ para este paso. Una comparación seria repite varios $h$; un único punto no demuestra orden ni superioridad general.

Fórmulas BDF y rigidez

Las fórmulas de diferenciación regresiva interpolan $y$ y aproximan su derivada en el extremo nuevo. BDF2 es

$$\frac{3y_{n+1}-4y_n+y_{n-1}}{2h} =f(t_{n+1},y_{n+1}).$$

Es implícito y de orden dos. Su región de estabilidad lo hace útil para muchos problemas rígidos. La rigidez aparece cuando modos rápidos estables obligan a un método explícito a usar pasos mucho menores que la escala de interés. Un método implícito permite pasos mayores por estabilidad, pero la precisión y la resolución de transitorios siguen imponiendo límites.

Para la ecuación de prueba $y'=\lambda y$, un método multipaso genera un polinomio característico dependiente de $z=h\lambda$. En $z=0$, la condición de raíces exige que toda raíz tenga módulo no mayor que uno y que las de módulo uno sean simples. Ésta es la cero-estabilidad. Consistencia más cero-estabilidad conduce a convergencia para problemas regulares. Un método puede tener error local de orden alto y fallar si sus modos parásitos crecen.

Arranque, cambio de paso y memoria

AB2 necesita $y_0,y_1$; AB3 necesita tres estados. Se generan con un método de un paso de orden al menos comparable, usando pasos internos si hace falta. Iniciar AB3 con Euler puede contaminar el error global y ocultar el orden esperado.

Las fórmulas escritas arriba suponen paso constante. Cambiar $h$ sin recalcular coeficientes destruye la interpolación que las originó. Los integradores de producción mantienen una representación de la historia, ajustan orden y paso y controlan errores. Implementar una fórmula fija sirve para aprender y validar, pero no reproduce automáticamente un BDF adaptativo de biblioteca.

Fuchs incluye el refinamiento y el control del integrador como parte de la práctica en sistemas dinámicos; una órbita plausible no sustituye una prueba de orden Fuchs, 2013 Armin Fuchs (2013) Nonlinear Dynamics in Complex Systems: Theory and Applications for the Life-, Neuro- and Natural Sciences Springer Ubicación consultada: apéndice B, pp. 205-212 .

Disparo para problemas de dos puntos

Considere

$$y''=F(x,y,y'), \qquad y(a)=\alpha, \qquad y(b)=\beta.$$

El método de disparo convierte el problema de frontera en una familia de problemas iniciales. Se propone una pendiente desconocida

$$y'(a)=s,$$

se integra hasta $b$ y se define el residuo

$$R(s)=y(b;s)-\beta.$$

La condición derecha se satisface cuando $R(s)=0$. Bisección requiere encerrar un cambio de signo y es robusta; secante evita derivadas; Newton usa

$$s_{k+1}=s_k-\frac{R(s_k)}{R'(s_k)}.$$

La derivada puede aproximarse por diferencias, pero es mejor integrar una ecuación de sensibilidad. Si $z=\partial y/\partial s$, al diferenciar la EDO se obtiene

$$z''=F_yz+F_{y'}z', \qquad z(a)=0, \qquad z'(a)=1,$$

y $R'(s)=z(b)$. Así Newton usa información consistente con la misma trayectoria.

Ejemplo trabajado: disparo lineal

Para

$$-y''=1, \qquad y(0)=0, \qquad y(1)=0,$$

introducimos $y'(0)=s$. La solución del problema inicial es

$$y(x;s)=sx-\frac{x^2}{2}.$$

El residuo derecho vale

$$R(s)=y(1;s)=s-\frac12.$$

Su raíz es $s=1/2$, y recuperamos $y=x(1-x)/2$. Como $R$ es lineal, Newton converge en una iteración desde cualquier valor. En un problema no lineal, $R$ puede tener varias raíces, tangencias o discontinuidades debidas a explosión de la solución inicial. Un solo disparo encontrado no demuestra unicidad.

El disparo simple puede estar mal condicionado en problemas con modos crecientes: un error diminuto en $s$ se amplifica antes de llegar a $b$. El disparo múltiple divide el intervalo, introduce estados en nodos intermedios y exige continuidad; reduce esa amplificación a cambio de un sistema algebraico mayor.

Diferencias finitas para problemas de frontera

En una malla uniforme $x_i=a+ih$, la segunda derivada se aproxima por

$$y''(x_i)=\frac{y_{i-1}-2y_i+y_{i+1}}{h^2}+O(h^2).$$

Para

$$-y''=f(x), \qquad y(a)=\alpha, \qquad y(b)=\beta,$$

los nodos interiores satisfacen

$$-y_{i-1}+2y_i-y_{i+1}=h^2f_i.$$

Los valores de frontera conocidos pasan al lado derecho. El resultado es un sistema tridiagonal. Conviene resolverlo con un algoritmo para matrices banda, no con una inversa densa.

Para el operador conservativo

$$-(p(x)y')'+q(x)y=f(x),$$

se discretiza el flujo en semipuntos:

$$-\frac1h\left[ p_{i+1/2}\frac{y_{i+1}-y_i}{h} -p_{i-1/2}\frac{y_i-y_{i-1}}{h} \right]+q_iy_i=f_i.$$

Esta forma conserva mejor balances y maneja coeficientes variables. Expandir primero las derivadas y discretizar término por término puede perder simetría o conservación.

Ejemplo trabajado: matriz tridiagonal

Retomemos $-y''=1$ en $[0,1]$ con fronteras cero y $h=1/4$. Hay tres incógnitas $y_1,y_2,y_3$:

$$\begin{pmatrix} 2&-1&0\\ -1&2&-1\\ 0&-1&2 \end{pmatrix} \begin{pmatrix}y_1\\y_2\\y_3\end{pmatrix} =\frac1{16} \begin{pmatrix}1\\1\\1\end{pmatrix}.$$

La solución es

$$(y_1,y_2,y_3)= \left(\frac3{32},\frac18,\frac3{32}\right).$$

Coincide exactamente en los nodos con $y=x(1-x)/2$ porque la diferencia central reproduce exactamente la segunda derivada de un polinomio cuadrático. Esta coincidencia especial no demuestra que el esquema sea exacto para otras cargas; una prueba con una solución no polinómica sigue siendo necesaria.

Condiciones de Neumann y Robin

Una condición $y'(a)=\gamma$ puede aproximarse con una diferencia unilateral de primer orden,

$$\frac{y_1-y_0}{h}=\gamma,$$

pero reduciría el orden global cerca del borde. Una fórmula unilateral de segundo orden es

$$\frac{-3y_0+4y_1-y_2}{2h}=\gamma.$$

También pueden usarse nodos fantasma y eliminarse con la condición. El tratamiento de frontera debe tener orden compatible con el interior. En un problema de Neumann puro, la matriz puede ser singular debido al modo constante; se verifica la compatibilidad de la carga y se fija una normalización, por ejemplo media cero.

Zill discute disparo y diferencias finitas como aproximaciones distintas del mismo problema de dos puntos Zill, 2009 Dennis G. Zill (2009) Ecuaciones diferenciales con aplicaciones de modelado 9.ª ed. · Cengage Learning Ubicación consultada: cap. 9, sección de solución numérica de problemas con valores en la frontera . Compararlas en un caso común es una validación cruzada útil.

Solución manufacturada y orden observado

Para validar código se elige una función conocida y se fabrica la carga. Si

$$y_*(x)=\sin\pi x,$$

entonces

$$-y_*''=\pi^2\sin\pi x, \qquad y_*(0)=y_*(1)=0.$$

Se resuelve numéricamente con $h, h/2, h/4,\ldots$ y se calcula

$$E_h=\max_i|y_i-y_*(x_i)|.$$

El orden observado es

$$p_{\mathrm{obs}}= \log_2\frac{E_h}{E_{h/2}}.$$

Para diferencias centrales y fronteras de segundo orden debe acercarse a dos antes de que dominen redondeo o tolerancias. Sin solución exacta puede compararse contra una referencia refinada, declarando que no es verdad exacta.

Protocolo de validación completo

Para un problema inicial multipaso:

1. prueba con una solución exacta no trivial;
2. usa un arranque de orden suficiente;
3. refina el paso sobre el mismo intervalo;
4. registra error máximo, final y costo;
5. comprueba cero-estabilidad y región de estabilidad en ecuaciones de prueba;
6. separa error de arranque del error del método principal.

Para un problema de frontera:

1. calcula residuo de la EDO en los nodos o mediante una reconstrucción;
2. reporta por separado los dos residuos de frontera;
3. refina malla y tolerancia del integrador;
4. compara disparo y diferencias finitas cuando ambos sean apropiados;
5. usa una solución manufacturada;
6. vigila condicionamiento, número de iteraciones y sensibilidad al valor inicial del disparo.

Ejercicios numéricos con evidencia

1. Implementa AB2 y AB3 para $y'=-2y+t$, $y(0)=1$. Inicia con RK4, refina $h$ y estima el orden de cada método.
2. Implementa un predictor AB2 y un corrector trapezoidal. Compara una corrección, iteración hasta tolerancia y solución implícita exacta para $y'=-20y$.
3. Aplica BDF2 a $y'=-100(y-\cos t)-\sin t$, cuya solución exacta puede elegirse como $\cos t$ con el dato correcto. Compara estabilidad y costo con RK4 fijo.
4. Resuelve por disparo $y''=-e^y$, $y(0)=y(1)=0$. Explora varios valores iniciales de la pendiente y documenta si el residuo sugiere más de una raíz.
5. Añade la ecuación de sensibilidad al disparo de un problema no lineal y compara Newton con secante en iteraciones, evaluaciones y robustez.
6. Construye la matriz de diferencias para $-y''+4y=x$ en $[0,1]$ con Dirichlet. Verifica orden dos mediante una solución manufacturada distinta de un polinomio cuadrático.
7. Discretiza $-(p y')'=f$ con $p(x)=1+x$ en forma de flujo. Comprueba simetría y positividad de la matriz bajo fronteras de Dirichlet.
8. Resuelve un mismo problema lineal de frontera por disparo y diferencias finitas. Compara error, residuo, condicionamiento y tiempo; explica en qué régimen preferirías cada método.

Un algoritmo numérico completo incluye método, arranque, resolución algebraica, criterio de parada y evidencia de refinamiento. Multipaso reduce evaluaciones reutilizando historia; disparo ajusta datos iniciales para alcanzar la frontera; diferencias finitas resuelve el perfil entero. La validación determina cuándo cualquiera de esas respuestas merece confianza.