5.7: Enfoques numéricos

Si la amplitud de las oscilaciones, por cualquier razón, llega a ser tan grande que los términos no lineales en la ecuación que describe un oscilador se vuelven comparables con sus términos lineales, los métodos numéricos son prácticamente la única vía disponible para sus estudios teóricos. En los sistemas Hamiltonianos 1D, dichos métodos pueden aplicarse directamente a la Ec. (3.26), pero los sistemas disipativos y/o paramétricos suelen carecer de tales primeras integrales de movimiento, por lo que la ecuación diferencial inicial tiene que ser resuelta.

Discutamos la idea general de tales métodos sobre el ejemplo de lo que los matemáticos llaman el problema de Cauchy (encontrar la solución para todos los momentos del tiempo, a partir de las condiciones iniciales conocidas) para la ecuación diferencial de primer orden\[\dot{q}=f(t, q) \text {. }\] (La generalización a un sistema de varias ecuaciones de este tipo es sencillo.) Rompiendo el eje de tiempo en pequeños pasos iguales\(h\) (Figura 11) podemos reducir el problema de integración de ecuaciones a encontrar el valor de la función en el siguiente punto temporal,\(q_{n+1} \equiv q\left(t_{n+1}\right) \equiv q\left(t_{n}+h\right)\) a partir del valor previamente encontrado\(q_{n}=q\left(t_{n}\right)-\) y, si es necesario, los valores de\(q\) otros pasos de tiempo anteriores.

Figura 5.11. Las nociones básicas utilizadas en la integración numérica de ecuaciones diferenciales ordinarias.

En el enfoque más simple (llamado el método de Euler),\(q_{n+1}\) se encuentra usando la siguiente fórmula:\[\begin{aligned} &q_{n+1}=q_{n}+k, \\ &k \equiv h f\left(t_{n}, q_{n}\right) . \end{aligned}\] Esta aproximación equivale a la sustitución de la función genuina\(q(t)\), en el segmento\(\left[t_{n}\right.\),\(t_{n+1}\)], con los dos primeros términos de su expansión Taylor en punto \(t_{n}\):\[q\left(t_{n}+h\right) \approx q\left(t_{n}\right)+\dot{q}\left(t_{n}\right) h \equiv q\left(t_{n}\right)+h f\left(t_{n}, q_{n}\right) .\] Esta aproximación tiene un error proporcional a\(h^{2}\). Se podría argumentar que al hacer el paso\(h\) lo suficientemente pequeño, el error del método Euler podría hacerse arbitrariamente pequeño, pero incluso con toda la potencia de cálculo de números de las plataformas informáticas modernas, el tiempo de CPU necesario para alcanzar la precisión suficiente puede ser demasiado grande para grandes problemas. \({ }^{35}\)Además de eso, el incremento del número de pasos de tiempo, que es necesario\(h \rightarrow 0\) en un intervalo de tiempo total fijo, incrementa los errores totales de redondeo y eventualmente puede provocar un aumento, en lugar de la reducción del error general del resultado calculado.

Una forma más eficiente es modificar la Ec. (96) para incluir los términos de segundo orden en\(h\). Hay varias formas de hacerlo, por ejemplo usando el método\(2^{\text {nd }}\) -order Runge-Kutta:\[\begin{aligned} &q_{n+1}=q_{n}+k_{2}, \\ &k_{2} \equiv h f\left(t_{n}+\frac{h}{2}, q_{n}+\frac{k_{1}}{2}\right), \quad k_{1} \equiv h f\left(t_{n}, q_{n}\right) . \end{aligned}\] Uno puede verificar fácilmente que este método da el resultado exacto si la función\(q(t)\) es un polinomio cuadrático, y por lo tanto en el caso general sus errores son del orden de\(h^{3}\). Vemos que la idea principal aquí es romper primero el segmento\(\left[t_{n}, t_{n+1}\right]\) por la mitad (ver Figura 11), evaluar el lado derecho de la ecuación diferencial (95) en el punto intermedio (en ambos\(t\) y\(q\)) entre el número de puntos\(n\) y\((n+1)\), y luego usar esta información para predecir\(q_{n+1}\).

La ventaja del enfoque Runge-Kutta sobre otros métodos de segundo orden es que puede extenderse fácilmente al\(4^{\text {th }}\) orden, sin una ruptura adicional del intervalo\(\left[t_{n}, t_{n+1}\right] .:\)\[\begin{aligned} &q_{n+1}=q_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right), \\ &k_{4} \equiv h f\left(t_{n}+h, q_{n}+k_{3}\right), k_{3} \equiv h f\left(t_{n}+\frac{h}{2}, q_{n}+\frac{k_{2}}{2}\right), k_{2} \equiv h f\left(t_{n}+\frac{h}{2}, q_{n}+\frac{k_{1}}{2}\right), k_{1} \equiv h f\left(t_{n}, q_{n}\right) . \end{aligned}\] Este método tiene un error mucho menor\(O\left(h^{5}\right)\),, sin ser demasiado engorroso. Estas características han hecho del\(4^{\text {th }}\) orden -order Runge-Kutta el método predeterminado en la mayoría de las bibliotecas numéricas. Su extensión a órdenes superiores es posible, pero requiere fórmulas más complejas, y se justifica solo para algunos casos especiales, por ejemplo, funciones muy abruptas\(q(t) .{ }^{36}\) La mejora más frecuente del método es un ajuste automático del paso\(h\) para alcanzar la precisión preespecificada, pero no hacer más cálculos de los necesarios.

La Figura 12 muestra un ejemplo típico de una aplicación de ese método al problema muy simple de un oscilador lineal amortiguado, para dos valores de paso de tiempo fijo\(h\) (expresados en términos del número\(N\) de tales pasos por periodo de oscilación). Las líneas rectas negras conectan los puntos adyacentes obtenidos por el método de\(4^{\text {th }}\) orden Runge-Kutta, mientras que los puntos conectados con las líneas rectas verdes representan la solución analítica exacta (22). Las gráficas muestran que algunos errores por ciento comienzan a aparecer solo en tan pocos como pasos de\(\sim 10\) tiempo por periodo, por lo que el método es realmente muy eficiente.

Figura 5.12. Resultados de la solución Runge-Kutta de la Ec. (6) (con\(\delta / \omega_{0}=0.03\)) para: (a) 30 y (b) 6 puntos por periodo de oscilación. Los resultados se muestran por puntos; las líneas negras y verdes son solo las guías para el ojo.

Permítanme esperar que la discusión en la siguiente sección haga muy claras las conveniencias y los handicaps del enfoque numérico a los problemas de la dinámica no lineal.

---

\({ }^{35}\)Además, el método de Euler no es reversible en el tiempo, el hándicap que puede ser esencial para la integración de los sistemas hamiltonianos descritos por sistemas de ecuaciones diferenciales de segundo orden. Sin embargo, este inconveniente puede superarse en parte por el llamado salto, el solapamiento de pasos de tiempo\(h\) para una coordenada generalizada y la velocidad generalizada correspondiente.

\({ }^{36}\)Los enfoques más populares en tales casos son la extrapolación de Richardson, el algoritmo Bulirsch-Stoer y un conjunto de las llamadas técnicas de predicción-corrección, por ejemplo el método Adams-Bashforth-Moulton - ver la literatura recomendada en MA Sec. 16 (iii).