Euler y el significado de aproximar una evolución continua

Del campo de velocidades a una trayectoria calculada

En los capítulos anteriores trabajamos con modelos discretos, en los que el estado se actualiza por rondas. Una ecuación como $x_{n+1}=F(x_n)$ especifica directamente qué estado sigue al anterior. Si el modelo está definido así, no hace falta inventar estados intermedios entre dos rondas. En el capítulo 23 vimos incluso que unir los puntos con segmentos no convierte esa dinámica discreta en una evolución continua.

Ahora el problema cambia. Consideremos una ecuación diferencial ordinaria de primer orden,

La función $f$ no entrega directamente el estado al cabo de una hora, de un día o de cualquier otro intervalo finito. Especifica la tasa instantánea de cambio en cada punto donde se evalúa. Si existe una solución, debe ser compatible con esas pendientes a lo largo de todo su recorrido. Algunas ecuaciones admiten una expresión analítica cerrada, mientras que muchas otras exigen aproximar la trayectoria mediante procedimientos numéricos.

El método de Euler explícito es la forma más simple de hacer esa aproximación. Rara vez sería el integrador elegido para un problema avanzado. Su valor en este libro es que hace visibles, con muy poca maquinaria, distinciones que acompañarán todo el estudio de los métodos numéricos. Separa el modelo del algoritmo y el error local del global. También permite distinguir la exactitud de la estabilidad, así como las propiedades de la ecuación de las de su discretización. Los textos de análisis numérico lo usan como caso básico para mostrar por qué la consistencia de un método necesita la estabilidad para asegurar la convergencia (Hairer et al., 1993; Iserles, 2008).

En todo el capítulo supondremos, salvo indicación contraria, que el problema de valor inicial está bien planteado y tiene solución única en el intervalo considerado. Así evitaremos atribuir al integrador dificultades que provengan de una ecuación sin solución única o de un modelo mal especificado.

Derivación desde la forma integral

Para una solución suficientemente regular, la ecuación diferencial equivale a

Sea $h=t_{n+1}-t_n>0$ el tamaño del paso. La dificultad está en la integral, porque calcularla exactamente exigiría conocer $x(s)$ dentro del intervalo, que es justamente lo que se busca. Euler explícito reemplaza la tasa variable dentro del intervalo por la tasa evaluada en su extremo izquierdo, es decir, $f(s,x(s))\approx f(t_n,x_n)$. La integral vale entonces aproximadamente $h f(t_n,x_n)$ y se obtiene la regla

La misma fórmula aparece al truncar el desarrollo de Taylor después del término lineal. Si la solución tiene dos derivadas continuas, existe un $\xi_n$ entre $t_n$ y $t_n+h$ tal que

Como $x'(t_n)=f(t_n,x(t_n))$, el método conserva los dos primeros términos y descarta el resto.

Llamaremos defecto local de un paso a la diferencia que resulta al insertar la solución exacta en una actualización de Euler. Bajo las condiciones anteriores, ese defecto es de orden $O(h^2)$. Algunos textos definen el error local de truncamiento dividiendo el defecto por $h$, de modo que lo expresan como $O(h)$. Ambas convenciones describen el mismo hecho. Usaremos la primera para que no parezca contradictorio afirmar después que el error global de Euler es de primer orden.

El procedimiento tiene una lectura geométrica. Desde $(t_n,x_n)$ se avanza durante un tiempo $h$ siguiendo la pendiente que el campo asigna al comienzo del intervalo, mientras que la solución verdadera cambia de pendiente a medida que avanza. Reducir $h$ permite actualizar esa pendiente con mayor frecuencia.

La afirmación de que un defecto local de orden $h^2$ produce un error global de orden $h$ puede precisarse. Sea $e_n=x(t_n)-x_n$ el error global y supongamos que $f$ es Lipschitz respecto de $x$, con constante $L$, en la región recorrida. Si además $|x''(t)|\leq M$, el defecto local satisface $|d_{n+1}|\leq Mh^2/2$. Al restar la actualización numérica de la identidad exacta se obtiene

Por la condición de Lipschitz,

Si $e_0=0$ y $t_n-t_0=nh\leq T$, la iteración de esta desigualdad da, para $L>0$,

