Métodos numéricos para ecuaciones diferenciales ordinarias

Ilustración de la integración numérica para la ecuación diferencial $y'=y,y(0)=1.$ Azul: el Método de Euler, verde: el método del punto medio, rojo: la solución exacta, $y=e^{t}.$ El tamaño del paso es $h=1.0$

La misma ilustración para $h=0.25$. Se ve que el método del punto medio converge más rápido que el método de Euler

Los métodos numéricos para ecuaciones diferenciales ordinarias son procedimientos utilizados para encontrar aproximaciones numericas a las soluciones de ecuaciones diferenciales ordinarias (EDO). Su uso también se conoce como integración numérica, aunque este término a veces se toma para significar el cálculo de una integración.

Muchas ecuaciones diferenciales no pueden resolverse usando funciones típicas ("análisis"). Sin embargo, a efectos prácticos, como en ingeniería, una aproximación numérica a la solución suele ser suficiente. Los algoritmos estudiados aquí pueden usarse para calcular tal aproximación. Un método alternativo es utilizar técnicas de cálculo infinitesimal para obtener una expansión en serie de la solución.

Las ecuaciones diferenciales ordinarias se presentan en muchas disciplinas científicas, por ejemplo, en física, química, biología y economía. Además, algunos métodos en ecuaciones diferenciales parciales numéricas convierten una ecuación diferencial parcial en una ecuación diferencial ordinaria, que luego debe resolverse. (solve by: Wellington Castillo)

El problema

Una ecuación diferencial de primer orden es un problema de valor inicial (PVI) de la forma,​

$y'(t)=f(t,y(t)),\qquad y(t_{0})=y_{0},\qquad \qquad (1)$

donde f es una función que asigna [t₀,∞) × R^d a R^d, con la condición inicial y₀ ∈ R^d es un vector dado. Primer orden significa que solo la primera derivada de ("y") aparece en la ecuación, y las derivadas más altas están ausentes.

Sin pérdida de generalidad en los sistemas de orden superior, en este artículo se restringe la explicación a las ecuaciones diferenciales de "primer orden", porque un EDO de orden superior se puede convertir en un sistema más grande de ecuaciones de primer orden mediante la introducción de variables adicionales. Por ejemplo, la ecuación de segundo orden y'' = -y puede reescribirse como dos ecuaciones de primer orden: y' = z y z' = -y.

En esta sección, se describen métodos numéricos para los PVI, teniendo en cuenta que los problemas de condición de frontera (PCF) requieren un conjunto diferente de herramientas. En un PCF, se definen valores o componentes de la solución y en más de un punto. Debido a esto, se deben usar diferentes métodos para resolverlo. Por ejemplo, el método de disparo (y sus variantes) o métodos globales como las diferencias finitas, el método de Galerkin o el método de colocación son apropiados para esa clase de problemas.

El teorema de Picard-Lindelöf establece que existe una solución única, siempre que f sea lipschitzianamente continua.

Métodos

Los métodos numéricos para resolver PVI de primer orden a menudo se dividen en una de estas dos grandes categorías: método lineal multipaso o método de Runge-Kutta. Se puede lograr una separación adicional dividiendo los métodos en aquellos que son explícitos y aquellos que son implícitos. Por ejemplo, los métodos lineales multipaso implícitos incluyen el método de Adams-Moulton y la fórmula de diferenciación hacia atrás (FDA), mientras que el método de Runge-Kutta​ incluye Runge-Kutta diagonalmente implícito (RKDI), Runge-Kutta diagonalmente implícito simple (RKDIS) y Gauss-Radau (basado en la cuadratura gaussiana). Los ejemplos explícitos de la familia lineal multipaso incluyen el método lineal multipaso, y cualquier método de Runge-Kutta con una diagonal inferior es explícito. Una regla general suelta dicta que las ecuaciones diferenciales rígidas requieren el uso de esquemas implícitos, mientras que los problemas no rígidos se pueden resolver de manera más eficiente con esquemas explícitos.

Los llamados métodos lineales generales (MLG) son una generalización de las dos grandes clases de métodos anteriores.

Método de Euler

Desde cualquier punto de una curva, se puede encontrar una aproximación de otro punto cercano en la curva moviéndose una corta distancia sobre una línea tangente a la curva.

Comenzando con la ecuación diferencial (1), se reemplaza la derivada y' por la aproximación respecto a una diferencia finita

