4: Ecuaciones diferenciales

Ahora discutimos la solución numérica de las ecuaciones diferenciales ordinarias. Incluiremos el problema del valor inicial y el problema del valor límite.

Problema de valor inicial

Comenzamos con el método simple de Euler, luego discutimos los métodos RungeKutta más sofisticados y concluimos con el método Runge-Kutta-Fehlberg, tal como se implementa en la función de MATLAB ode\(45 . \mathrm{m}\). Nuestras ecuaciones diferenciales son para\(x=x(t)\), donde el tiempo\(t\) es la variable independiente.

Método Euler

El método de Euler es el método más sencillo para integrar una ecuación diferencial. Considerar la ecuación diferencial de primer orden

\[\dot{x}=f(t, x), \nonumber \]

con la condición inicial\(x(0)=x_{0}\). Definir\(t_{n}=n \Delta t\) y\(x_{n}=x\left(t_{n}\right)\). Una expansión de los\(x_{n+1}\) resultados de la serie Taylor

\[\begin{aligned} x_{n+1} &=x\left(t_{n}+\Delta t\right) \\ &=x\left(t_{n}\right)+\Delta t \dot{x}\left(t_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) \\ &=x\left(t_{n}\right)+\Delta t f\left(t_{n}, x_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) . \end{aligned} \nonumber \]

Por lo tanto, el Método Euler está escrito como

\[x_{n+1}=x\left(t_{n}\right)+\Delta t f\left(t_{n}, x_{n}\right) . \nonumber \]

Decimos que el método de Euler avanza en el tiempo usando un paso de tiempo\(\Delta t\), a partir del valor inicial\(x_{0}=x(0)\). El error local del Método Euler es\(\mathrm{O}\left(\Delta t^{2}\right)\). El error global, sin embargo, incurrido al integrarse a un tiempo\(T\), es un factor de\(1 / \Delta t\) mayor tamaño y viene dado por\(\mathrm{O}(\Delta t)\). Por lo tanto, es costumbre llamar al Método Euler un método de primer orden.

Método Euler modificado

Este método es un método llamado predictor-corrector. También es el primero de lo que veremos son los métodos Runge-Kutta. Como antes, queremos resolver (4.1). La idea es promediar el valor de\(\dot{x}\) al inicio y al final del paso de tiempo. Es decir, nos gustaría modificar el método de Euler y escribir

\[x_{n+1}=x_{n}+\frac{1}{2} \Delta t\left(f\left(t_{n}, x_{n}\right)+f\left(t_{n}+\Delta t, x_{n+1}\right)\right) . \nonumber \]

El problema obvio con esta fórmula es que el valor desconocido\(x_{n+1}\) aparece en el lado derecho. Podemos, sin embargo, estimar este valor, en lo que se llama el paso predictor. Para el paso predictor, utilizamos el método de Euler para encontrar

\[x_{n+1}^{p}=x_{n}+\Delta t f\left(t_{n}, x_{n}\right) \text {. } \nonumber \]

El paso corrector se convierte entonces

\[x_{n+1}=x_{n}+\frac{1}{2} \Delta t\left(f\left(t_{n}, x_{n}\right)+f\left(t_{n}+\Delta t, x_{n+1}^{p}\right)\right) . \nonumber \]

El Método Euler Modificado se puede reescribir de la siguiente forma que posteriormente identificaremos como un método Runge-Kutta

\[\begin{align}\nonumber &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{1}\right), \\ &x_{n+1}=x_{n}+\frac{1}{2}\left(k_{1}+k_{2}\right) . \nonumber \end{align} \nonumber \]

Métodos Runge-Kutta de segundo orden

Ahora derivamos la familia completa de métodos Runge-Kutta de segundo orden. Los métodos de orden superior pueden derivarse de manera similar, pero requieren sustancialmente más álgebra.

De nuevo consideramos la ecuación diferencial dada por (4.1). Un método general Runge-Kutta de segundo orden puede escribirse en la forma

\[\begin{align} \nonumber &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta k_{1}\right), \\ &x_{n+1}=x_{n}+a k_{1}+b k_{2}, \nonumber \end{align} \nonumber \]