El factor que multiplica a $h$ depende del horizonte y de cuánto amplifica el campo las diferencias entre trayectorias. La conclusión $O(h)$ exige, por lo tanto, algo más que contar unos $T/h$ defectos locales, porque también hay que controlar cómo se propagan. Esta observación anticipa la relación entre estabilidad y convergencia. En el caso $L=0$, la suma directa de los defectos da una cota proporcional a $Th$. Estas cotas suelen ser conservadoras. En el decaimiento de la sección siguiente el campo contrae las diferencias y el error real queda muy por debajo de la cota, como muestra el ejercicio 24.2.

Un problema cuya solución conocemos exactamente

Tomemos el decaimiento lineal

cuya solución exacta es $x(t)=x_0e^{-kt}$. Si $x_0>0$, la solución permanece positiva y decrece monótonamente hacia cero. Euler produce $x_{n+1}=(1-hk)x_n$ y, por iteración,

La comparación permite separar con claridad el comportamiento de la ecuación del que introduce el algoritmo.

Con $k=1$ y $h=0{,}5$, el multiplicador de Euler es $0{,}5$ y la secuencia es $1$; $0{,}5$; $0{,}25$; $0{,}125$. En esos mismos tiempos la solución exacta vale $1$; $e^{-0{,}5}\approx0{,}607$; $e^{-1}\approx0{,}368$; $e^{-1{,}5}\approx0{,}223$. El método sigue la tendencia correcta, aunque decae demasiado rápido.

Con $h=1{,}5$, en cambio, el multiplicador es $1-h=-0{,}5$ y Euler genera $1$; $-0{,}5$; $0{,}25$; $-0{,}125$. La magnitud disminuye, pero el signo alterna. Esa oscilación no pertenece a la ecuación $x'=-x$, cuya solución positiva nunca cruza el cero. Es un artefacto numérico que produce el tamaño del paso.

El ejemplo deshace una confusión frecuente. Una trayectoria numérica puede mostrar un rasgo cualitativo que no existe en el modelo continuo. Quien examinara solo la salida del programa, sin conocer el algoritmo, podría atribuir al sistema una oscilación que introdujo la discretización.

Error local y acumulación del error global

En el decaimiento con $k=1$, la evolución exacta durante un paso multiplica por

mientras que Euler multiplica por $1-h$. La discrepancia de un paso comienza, por tanto, en $h^2/2$. Para llegar a un tiempo fijo $T$, sin embargo, se necesitan unos $N=T/h$ pasos. La acumulación de muchos defectos locales produce un error global que, bajo hipótesis de suavidad y estabilidad, es de orden $O(h)$.

Por eso Euler es un método de primer orden. La clasificación es asintótica, porque describe el comportamiento del error cuando $h$ es suficientemente pequeño. No garantiza que dividir por dos un paso concreto divida exactamente el error por dos. El cuadro 24.1 reúne los resultados para $x'=-x$ con $x(0)=1$ y $T=2$.

Cuadro 24.1. Método de Euler para $x'=-x$ con $x(0)=1$ en $T=2$, cuyo valor exacto es $e^{-2}\approx0{,}135335283$. El cociente divide el error de la fila anterior por el de la fila actual.

Con los pasos más pequeños, dividir $h$ por dos divide el error por un factor cercano a dos. En una representación log-log, un error proporcional a $h$ se aproxima a una recta de pendiente uno (figura 24.2).

Esa precisión tiene un costo. Para un horizonte fijo, dividir $h$ por dos duplica aproximadamente el número de evaluaciones de $f$. Un método de mayor orden puede alcanzar una precisión dada con menos pasos, aunque cada paso cueste más. El capítulo 25 hará esa comparación.

Reducir indefinidamente el paso tampoco mejora sin límite el resultado en aritmética finita. El error de truncamiento disminuye, pero crece el número de operaciones y los errores de redondeo pueden acumularse (Higham, 2002). En este capítulo trabajamos lejos de ese régimen. Un paso más pequeño mejora el resultado mientras el error de truncamiento domine sobre el de redondeo.

Estabilidad absoluta y decaimiento numérico

Para estudiar la estabilidad numérica se usa la ecuación de prueba $y'=\lambda y$, con $\lambda\in\mathbb C$. Cuando $\operatorname{Re}(\lambda)<0$, la solución exacta $y(t)=y_0e^{\lambda t}$ decae en magnitud. Euler produce

