¿Por qué la solución numérica de esta ecuación es inestable? ¿Es esta ecuación rígida?

Question

Estoy tratando de resolver la siguiente ecuación con un cuarto orden explícito utilizando el método Runge-Kutta: $$y' = t(y - t \sin t)$$ con condiciones iniciales $y(0) = 1$ en el intervalo $[0, 10]$ .

La solución analítica es $y = t \sin t + \cos t$ pero el método numérico explota incluso para un paso de tiempo de $0.0001$ .

¿Alguna idea?

Answer 1

Sí, este problema de valor inicial está hecho para hacer esas cosas. A petición de Robert Lewis, aquí hay una imagen: la curva roja es la solución analítica, la curva azul es la implementación de libro de texto de RK4 con los parámetros de su puesto.

(La curva azul es bastante vertical a partir de ahí.) No creo que esta ecuación pueda calificarse como rígido (término que nunca he entendido bien), es simplemente hipersensible a la condición inicial. De hecho, la solución exacta con la condición inicial $y(0)=1+\delta$ es $$y(t) = t\sin t + \cos t + \delta\, e^{t^2/2} \tag{1}$$ donde el último término significa problemas. En cuanto nos salimos de la trayectoria exacta lo más mínimo, el $e^{t^2/2}$ y condena la solución, incluso si los cálculos posteriores no tienen ningún error.

Esto no es realmente un fallo del método. Digamos que el valor inicial $y(0)$ se conoce hasta $2^{-53}$ (en doble precisión). Entonces, ¿qué ordenador sabe realmente $\delta$ en (1) es que $|\delta|\le 2^{-53}$ . Pero $2^{-53} e^{t^2/2} > 500,000$ cuando $t=10$ . Esto es para el analítica solución con valor inicial de precisión de máquina. Así que el problema no está (sólo) en el método RK4, no tenemos suficiente precisión al principio para obtener un resultado decente al final. Realizar los cálculos en precisión cuádruple debería ayudar.

En cierto modo, lo de la trigonometría aquí es más que nada un adorno. La ecuación más simple $y'=ty$ ya tiene la misma característica. Con el valor inicial $0$ la solución es $y\equiv 0$ . Pero con un valor inicial no nulo $C$ es $y= C \exp (t^2/2)$ que es enorme.

Por cierto, he probado con solucionadores de ODE rígidos como ode15s en Matlab; eso no ayudó en absoluto.

---

A modo de apunte: considere lo que ocurre si invertimos el orden de los ejes temporales. Es decir, dejemos que $t = 10-s$ por lo que la función desconocida es $z(s) = y(10-t)$ . Entonces resuelve $$z' = - (10-s) (z-(10-s)\sin(10-s))$$ con alguna condición inicial $z(0)=C$ . Entonces el método numérico funciona bien; no se necesitan pasos muy pequeños. Por ejemplo, a continuación he resuelto con $C = -6.1$ utilizando sólo $40$ pasos de RK4 y trazamos la solución a la inversa:

En cierto modo, esto es una trampa: ¿cómo sabía yo que $C=-6.1$ ¿es el valor correcto a utilizar? (Lo tomé del valor del punto final de la solución exacta).

Answer 2

Sé que esta pregunta es vieja, pero pensé en publicar algo de información adicional sobre esto.

Una de las razones por las que la solución numérica de la ecuación diferencial ordinaria (EDO) explota, independientemente del método numérico utilizado, es por la propia EDO. Si usted va y estudia alguna teoría sobre la integración numérica de las EDOs, encontrará un problema modelo común utilizado para estudiar los métodos numéricos es:

$y' = \lambda y$

Dónde $\lambda$ puede ser un número complejo. La solución general de esta EDO es obviamente

$y(t) = A e^{\lambda t}$

Lo que encontrarás es que los métodos numéricos básicamente siempre divergen de la solución real cuando $Re(\lambda) > 0$ .

Si miras hacia atrás en tu ODE, el $ty$ es análogo al término $\lambda y$ en el modelo ODE. Una cosa interesante a tener en cuenta es que desde $t$ es positivo durante la vida de esta simulación, esto significa, análogamente, que su ODE tiene un $Re(\lambda) > 0$ lo que significa que debes esperar que la solución diverja y por lo tanto explote básicamente sin importar el método numérico que utilices.