con\(\alpha, \beta, a\) y\(b\) constantes que definen el método Runge-Kutta de segundo orden particular. Estas constantes deben restringirse estableciendo el error local del método Runge-Kutta de segundo orden para ser\(\mathrm{O}\left(\Delta t^{3}\right)\). Intuitivamente, podríamos adivinar que dos de las limitaciones serán\(a+b=1\) y\(\alpha=\beta\).

Calculamos la serie Taylor\(x_{n+1}\) directamente, y a partir del método Runge-Kutta, y requerimos que sean los mismos a pedido\(\Delta t^{2}\). Primero, calculamos la serie Taylor de\(x_{n+1}\). Tenemos

\[\begin{aligned} x_{n+1} &=x\left(t_{n}+\Delta t\right) \\ &=x\left(t_{n}\right)+\Delta t \dot{x}\left(t_{n}\right)+\frac{1}{2}(\Delta t)^{2} \ddot{x}\left(t_{n}\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{aligned} \nonumber \]

Ahora,

\[\dot{x}\left(t_{n}\right)=f\left(t_{n}, x_{n}\right) . \nonumber \]

La segunda derivada es más complicada y requiere derivadas parciales. Tenemos

\[\begin{aligned} \ddot{x}\left(t_{n}\right) &\left.=\frac{d}{d t} f(t, x(t))\right]_{t=t_{n}} \\ &=f_{t}\left(t_{n}, x_{n}\right)+\dot{x}\left(t_{n}\right) f_{x}\left(t_{n}, x_{n}\right) \\ &=f_{t}\left(t_{n}, x_{n}\right)+f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right) . \end{aligned} \nonumber \]

Por lo tanto,

\[\begin{align} \nonumber x_{n+1}=x_{n}+\Delta t f\left(t_{n}, x_{n}\right) & \\ &+\frac{1}{2}(\Delta t)^{2}\left(f_{t}\left(t_{n}, x_{n}\right)+f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{align} \nonumber \]

Segundo, calculamos\(x_{n+1}\) a partir del método Runge-Kutta dado por (4.7). Combinando (4.7) en una sola expresión, tenemos

\[x_{n+1}=x_{n}+a \Delta t f\left(t_{n}, x_{n}\right) \nonumber \]

\[+b \Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta \Delta t f\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \nonumber \]

Nosotros la serie Taylor expandimos usando

\[\begin{aligned} f\left(t_{n}+\alpha \Delta t, x_{n}+\beta \Delta t f\left(t_{n}, x_{n}\right)\right) & \\ &=f\left(t_{n}, x_{n}\right)+\alpha \Delta t f_{t}\left(t_{n}, x_{n}\right)+\beta \Delta t f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) . \end{aligned} \nonumber \]

Por lo tanto, la fórmula Runge-Kutta es

\[\begin{aligned} x_{n+1}=x_{n}+(a+b) \Delta t f &\left(t_{n}, x_{n}\right) \\ &+(\Delta t)^{2}\left(\alpha b f_{t}\left(t_{n}, x_{n}\right)+\beta b f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{aligned} \nonumber \]

Comparando (4.9) y (4.10), encontramos

\[a+b=1, \quad \alpha b=1 / 2, \quad \beta b=1 / 2 . \nonumber \]

Dado que solo hay tres ecuaciones para cuatro parámetros, existe una familia de métodos Runge-Kutta de segundo orden.

El método de Euler modificado dado por (4.6) corresponde a\(\alpha=\beta=1\) y\(a=\)\(b=1 / 2\). Otro método Runge-Kutta de segundo orden, llamado método de punto medio, corresponde a\(\alpha=\beta=1 / 2, a=0\) y\(b=1\). Este método está escrito como

\[\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}\right), \\ &x_{n+1}=x_{n}+k_{2} . \end{aligned} \nonumber \]

Métodos Runge-Kutta de orden superior

El método general de Runge-Kutta de segundo orden fue dado por (4.7). La forma general del método de tercer orden viene dada por