La función $R$ se denomina función o factor de amplificación de Euler explícito. Para que la aproximación tienda a cero en este problema lineal se necesita $|1+h\lambda|<1$. El conjunto de los $z$ que cumplen $|1+z|<1$ es el interior del círculo de radio uno centrado en $-1$ y forma la región de estabilidad absoluta del método para la ecuación de prueba (Ascher y Petzold, 1998; Butcher, 2016; Iserles, 2008).

Con $\lambda=-1$, la condición se reduce a $|1-h|<1$, es decir, a $0<h<2$. En $h=2$ el factor vale $-1$, así que la magnitud no disminuye y la solución numérica alterna indefinidamente. Para $h>2$ la magnitud crece, aunque la solución exacta decaiga.

Con $h=1{,}5$ el método es estable en ese sentido, porque $|-0{,}5|<1$. Aun así, produce valores negativos y oscilaciones que la solución exacta no tiene. La estabilidad impide el crecimiento artificial de ese modo, pero no garantiza la positividad ni la precisión.

En el decaimiento escalar, la no negatividad exige $1-hk\geq0$, es decir, $hk\leq1$. Si $hk=1$, la aproximación llega a cero en un solo paso, de modo que conservar la positividad estricta de la solución exacta exige $hk<1$. Ambas condiciones son más estrictas que la estabilidad absoluta, $0<hk<2$. En problemas aplicados, la propiedad que interesa puede restringir el paso aunque el integrador sea formalmente estable.

Estabilidad del método y estabilidad del sistema

El capítulo 10 examinó la estabilidad de un sistema frente a perturbaciones. Aquí aparece una noción relacionada pero distinta, la estabilidad numérica de un esquema. Si la ecuación continua tiene un equilibrio estable y el paso es demasiado grande, Euler puede generar una secuencia divergente. Esa divergencia no demuestra que el equilibrio continuo sea inestable. Indica que el método, con ese paso, no reproduce bien la dinámica cerca del equilibrio.

También puede ocurrir lo contrario. Euler implícito, que estudiaremos en el capítulo 26, reduce paso a paso la amplitud de un oscilador sin rozamiento, aunque el modelo la conserva. Por eso no deben inferirse propiedades del sistema original a partir de una implementación sin haber analizado el algoritmo.

La distinción abarca varios niveles. La estabilidad del modelo continuo depende de sus ecuaciones, mientras que la del método numérico se define respecto de ciertos problemas de prueba y de ciertos pasos. Una implementación puede, además, volverse inestable por errores de código o por efectos de la aritmética finita.

En problemas no lineales, el análisis suele combinar linealizaciones locales con experimentos de refinamiento del paso. Ninguno de esos procedimientos sustituye la comprensión del modelo. Un resultado estable tampoco garantiza que el método haya captado todas las propiedades cualitativas relevantes.

El modelo de dos compartimentos visto como Euler

Volvamos al modelo discreto de dos compartimentos del capítulo 21,

Consideremos ahora un modelo continuo de intercambio,

con tasas $a$ y $b$ en unidades de tiempo inverso. En forma matricial, $\dot z=Az$ con

Un paso de Euler de tamaño $h$ produce $z_{n+1}=(I+hA)z_n$, es decir,

Si se identifican $\alpha=ha$ y $\beta=hb$, se obtiene exactamente la forma algebraica del modelo discreto.

La coincidencia formal no obliga a dar a ambos la misma interpretación. En un modelo intrínsecamente discreto, $\alpha$ puede ser la fracción transferida durante una ronda que el propio sistema define. En el modelo continuo discretizado, $\alpha=ha$ aparece porque se aproximó una tasa continua durante un intervalo numérico $h$.

La diferencia se aclara al refinar el paso. Para aproximar el mismo modelo continuo con la mitad del paso, $\alpha=ha$ y $\beta=hb$ también deben reducirse a la mitad. Si $\alpha$ y $\beta$ se mantienen fijos mientras se declara que $h$ cambió, cambian las tasas implícitas $a=\alpha/h$ y $b=\beta/h$, de modo que ya no se refina la aproximación del mismo sistema.

