Calcular experimentalmente el orden de error de Runge-Kutta 4

Question

El problema

Utilice el método Runge-Kutta de orden 4 para resolver la ecuación diferencial

$ \frac{\partial^2 y}{\partial t^2} = -g + \beta e^{-y/\alpha }*\left | \frac{\partial y}{\partial t} \right |^{2} $

Y corroborar que su error global es O (h^4)

El modelo matemático

Convierto el problema en un sistema de ecuaciones diferenciales de orden 1:

• $ \frac{\partial y}{\partial t} = v $
• $ \frac{\partial v}{\partial t} = -g + \beta e^{-y/\alpha }*\left | v \right |^{2} $

Por lo tanto, defino las variables de discretización u (para la posición) y v (para la velocidad) como:

• v = f(v, u, t)
• u = g(v, t)

Y utilizar los siguientes incrementos para el método Runge-Kutta de orden 4:

Para u

• k1v = h f(vn, un, tn)
• k2v = h f(vn + 0,5 k1v, un + 0,5 k1u, tn + 0,5 h)
• k3v = h f(vn + 0,5 k2v, un + 0,5 k2u, tn + 0,5 h)
• k4v = h f(vn + k3v, un + k3u, tn + h)

Para v

• k1u = h f(vn, tn)
• k2u = h f(vn + 0,5 k1v, tn + 0,5 h)
• k3u = h f(vn + 0,5 k2v, tn + 0,5 h)
• k4u = h f(vn + k3v, tn + h)

Y utilizarlos en la expresión RK4 para cada uno de ellos:

$ u_{n+1} = u_{n} + \frac{1}{6} (k_{1u} + 2 k_{2u} + 2 k_{3u} + k_{4u}) $

$ v_{n+1} = v_{n} + \frac{1}{6} (k_{1v} + 2 k_{2v} + 2 k_{3v} + k_{4v}) $

NOTA: Primero resuelvo para v . Para calcular el orden del error, resolveré 120 = h i veces con h = 0,1, h = 0,05 y utilizar el resultado dado para h = 0,001 como valor "real", ya que no conozco la función que resuelve la EDO. Entonces debo corroborar que el valor absoluto del "real" menos el resultado que obtuve de h = 0,1 debe ser 16 veces mayor que lo que obtengo al restar el valor que obtuve de h = 0,05 del valor "real".

El programa

Estoy usando C++ para resolver esto.

#include <iostream>
#include <math.h>
#include <cmath>
#include <sstream>
#include <fstream>
#include <vector>
#include <cstdlib>

long double rungeKutta(long double h)
{
    long double alpha = 6629;
    long double beta = 0.0047;

    long double pos = 39068;
    long double speed = 0;

    for (int i = 1; h*i < 120; i++)
    {
        long double k1v = h * (-9.8 + beta * exp(-pos/alpha) * pow(speed, 2));
        long double k1y = h * speed;
        long double k2v = h * (-9.8 + beta * exp(-(pos + 0.5*k1y)/alpha) * pow(speed + 0.5*k1v, 2));
        long double k2y = h * (speed + 0.5*k1v);
        long double k3v = h * (-9.8 + beta * exp(-(pos + 0.5*k2y)/alpha) * pow(speed + 0.5*k2v, 2));
        long double k3y = h * (speed + 0.5*k2v);
        long double k4v = h * (-9.8 + beta * exp(-(pos + k3y)/alpha) * pow(speed  + k3v, 2));
        long double k4y = h * (speed + k3v);

        speed = speed + (k1v + 2.0*(k2v + k3v) + k4v)/6;
        pos = pos + (k1y + 2.0*(k2y + k3y) + k4y)/6;
    }

    return pos;
}

int _tmain(int argc, _TCHAR* argv[])
{    
    long double errorOne = rungeKutta(0.01);
    long double errorTwo = rungeKutta(0.005);
    long double real = rungeKutta(0.0001);

    cout << fabs(real-errorOne) << endl << fabs(real - errorTwo) << endl;
    system("pause");
    return 0;
}

Los resultados

Me parece que el error sólo se reduce a la MITAD y no a la 1/16 parte del primer resultado.

¿Qué estoy haciendo mal? Me he quedado sin ideas.

Gracias.

Answer 1

Su problema puede ser que detenga el bucle en el momento en que $t=h*i\ge 120$ a partir de i=1 . Dado que los valores de h utilizado divide 120 (o cualquier número entero), esto significa que el número de pasos realizados puede ser uno fuera del número requerido debido a errores numéricos en el cálculo i*h . Esto puede dar un error de tamaño $h$ a la hora de finalización deseada $120$ que se refleja en un error de tamaño $h$ en los valores de la solución.

Por lo tanto, para estar absolutamente seguro de que se integra la hora correcta, defina N=floor(120/h) para que $Nh\le 120<N+1$ , bucle para i=0 a N con dt=h para i=0 a N-1 y en el último paso i==N , set dt=120-N*h .

---

Y efectivamente, si se sigue el tiempo real durante la integración introduciendo t=0 antes del bucle y actualizando t+=h dentro del bucle, se encontrará que el bucle termina en t==120-h . Una alternativa al uso del número real de pasos en el bucle es cambiar la condición del bucle por (i-0.5)*h<120 para que los errores de redondeo en h ser compensado.

Y entonces descubrirás que con h=0.01 y h=0.005 ya has superado el punto de equilibrio en el que los errores de punto flotante acumulados tienen un peso mayor que los errores de discretización. Comparando h=0.4 y h=0.2 resulta en el factor esperado de 16.

Answer 2

Usted ha definido:

• u' = g(v, t )

pero lo usas como

• k1u = h f(vn, un)

en lugar de

• k1u = h g(vn, tn )

Answer 3

Esta es una vieja pregunta mía, pero la contestaré ya que veo que ha tenido alguna actividad reciente. El problema surgió debido a una mala interpretación del problema y de los datos dados. Después de volver a ejecutar todas las simulaciones que alimentaron los datos a mi función, funcionó correctamente.

No recuerdo si fue exactamente la orden 4 la que obtuve, así que puede haber algún otro problema con mi función, pero como era lo suficientemente cercana no seguí buscando.