\[\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta k_{1}\right), \\ &k_{3}=\Delta t f\left(t_{n}+\gamma \Delta t, x_{n}+\delta k_{1}+\epsilon k_{2}\right), \\ &x_{n+1}=x_{n}+a k_{1}+b k_{2}+c k_{3} . \end{aligned} \nonumber \]

Se pueden adivinar las siguientes restricciones sobre las constantes:\(\alpha=\beta, \gamma=\delta+\epsilon\), y\(a+b+c=1\). Es necesario derivar las restricciones restantes.

El método de cuarto orden tiene un\(k_{1}, k_{2}, k_{3}\) y\(k_{4}\). El método de quinto orden requiere al menos para\(k_{6}\). La siguiente tabla da el orden mínimo del método y el número de etapas requeridas.

ordenar	2	3	4	5	6	7	8
# etapas mínimas	2	3	4	6	7	9	11

Debido al salto en el número de etapas requeridas entre los métodos de cuarto y quinto orden, el método Runge-Kutta de cuarto orden tiene cierto atractivo. El método general de cuarto orden con cuatro etapas tiene 13 constantes y 11 restricciones. Un método particularmente simple de cuarto orden que ha sido ampliamente utilizado en el pasado por los físicos es dado por

\[\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right) \\ &k_{2}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}\right) \\ &k_{3}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}\right) \\ &k_{4}=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{3}\right) \\ &x_{n+1}=x_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right) . \end{aligned} \nonumber \]

Métodos adaptativos Runge-Kutta

Al igual que en la integración adaptativa, es útil idear un integrador oda que encuentre automáticamente lo apropiado\(\Delta t\). El método del príncipe latente, que se implementa en la oda de MATLAB\(45 . \mathrm{m}\), encuentra el tamaño de paso apropiado comparando los resultados de un método de quinto y cuarto orden. Requiere seis evaluaciones de función por paso de tiempo y construye un método de quinto y cuarto orden a partir de las mismas evaluaciones de función.

Supongamos que el método de quinto orden encuentra\(x_{n+1}\) con error local\(\mathrm{O}\left(\Delta t^{6}\right)\), y el método de cuarto orden encuentra\(x_{n+1}^{\prime}\) con error local\(\mathrm{O}\left(\Delta t^{5}\right)\). Dejar\(\varepsilon\) ser la tolerancia de error deseada del método, y dejar\(e\) ser el error real. Podemos estimar\(e\) a partir de la diferencia entre los métodos de quinto y cuarto orden; es decir,

\[e=\left|x_{n+1}-x_{n+1}^{\prime}\right| \text {. } \nonumber \]

Ahora\(e\) es de\(\mathrm{O}\left(\Delta t^{5}\right)\), donde\(\Delta t\) se toma el tamaño del paso. \(\Delta \tau\)Sea el tamaño estimado de paso requerido para obtener el error deseado\(\varepsilon\). Entonces tenemos

\[e / \varepsilon=(\Delta t)^{5} /(\Delta \tau)^{5}, \nonumber \]

o resolviendo para\(\Delta \tau\),

\[\Delta \tau=\Delta t\left(\frac{\varepsilon}{e}\right)^{1 / 5} \nonumber \]

Por un lado, si el error real es menor que el error deseado,\(e<\varepsilon\), entonces aceptamos\(x_{n+1}\) y hacemos el siguiente paso de tiempo usando el valor mayor de\(\Delta \tau\). Por otro lado, si el error real es mayor que el error deseado\(e>\varepsilon\), entonces rechazamos el paso de integración y rehacemos el paso de tiempo usando el valor menor de\(\Delta \tau\). En la práctica, uno suele aumentar el paso de tiempo un poco menos y disminuye el paso de tiempo un poco más para evitar el desperdicio de demasiados pasos de tiempo fallidos.

Sistema de ecuaciones diferenciales

Nuestros métodos numéricos se pueden adaptar fácilmente para resolver ecuaciones diferenciales de orden superior, o equivalentemente, un sistema de ecuaciones diferenciales. Primero, mostramos cómo una ecuación diferencial de segundo orden se puede reducir a dos ecuaciones de primer orden. Considerar

\[\ddot{x}=f(t, x, \dot{x}) . \nonumber \]