Los capítulos 21 y 23 dan un ejemplo concreto. Con $a=0{,}2$ y $b=0{,}1$ por día, el paso $h=1$ día reproduce la configuración A, con fracciones $(0{,}20;\ 0{,}10)$. El paso $h=0{,}5$ día reproduce la configuración B, con fracciones $(0{,}10;\ 0{,}05)$, siempre que cada ronda de B dure medio día. Vistas así, A y B dejan de ser dos modelos distintos y pasan a ser dos aproximaciones del mismo modelo continuo. Durante un día, la desviación respecto del equilibrio se multiplica por $0{,}7$ con el primer paso y por $0{,}85^2=0{,}7225$ con el segundo, mientras que la solución exacta la multiplica por $e^{-0{,}3}\approx0{,}7408$. Al dividir el paso por dos, el error de ese factor pasa de $0{,}0408$ a $0{,}0183$.

Euler conserva exactamente el total, porque $d(x+y)/dt=0$ y, en la actualización,

Una aproximación de primer orden puede preservar exactamente un invariante lineal, aunque eso no significa que Euler preserve todas las propiedades del flujo continuo.

La matriz $I+hA$ transforma todo estado no negativo en otro no negativo si y solo si sus entradas son no negativas, es decir, si $ha\leq1$ y $hb\leq1$. El modo no conservado tiene valor propio continuo $-(a+b)$ y multiplicador de Euler $1-h(a+b)$, cuyo decaimiento exige $|1-h(a+b)|<1$. Conservar el total no asegura la estabilidad de ese modo ni la no negatividad. Cada una de esas propiedades debe comprobarse por separado.

Paso numérico e intervalo de observación

En el capítulo 23 vimos que una trayectoria bien calculada puede guardarse con una frecuencia que oculta parte de su comportamiento. El paso del integrador y el intervalo de salida no tienen por qué coincidir.

Podríamos integrar una ecuación continua con un paso interno de $h=0{,}01$ días y guardar una observación por día. El archivo resultante tendría un intervalo de salida de un día, aunque el algoritmo habría dado cien pasos internos entre dos filas consecutivas. A la inversa, disponer de datos diarios no autoriza a usar un paso numérico de un día, porque la estabilidad y la exactitud dependen de las escalas temporales del modelo. La frecuencia con que se midieron los datos es otra cuestión.

Tampoco hay que identificar automáticamente $h$ con un proceso físico de muestreo. En Euler, $h$ es una decisión numérica para aproximar una ecuación continua. El modelo puede contener además retardos o mediciones discretas con escalas temporales propias.

La separación importa especialmente en modelos clínicos o biológicos. Un proceso puede combinar una variable rápida con otra lenta, mientras las mediciones se hacen con mucha menor frecuencia que ambas. El integrador debe resolver la dinámica que especifican las ecuaciones. La frecuencia de observación responde a otra pregunta, la de qué información hay disponible para contrastar el modelo.

Implementación mínima y prueba de refinamiento

Una implementación vectorial de Euler puede escribirse así.

import numpy as np


def euler(f, t0, x0, h, pasos):
    if not (np.isfinite(h) and h > 0):
        raise ValueError("h debe ser un número finito y positivo.")
    if (isinstance(pasos, (bool, np.bool_))
            or not isinstance(pasos, (int, np.integer))
            or pasos < 0):
        raise ValueError("pasos debe ser un entero no negativo.")

    x0 = np.array(x0, dtype=float, copy=True, ndmin=1)
    t = t0 + h * np.arange(pasos + 1)
    x = np.empty((pasos + 1, x0.size), dtype=float)
    x[0] = x0

    for n in range(pasos):
        pendiente = np.asarray(f(t[n], x[n]), dtype=float)
        if pendiente.shape != x0.shape:
            raise ValueError("f devolvió una forma incompatible.")
        x[n + 1] = x[n] + h * pendiente

    return t, x


t, x = euler(lambda t, x: -x, 0.0, 1.0, 0.25, 8)
print(t[-1], x[-1, 0])