$y'(t)\approx {\frac {y(t+h)-y(t)}{h}},\qquad \qquad (2)$

que cuando se reorganiza produce la siguiente fórmula

$y(t+h)\approx y(t)+hy'(t)\qquad \qquad$

y usando (1) da:

$y(t+h)\approx y(t)+hf(t,y(t)).\qquad \qquad (3)$

Esta fórmula generalmente se aplica de la manera que se explica a continuación.

Se elige el tamaño de paso h y se construye la secuencia t₀, t₁ = t₀ + h, t₂ = t₀ + 2h, ... Denotando por y_n una estimación numérica de la solución exacta y(t_n).

De acuerdo con (3), se calculan estas estimaciones mediante el siguiente esquema recursivo:$y_{n+1}=y_{n}+hf(t_{n},y_{n}).\qquad \qquad (4)$ Este es el método de Euler (en contraste con el método de Euler hacia atrás, que se describe a continuación). El método lleva el nombre de Leonhard Euler que lo describió en 1768.

Es un ejemplo de un método explícito. Esto significa que el nuevo valor y_n+1 se define en términos de datos que ya se conocen, como y_n.

Método de Euler hacia atrás

Si, en lugar de (2), se usa la aproximación

$y'(t)\approx {\frac {y(t)-y(t-h)}{h}},\qquad \qquad (5)$

se obtiene el "método de Euler hacia atrás":

$y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1}).\qquad \qquad (6)$

Esto implica que se trata de un método implicito, lo que significa que previamente se debe resolver otra ecuación con el fin de encontrar y_n+1. Para ello, a menudo se usa el método del punto fijo o alguna modificación del método de Newton-Raphson.

Sin embargo, normalmente cuesta más tiempo resolver esta ecuación que los cálculos de los métodos explícitos. Este costo debe tenerse en cuenta cuando se selecciona el método a utilizar. La ventaja de los métodos implícitos como (6) es que generalmente son más estables para resolver una ecuación rígida, lo que significa que se puede usar un tamaño de paso h más grande.

Método integrador exponencial de primer orden

Los integradores exponenciales describen una gran clase de métodos que han experimentado un gran desarrollo.​ Si origen se remonta al menos a la década de 1960.

En lugar de (1), se asume que la ecuación diferencial es cualquiera de la forma

$y'(t)=-A\,y+{\mathcal {N}}(y),\qquad \qquad \qquad (7)$

o se ha linealizado localmente sobre una forma original para producir un término lineal $-Ay$ y un término no lineal ${\mathcal {N}}(y)$.

Los integradores exponenciales se construyen multiplicando (7) por $e^{At}$ e integrando exactamente el resultado sobre un intervalo de tiempo $[t_{n},t_{n+1}=t_{n}+h]$:

$y_{n+1}=e^{-Ah}y_{n}+\int _{0}^{h}e^{-(h-\tau )A}{\mathcal {N}}\left(y\left(t_{n}+\tau \right)\right)\,d\tau .$

Esta ecuación integral es exacta, pero no define la integral.

El integrador exponencial de primer orden se puede determinar manteniendo ${\mathcal {N}}(y(t_{n}+\tau ))$ constante durante todo el intervalo:

$y_{n+1}=e^{-Ah}y_{n}+A^{-1}(1-e^{-Ah}){\mathcal {N}}(y(t_{n}))\ .\qquad \qquad (8)$

Generalizaciones

El método de Euler a menudo no es lo suficientemente exacto. En términos más precisos, solo tiene orden uno (el concepto de orden se explica a continuación). Esto hizo que los matemáticos buscaran métodos de orden superior.

Una posibilidad es usar no solo el valor previamente calculado y_n para determinar y_n+1, sino hacer que la solución dependa de más valores calculados previamente. Esto produce el llamado "método de varios pasos". Quizás el más simple es el método del salto de rana, que es de segundo orden y (más o menos) se basa en dos valores previos cada vez.

Casi todos los métodos prácticos de varios pasos pertenecen a la familia del método lineal multipaso, que tienen la forma

$\alpha _{k}y_{n+k}+\alpha _{k-1}y_{n+k-1}+\cdots +\alpha _{0}y_{n}$

$=h\left[\beta _{k}f(t_{n+k},y_{n+k})+\beta _{k-1}f(t_{n+k-1},y_{n+k-1})+\cdots +\beta _{0}f(t_{n},y_{n})\right].$