Esta ecuación de segundo orden se puede reescribir como dos ecuaciones de primer orden definiendo\(u=\dot{x}\). Entonces tenemos el sistema

\[\dot{x}=u, \quad \dot{u}=f(t, x, u) . \nonumber \]

Este truco también funciona para ecuaciones de orden superior. Por ejemplo, la ecuación de tercer orden

\[\dddot x=f(t, x, \dot{x}, \ddot{x}), \nonumber \]

se puede escribir como

\[\dot{x}=u, \quad \dot{u}=v, \quad \dot{v}=f(t, x, u, v) . \nonumber \]

Podemos generalizar el método Runge-Kutta para resolver un sistema de ecuaciones diferenciales. A modo de ejemplo, consideremos el siguiente sistema de dos odas,

\[\dot{x}=f(t, x, y), \quad \dot{y}=g(t, x, y), \nonumber \]

con las condiciones iniciales\(x(0)=x_{0}\) y\(y(0)=y_{0}\). La generalización del método Runge-Kutta de cuarto orden comúnmente utilizado sería

\[\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}, y_{n}\right), \\ l_{1} &=\Delta t g\left(t_{n}, x_{n}, y_{n}\right), \\ k_{2} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}, y_{n}+\frac{1}{2} l_{1}\right), \\ l_{2} &=\Delta t g\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}, y_{n}+\frac{1}{2} l_{1}\right), \\ k_{3} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}, y_{n}+\frac{1}{2} l_{2}\right), \\ l_{3} &=\Delta t g\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}, y_{n}+\frac{1}{2} l_{2}\right), \\ k_{4} &=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{3}, y_{n}+l_{3}\right), \\ l_{4} &=\Delta t g\left(t_{n}+\Delta t, x_{n}+k_{3}, y_{n}+l_{3}\right), \\ x_{n+1} &=x_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right), \\ y_{n+1} &=y_{n}+\frac{1}{6}\left(l_{1}+2 l_{2}+2 l_{3}+l_{4}\right) . \end{aligned} \nonumber \]

Problemas de valor límite

Método de disparo

Consideramos la oda general de la forma

\[\frac{d^{2} y}{d x^{2}}=f(x, y, d y / d x), \nonumber \]

con condiciones de límite de dos puntos\(y(0)=A\) y\(y(1)=B\). Primero formularemos la oda como un problema de valor inicial. Tenemos

\[\frac{d y}{d x}=z, \quad \frac{d z}{d x}=f(x, y, z) . \nonumber \]

Se conoce la condición inicial, pero\(y(0)=A\) se\(z(0)=b\) desconoce la segunda condición inicial. Nuestro objetivo es determinar\(b\) tal que\(y(1)=B\).

De hecho, este es un problema de búsqueda de raíces para una función adecuadamente definida. Definimos la función de\(F=F(b)\) tal manera que

\[F(b)=y(1)-B, \nonumber \]

donde\(y(1)\) es el valor numérico obtenido de integrar las ecuaciones diferenciales acopladas de primer orden con\(y(0)=A\) y\(z(0)=b\). Nuestra rutina de búsqueda de raíces querrá resolver\(F(b)=0\). (El método se llama disparo porque la pendiente de la curva de solución para\(y=y(x)\) at\(x=0\) viene dada por\(b\), y la solución alcanza el valor\(y(1)\) at\(x=1\). Esto parece apuntar con un arma y tratar de disparar al objetivo, que es\(B\).)

Para determinar el valor de\(b\) que resuelve\(F(b)=0\), iteramos usando el método secante, dado por

\[b_{n+1}=b_{n}-F\left(b_{n}\right) \frac{b_{n}-b_{n-1}}{F\left(b_{n}\right)-F\left(b_{n-1}\right)} . \nonumber \]

Tenemos que comenzar con dos conjeturas iniciales para\(b\), resolviendo la oda para los dos valores correspondientes de\(y(1)\). Entonces el método secante nos dará el siguiente valor de\(b\) para probar, e iteramos hasta\(|y(1)-B|<\) tol, donde tol es alguna tolerancia especificada para el error.