La función comprueba que el paso sea finito y positivo y que el número de pasos sea un entero. Los tiempos se calculan como $t_0+nh$, porque sumar $h$ repetidamente acumula errores de redondeo. Con $h=0{,}1$, diez sumas sucesivas dan $0{,}9999999999999999$, mientras que el producto $10\times0{,}1$ da exactamente 1. Las dos últimas líneas aplican el método al decaimiento con $h=0{,}25$ e imprimen 2.0 0.1001129150390625, el valor del cuadro 24.1 con todas sus cifras.

Para evaluar un integrador no basta con ejecutarlo una vez y observar que la figura parece razonable. Cuando existe una solución conocida, hay que compararla directamente. Cuando no existe, se hacen estudios de refinamiento, que repiten el cálculo con pasos $h$, $h/2$ y $h/4$ y comprueban si las diferencias disminuyen según el orden esperado, dentro de un régimen donde domina el error de discretización.

En un método de orden $p$, el error se comporta como $E(h)\approx Ch^p$. Si se conoce la solución, el cociente $E(h)/E(h/2)$ se aproxima a $2^p$, que vale 2 en Euler. Sin solución conocida, las diferencias entre aproximaciones sucesivas permiten estimar el orden observado,

Con los valores del cuadro 24.1 para $h=0{,}25$, $0{,}125$ y $0{,}0625$, las diferencias valen $0{,}017954$ y $0{,}008722$, lo que da $p_{\mathrm{obs}}\approx1{,}04$.

Un orden observado correcto no prueba que el código sea correcto. Dos implementaciones equivocadas pueden mostrar patrones de convergencia parecidos. Por eso el refinamiento se combina con casos de referencia, invariantes, límites conocidos y pruebas unitarias, siguiendo la estrategia de verificación iniciada en los capítulos 21 y 22. El refinamiento tampoco valida el modelo frente a datos reales. Solo estudia si el cálculo converge hacia la solución de las ecuaciones que decidimos resolver.

Qué significa afirmar que Euler «aproxima la solución»

La frase supone varias condiciones que suelen quedar implícitas. Debe existir un problema continuo suficientemente regular y la implementación debe reproducir correctamente el algoritmo. El paso tiene que pertenecer a un régimen compatible con la estabilidad y la exactitud. La comparación, por último, debe hacerse sobre un horizonte y con una norma definidos.

La convergencia es una afirmación sobre un límite. Bajo hipótesis apropiadas, cuando $h\to0$ con el horizonte fijo, la trayectoria numérica se acerca a la solución exacta. La afirmación nada dice de que una discretización finita concreta sea una versión físicamente fiel del proceso.

La diferencia tiene alcance filosófico, porque el objeto numérico y el objeto matemático no son idénticos. La ecuación diferencial define una dinámica continua idealizada, mientras que Euler define otra dinámica, discreta y dependiente de $h$, destinada a aproximarla. Podemos demostrar relaciones entre ambas. Lo que no debemos hacer es borrar la diferencia después de obtener una gráfica convincente.

Este problema también debe separarse de la relación entre el modelo y la realidad. Aunque la ecuación diferencial se resolviera con un error numérico despreciable, podría seguir describiendo mal el fenómeno empírico. Reducir el error del integrador no corrige parámetros mal estimados ni mecanismos equivocados.

Una simulación de un modelo continuo combina, por tanto, al menos tres fuentes de discrepancia. Las ecuaciones pueden describir mal el fenómeno, lo que llamamos error de modelamiento. Los parámetros estimados a partir de datos agregan un error paramétrico. Aproximar la solución introduce, finalmente, el error numérico.

Confundir estas fuentes conduce a diagnósticos equivocados. Si una solución numérica oscila porque el paso es demasiado grande, cambiar la teoría del fenómeno sería una mala respuesta. Si la integración converge pero no reproduce observaciones relevantes, seguir reduciendo $h$ tampoco resolverá el problema.

Consolidación

El método de Euler explícito parte de una ecuación $x'=f(t,x)$ y reemplaza la evolución continua durante un intervalo $h$ por una actualización basada en la pendiente al comienzo del intervalo, $x_{n+1}=x_n+h f(t_n,x_n)$. Con suficiente regularidad, su defecto local es $O(h^2)$ y su error global en un horizonte fijo es $O(h)$, de modo que es un método de primer orden.