Otra posibilidad es usar más puntos en el intervalo [t_n, t_n+1]. Esto lleva a la familia del método de Runge-Kutta, llamada así por Carl Runge y Martin Wilhelm Kutta. Uno de sus métodos de cuarto orden es especialmente popular.

Funciones avanzadas

Una buena implementación de uno de estos métodos para resolver una EDO a menudo implica algo más que elegir una fórmula y un paso de intervalo adecuados.

A menudo es ineficiente usar el mismo tamaño de paso todo el tiempo, por lo que se han desarrollado "métodos de tamaño de paso variable". Por lo general, el tamaño del paso se elige de modo que el error (local) por paso esté por debajo de cierto nivel de tolerancia. Esto significa que los métodos también deben calcular un "indicador de error", una estimación del error local.

Una extensión de esta idea es elegir dinámicamente entre diferentes métodos de diferentes órdenes (esto se llama un "método de orden variable"). Los métodos basados en la extrapolación de Richardson, como el algoritmo de Bulirsch-Stoer, a menudo se utilizan para construir varios métodos de diferentes órdenes.

Otras características deseables incluyen:

• Salida densa: aproximaciones numéricas asequibles en todo el intervalo de integración, y no solo en los puntos t₀, t₁, t₂, ...
• Ubicación del evento: encontrar los momentos en que, por ejemplo, una función particular desaparece. Esto normalmente requiere el uso de un estimador de resolución numérica de ecuaciones no lineales.
• Soporte para computación paralela.
• Cuando se usa para integrar con respecto al tiempo, reversibilidad del tiempo.

Métodos alternativos

Algunos otros métodos no se detallan el presente artículo. Entre estos métodos alternativos, se encuentran:

• Métodos multiderivadas, que utilizan no solo la función f sino también sus derivadas. Esta clase incluye los métodos de Hermite-Obreschkoff y de Fehlberg, así como procedimientos como el método de Parker-Sochacki o el método de Bychkov-Scherbakov, que calculan los coeficientes de una serie de Taylor de la solución y de forma recursiva.
• Métodos para EDO de segundo orden". Anteriormente se dijo que todas las EDO de orden superior pueden transformarse en EDO de primer orden de la forma (1). Si bien esto es cierto, puede que no siempre sea la mejor manera de proceder. En particular, el método de Nyström trabaja directamente con ecuaciones de segundo orden.
• Los métodos de integración geométrica están especialmente diseñados para clases especiales de EDO (por ejemplo, el integrador simpléctico, adecuado para determinar la solución de ecuaciones hamiltonianas), ya que tienen en cuenta que la solución numérica respete la estructura o la geometría subyacente de estas clases.
• Los métodos de sistemas de estado cuantificados forman una familia de procedimientos de integración de EDO basada en la idea de la cuantificación de estado. Son eficientes al simular sistemas dispersos con discontinuidades frecuentes.

Métodos paralelos en el tiempo

Para aplicaciones que requieren computación paralela en supercomputadoras, el grado de complejidad requerido por un método numérico se vuelve relevante. En vista de los desafíos de los sistemas informáticos que trabajan a exaescala, se están estudiando métodos numéricos para problemas de valor inicial que puedan proporcionar concurrencia en la dirección temporal.​ Parareal es un ejemplo relativamente conocido de tal método de integración paralelo en el tiempo, pero las primeras ideas al respecto se remontan a la década de 1960.​

Análisis

El análisis numérico no es solo el diseño de métodos numéricos, sino que también implica el estudio de cómo funcionan. Tres conceptos son centrales en este estudio:

• Convergencia: si el método se aproxima a la solución,
• Orden: la eficiencia con la que se aproxima a la solución, y
• Estabilidad: si los errores permanecen acotados.

Convergencia

Se dice que un método numérico es convergente si la solución numérica se aproxima a la solución exacta cuando el tamaño de paso h tiende a 0. Más precisamente, se requiere que para cada EDO (1) con una función lipschitziana f y cada t^* > 0,

$\lim _{h\to 0+}\max _{n=0,1,\dots ,\lfloor t^{*}/h\rfloor }\|y_{n,h}-y(t_{n})\|=0.$

Todos los métodos mencionados anteriormente son convergentes.

Consistencia y orden

Supóngase que el método numérico es

