¿Solucionador de EDO numérico robusto?

Question

Hice un poco explícita de los métodos de Runge-Kutta de 4º orden solver hace un par de días, pero cuando las pruebas en contra de varios de 1er y 2º orden Odas elegido al azar (por ejemplo $d^{2}y/dt^{2} = -y \sin(y)$, $d^{2}y/dt^{2} = -yt$ o $d^{2}y/dt^{2} = -y + t^{2}$) parecía que la mayoría eran rígidos Odas (a menos que el algoritmo estoy usando es incorrecta), mediante la comparación de mi salida para que de Mathematica NDSolve, y por lo tanto prestados mi RK4 solver ser bastante inútil. Como tal, me he decidido a probar y encontrar un solucionador numérico que puedo crear, que es robusto y puede resolver rígido y no rígido Odas. Qué tal un algoritmo existe, o es un caso de los más robustos un solver el más abstruso su algoritmo se convierte en. Mejor aún, ¿existe tal cosa como un universal solver que es capaz de resolver cualquier ODE que lanzar en él?

EDIT: he Aquí un ejemplo de mi RK4 solver de salida para $d^{2}y/dt^{2} = -y \sin(y)$ usando un tamaño de paso de $h=0.005$:

Y he aquí lo que me sale en NDSolve:

Answer 1

He sido capaz de reproducir tanto de sus fotos con NDSolve. La segunda, suavemente-buscando uno es la solución de $y''(t)=-y(t)\sin(y(t))$ con condiciones iniciales $y(0)=0$$y'(0)=50$.

Llego a la primera, si me parcela de la derivada de la solución:

Así, parece que usted está tomando mal la salida de su correcta solver. Como es un método de Runge-Kutta, que es más probable la división de la ecuación en el sistema de dos ecuaciones, una para $y'(t)$ y otro para $y(t)$. Si usted está tomando la primera como la solución a la ecuación original, mientras que usted tiene que tomar el último.

Answer 2

Sí, este tipo de algoritmos existen.

Implícito de varios pasos se han desarrollado métodos deliberadamente para resolver rígido ODE problemas. Se puede requerir un método explícito para inicializar, pero este no es un problema porque puede ampliar el tamaño del paso de su explícito de solver en la rigidez de la región, más discretamente.

En la no-rígido escalas de un problema, el implícito de varios pasos son los métodos computacionalmente un poco más caro, pero en última instancia será el rendimiento como una precisa una solución como algo parecido a un método RK4.

Alternativamente, usted puede utilizar una adaptación de solver, tales como el método de Runge-Kutta-Fehlberg o un Dormand-Prince 4(5) par.

Answer 3

He probado a hacer un segundo pedido ODE solver en ocaml:

let secondode n=
    let rec iode i a b c=
        if i=n then a
        else iode (i+1) (a+.0.001*.b) (b+.0.001*.c) (-.a*.(sin a))
    in iode 0 0. 1.75 0.
;;

secondode 200000;;

y el resultado es $\textrm{float} = -75208.630048495936$, la cual difiere significativamente del valor real $\approx -1.75$. Así que asumo que mi ingenua algoritmo iterativo no es estable. Aquí me eligen deliberadamente $0.001$ ya que no es demasiado pequeño y no demasiado grande. Si puedo elegir exactitud a ser $10^{-5}$, entonces el resultado es $\textrm{float} = -2.050306346414887$, que todavía se distingue mucho. Y precisión en $10^{-6}$, mi resultado es de alrededor de $-1.80$.

Por eso creo que el método de Runge-Kutta es bastante estable (en comparación con mi ingenuo método). Pero yo realmente no sé cómo Mathematica que....