En la ecuación de prueba $y'=\lambda y$, cada paso multiplica por $R(h\lambda)=1+h\lambda$. La condición $|R|<1$ define la región de estabilidad absoluta. Para $y'=-y$, Euler puede ser estable y producir, sin embargo, signos alternantes que la solución exacta no tiene. Ser estable no asegura la precisión ni la positividad de la aproximación.

El modelo discreto de dos compartimentos coincide algebraicamente con un paso de Euler del modelo continuo de intercambio cuando $\alpha=ha$ y $\beta=hb$. La equivalencia depende de la interpretación. Refinar $h$ con $a$ y $b$ fijos exige cambiar $\alpha$ y $\beta$, mientras que mantenerlos fijos define otra familia de modelos. Euler conserva además exactamente el total, lo que muestra que un esquema numérico puede preservar algunos invariantes y violar otras propiedades del flujo.

El capítulo 25 comparará esta aproximación de primer orden con los métodos de Runge-Kutta y con la selección adaptativa del paso. La pregunta ya no será solo si el resultado converge al reducir $h$. También importará cuánto trabajo computacional exige alcanzar una precisión especificada.

Ejercicios

Ejercicio 24.1. Aplique el método de Euler a $x'=-x$ con $x(0)=1$ hasta $t=3$, con pasos $h=0{,}5$ y $h=1{,}5$. (a) Compare los valores en $t=1{,}5$ y $t=3$ con la solución exacta. (b) ¿Cuál de los dos pasos conserva el signo y la monotonía de la solución? (c) ¿Para qué pasos conserva Euler la positividad estricta en este problema?

Ejercicio 24.2. Para $x'=-x$ en $[0,2]$ con $x(0)=1$, (a) identifique una constante de Lipschitz $L$ y una cota $M$ de $|x''|$. (b) Calcule con $h=0{,}25$ las dos cotas del error global que da el texto y compárelas con el error real del cuadro 24.1. (c) ¿Por qué las cotas son tan conservadoras en este problema?

Ejercicio 24.3. En el modelo continuo de intercambio con $a=0{,}2$ y $b=0{,}1$ por día, (a) calcule el factor exacto por día de la desviación respecto del equilibrio. (b) Calcule el factor de Euler para un día con $h$ igual a 1; 0,5; 0,25 y 0,1, junto con su error. (c) ¿Para qué pasos es no negativa la matriz $I+hA$? ¿Para cuáles es estable el modo no conservado? (d) Describa qué ocurre con $h=6$ días desde el estado $(100,0)$.

Ejercicio 24.4. (a) Determine para qué pasos es estable Euler explícito con $\lambda=-4$ y para cuáles conserva la no negatividad. (b) Repita el análisis de estabilidad con $\lambda=-1+i$. (c) Muestre que con $\lambda=i$ el método es inestable para todo $h>0$.

Ejercicio 24.5. Con los valores del cuadro 24.1 para $h=0{,}25$, $0{,}125$ y $0{,}0625$, (a) estime el orden observado sin usar la solución exacta. (b) Calcule la extrapolación $2x_{h/2}-x_h$ con $h=0{,}125$ y compare su error con el de la aproximación con $h=0{,}0625$. (c) ¿Por qué funciona esa extrapolación?

Ejercicio 24.6. (a) Suponga que la validación de pasos se escribiera como pasos < 0 or int(pasos) != pasos. ¿Qué ocurriría con pasos=30.0? (b) Compare el último tiempo al calcular los tiempos con diez sumas sucesivas de 0,1 y con t0 + h * np.arange(11). (c) ¿Qué comprobación de la función detecta un campo f que devuelve un vector de longitud equivocada?

Referencias

Ascher, U. M. y Petzold, L. R. (1998). Computer methods for ordinary differential equations and differential-algebraic equations. Society for Industrial and Applied Mathematics.

Butcher, J. C. (2016). Numerical methods for ordinary differential equations (3.ª ed.). Wiley. https://doi.org/10.1002/9781119121534

Hairer, E., Nørsett, S. P. y Wanner, G. (1993). Solving ordinary differential equations I: Nonstiff problems (2.ª ed. rev.). Springer.

Higham, N. J. (2002). Accuracy and stability of numerical algorithms (2.ª ed.). Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898718027

Iserles, A. (2008). A first course in the numerical analysis of differential equations (2.ª ed.). Cambridge University Press.