$y_{n+k}=\Psi (t_{n+k};y_{n},y_{n+1},\dots ,y_{n+k-1};h).\,$

El "error local (de truncamiento)" del método es el error cometido en un paso del método. Es decir, es la diferencia entre el resultado dado por el método, suponiendo que no se cometió ningún error en los pasos anteriores, y la solución exacta:

$\delta _{n+k}^{h}=\Psi \left(t_{n+k};y(t_{n}),y(t_{n+1}),\dots ,y(t_{n+k-1});h\right)-y(t_{n+k}).$

Se dice que el método es consistente si

$\lim _{h\to 0}{\frac {\delta _{n+k}^{h}}{h}}=0.$

El método tiene orden $p$ si

$\delta _{n+k}^{h}=O(h^{p+1})\quad {\mbox{as }}h\to 0.$

Por lo tanto, un método es consistente si tiene un orden mayor que 0. El método de Euler (hacia adelante) (4) y el método de Euler hacia atrás (6) introducido anteriormente tienen orden 1, por lo que son consistentes. La mayoría de los métodos utilizados en la práctica alcanzan un orden superior. La consistencia es una condición necesaria para la convergencia ^{[cita requerida]}, pero no suficiente; Para que un método sea convergente, debe ser consistente y cero-estable.

Un concepto relacionado es el de "error global (de truncamiento)", el error acumulado en todos los pasos necesarios para alcanzar un tiempo fijo t. Explícitamente, el error global en el tiempo t es y_N − y(t) donde N = (t−t₀)/h. El error global de un método de un solo paso p es O (h^p); en particular, tal método es convergente. Esta afirmación no es necesariamente cierta para los métodos de varios pasos.

Estabilidad y rigidez

Para algunas ecuaciones diferenciales, la aplicación de métodos estándar, como el método de Euler, el método de Runge-Kutta explícito o el método lineal multipaso (por ejemplo, el método de Adams-Bashforth), exhiben inestabilidad en las soluciones, aunque otros métodos pueden producir soluciones estables. Este "comportamiento difícil" en la ecuación (que puede no ser necesariamente complejo en sí mismo) se describe como "rigidez", y a menudo es causado por la presencia de diferentes escalas de tiempo en el problema subyacente. Por ejemplo, una colisión en un sistema mecánico como en un oscilador de impacto generalmente ocurre en una escala de tiempo mucho más pequeña que la necesaria para el caracterizar el movimiento de los objetos. Esta discrepancia genera "virajes bruscos" en las curvas de los parámetros de estado.

Los problemas rígidos son ubicuos en cinética química, teoría del control, mecánica de sólidos deformables, pronóstico del tiempo, biología, plasma y electrónica. Una forma de superar la rigidez es extender la noción de ecuación diferencial a la de inclusión diferencial, que permite y modela la falta de suavidad.​​

Historia

A continuación se muestra una cronología de algunos desarrollos importantes en este campo:

• 1768 - Leonhard Euler publica su método.
• 1824 - Augustin Louis Cauchy prueba la convergencia del método de Euler, para lo que usa el método implícito de Euler.
• 1855 - Primera mención del método lineal multipaso de John Couch Adams en una carta escrita por Francis Bashforth.
• 1895 - Carl Runge publica el primer método de Runge-Kutta.
• 1901 - Martin Wilhelm Kutta describe el popular método de Runge-Kutta de cuarto orden.
• 1910 - Lewis Fry Richardson anuncia el método de extrapolación que lleva su nombre, el denominado método de extrapolación de Richardson.
• 1952 - Charles F. Curtiss y Joseph Oakland Hirschfelder acuñan el término ecuación rígida.
• 1963 - Germund Dahlquist introduce el concepto de A-estabilidad de los métodos de integración.

Soluciones numéricas a problemas de valor límite unidimensional de segundo orden

Los problemas de valor límite (PVL) generalmente se resuelven numéricamente solucionando un problema de matriz aproximadamente equivalente, obtenido al discretizar el PVL original. El procedimiento más utilizado para resolver numéricamente estos problemas en una dimensión se llama método de las diferencias finitas. Este método aprovecha las combinaciones lineales de valores de puntos para construir coeficientes de diferencias finitas que describen derivadas de la función. Por ejemplo, la aproximación en diferencia finita de segundo orden a la primera derivada viene dada por:

${\frac {u_{i+1}-u_{i-1}}{2h}}=u'(x_{i})+{\mathcal {O}}(h^{2}),$

y la diferencia finita de segundo orden para la segunda derivada viene dado por:

${\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}=u''(x_{i})+{\mathcal {O}}(h^{2}).$

En ambas fórmulas, $h=x_{i}-x_{i-1}$ es la distancia entre los valores vecinos x en el dominio discretizado. Luego se construye un sistema lineal que luego se puede resolver con el estándar matrix methods. Por ejemplo, supongamos que la ecuación a resolver es:

${\frac {d^{2}u}{dx^{2}}}-u=0,$

$u(0)=0,$

$u(1)=1.$

El siguiente paso sería discretizar el problema y usar aproximaciones lineales derivadas como

$u''_{i}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}$

y resolver el sistema resultante de ecuaciones lineales. Esto llevaría a ecuaciones como:

${\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}-u_{i}=0,\quad \forall i={1,2,3,...,n-1}.$

A primera vista, este sistema de ecuaciones parece tener dificultades asociadas con el hecho de que la ecuación no implica términos que no se multipliquen por variables, pero en realidad esto es falso. En i=1 y nminus&;1 existe un término que involucra los valores límite $u(0)=u_{0}$ y $u(1)=u_{n}$ y dado que estos dos valores son conocidos, simplemente se puede sustituirlos en esta ecuación y, como resultado, obtener un sistema lineal de ecuaciones no homogéneo que tiene soluciones no triviales.

Véase también

• Número de Courant-Friedrichs-Levy
• Deriva energética
• Métodos lineales generales
• Anexo:Lista de temas de análisis numérico
• Algoritmo de propagación reversible del sistema de referencia
• Lenguaje Modelica y software OpenModelica

Bibliografía

• J. Arrieta, R. Ferreira, R. Pardo y A. Rodríguez-Bernal. "Análisis Numérico de Ecuaciones Diferenciales Ordinarias". Paraninfo, Madrid, 2020. ISBN 9788428344418, ISBN 8428344418.
• Bradie, Brian (2006). A Friendly Introduction to Numerical Analysis. Upper Saddle River, New Jersey: Pearson Prentice Hall. ISBN 978-0-13-013054-9. 
• J. C. Butcher, Métodos numéricos para ecuaciones diferenciales ordinarias , ISBN 0-471-96758-0
• Ernst Hairer, Syvert Paul Nørsett y Gerhard Wanner, "Resolviendo ecuaciones diferenciales ordinarias I: problemas no rígidos", segunda edición, Springer Verlag, Berlín, 1993. ISBN 3-540-56670-8.
• Ernst Hairer y Gerhard Wanner, "Resolviendo ecuaciones diferenciales ordinarias II: problemas rígidos y diferenciales-algebraicos", segunda edición, Springer Verlag, Berlín, 1996. ISBN 3-540-60452-9.
  (This two-volume monograph systematically covers all aspects of the field.)
• Hochbruck, Marlis; Ostermann, Alexander (May 2010). «Exponential integrators». Acta Numerica 19: 209-286. Bibcode:2010AcNum..19..209H. doi:10.1017/S0962492910000048. «citeseerx 10.1.1.187.6794». 
• Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996. ISBN 0-521-55376-8 (hardback), ISBN 0-521-55655-4 (paperback).
  (Textbook, targeting advanced undergraduate and postgraduate students in mathematics, which also discusses numerical partial differential equations.)
• John Denholm Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, Chichester, 1991. ISBN 0-471-92990-5.
  (Textbook, slightly more demanding than the book by Iserles.)

Enlaces externos

• Joseph W. Rudmin, Application of the Parker–Sochacki Method to Celestial Mechanics, 1998.
• Dominique Tournès, L'intégration approchée des équations différentielles ordinaires (1671-1914), thèse de doctorat de l'université Paris 7 - Denis Diderot, juin 1996. Réimp. Villeneuve d'Ascq : Presses universitaires du Septentrion, 1997, 468 p. (Extensive online material on ODE numerical analysis history, for English-language material on the history of ODE numerical analysis, see e.g. the paper books by Chabert and Goldstine quoted by him.)
• kv en GitHub (Librería C++ con resolución rigurosa de EDO)
• INTLAB Archivado el 25 de noviembre de 2019 en Wayback Machine. (Librería MATLAB/GNU Octave con resolución rigurosa de EDO)

• Datos: Q541182