La mayoría de ecuaciones que aparecen en problemas científicos y técnicos no pueden resolverse de forma exacta, podemos plantearnos la obtención de una aproximación suficientemente buena de la solución.
Esta solución, denominada solución numérica, se determina mediante un método o algoritmo numérico que se ejecutará con la ayuda de un ordenador.
Este es el camino habitual que se sigue en las aplicaciones científicas y técnicas, donde surgen ecuaciones no lineales de difícil o imposible resolución analítica.
Nos centraremos en el estudio de un problema de valor inicial de primer orden:
y'(x) = f (x, y(x)), x en [a, b],
y(a) = y0.
Consideremos una discretización o mallado del intervalo [a, b], esto es, una particióndel mismo de la forma
a = x0 < x1 < x2 < < xN-1 < xN = b.
Por simplicidad, supondremos que el mallado es uniforme:
siendo h = la longitud de cada subintervalo:
[xk, xk+1] con k = 0, . . . , N-1
Es igual a un cierto valor h > 0 fijo.De esta forma, se verifica la siguiente relación:
xk+1 = xk + h, con k = 0, . . . , N-1.
El valor h se denomina paso de malla y los puntos xk son los nodos del mallado
La relación entre el paso de malla h y el número de subintervalos N es:
h = (b - a) / N
- Mientras menor sea el paso de malla (esto es, mientras más puntos tenga el mallado) mejor será la aproximación de la solución.
- Primer se fija el número de subintervalos N, que debe ser un número entero, y después se define el paso de malla usando la relación h = (b-a)/ N.
- Si damos directamente el paso de malla h y definimos N = (b-a) / h , podríamos obtener un valor no entero.
- Un método numérico para la resolución del problema de valor inicial es un algoritmo que permite obtener, para cada nodo xk del mallado, una aproximación yk del valor de la solución exacta en el nodo,
- De y(xk) puedo obtener una aproximación y(xk), con k = 0, . . . , N.
- Obtenemos el conjunto de puntos {y1, y2, . . . , yN}
Una interpretación física (interpretación cinética) del método de Euler la daría una partícula cuya posición en el instante x denotaremos por posición = y(x). Supongamos que la velocidad de la partícula puede expresarse como una
función f (x, y), dependiente del tiempo x y la posición y. En tal caso, y(x) es solución de la ecuación diferencial velocidad : y' = f (x, y).
La primera iteración del método de Euler se puede escribir como:
- y1 = y0 + (x1 - x0) f (x0, y0)
- Donde
- x1 - x0 = h es el paso de malla.
- Posición inicial de la parcicula: y0
- Posicón final de la partícula: y1
- Velocidad inicial de la partícula: f (x0, y0)
- Tiempo transcurrido x1 - x0
- Espacio recorrido por la partícula: (x1 - x0) f (x0, y0)
- entre los instantes x0 y x1, viajando a velocidad f (x0, y0).
Estamos pues aproximando la posición final de la partícula como la posición inicial más el espacio recorrido al viajar a una velocidad constante igual a la velocidad inicial.
En lugar de considerar la velocidad inicial f (x0, y0) como estimación de la velocidad en el intervalo [x0, x1], podríamos considerar una aproximación distinta. Por ejemplo, podríamos pensar en utilizar la media de las velocidades inicial y final, la velocidad en el punto medio, o incluso un promedio general de velocidades. Al considerar cada una de estas posibilidades se obtienen diferentes métodos numérico:
- Considerando la media entre las velocidades inicial f (x0, y0) y final f (x1, y1) de la partícula, obtenemos la expresión y1 = y0 +h/2·( f (x0, y0) + f (x1, p1)), que es la primera iteración del método de Heun.
- Considerado un promedio de cuatro velocidades diferentes: una en el nodo inicial xk, dos correcciones en el nodo intermedio xk + h/2 , y otra en el nodo final xk + h, obtenemos el método de Runge-Kutta.
El método de Euler se fundamenta en una sencilla idea geométrica: aproximar el valor de la solución en cada nodo por el valor que proporciona la recta tangente a la solución trazada desde el nodo anterior.
Ejemplo 5.1. Consideremos el problema de valor inicial
- y' = 1/2· (x^2 - y), con x en [0, 1],
- y(0) = 1
La solución exacta, obtenida mediante separación de variables, es
- y(x) = x^2 - 4x + 8 - 7e^(-x/2)
El método de Euler tiene en este caso la siguiente forma:
- yk+1 = yk + h/2·((xk)^2 - yk), con k = 0, 1, . . . , N - 1,
donde h = 1/N.
En la siguiente tabla se muestran, con una precisión de seis cifras decimales, los valores obtenidos con el método de Euler para N = 10 (paso de malla h = 0,1), así como los valores exactos de la solución y los errores cometidos. En este caso, el error global es igual a e = 0,025697.
En la figura siguiente se representan la solución exacta del problema y la solución numérica obtenida.
Ejemplo del método de Euler en Sage:
sage : var(’x, y’)
sage : fun(x, y) = (x^2-y)/2 # segundo miembro de la ecuacion
sage : a = 0. # extremo inicial intervalo
sage : b = 1. # extremo final intervalo
sage : y0 = 1. # valor inicial
sage : N = 10 # numero de particiones
sage : h = (b-a)/N # paso de malla
sage : eulers_method (fun , a, y0 , h, b)
x y h*f(x,y)
0.000000000000000 1.00000000000000 -0.0500000000000000
0.100000000000000 0.950000000000000 -0.0470000000000000
0.200000000000000 0.903000000000000 -0.0431500000000000
0.300000000000000 0.859850000000000 -0.0384925000000000
0.400000000000000 0.821357500000000 -0.0330678750000000
0.500000000000000 0.788289625000000 -0.0269144812500000
0.600000000000000 0.761375143750000 -0.0200687571875000
0.700000000000000 0.741306386562500 -0.0125653193281250
0.800000000000000 0.728741067234375 -0.00443705336171875
0.900000000000000 0.724304013872656 0.00428479930636718
1.00000000000000 0.728588813179023 0.0135705593410488
MÉTODO DE RUNGE-KUTTA DE CUARTO ORDEN RK4
Tiene la siguiente expresión:
- yk+1 = yk + h/6· (K1 + 2·K2 + 2·K3 + K4), con k = 0, 1, . . . , N - 1,
donde 8
- K1 = f (xk, yk)
- K2 = f (xk + h/2, yk + h/2·K1)
- K3 = f (xk + h/2 , yk + h/2·K2)
- K4 = f (xk + h, yk + h·K3)
En este caso se ha considerado un promedio de cuatro velocidades diferentes: una en el nodo inicial xk, dos correcciones en el nodo intermedio xk + h 2 , y otra en el nodo final xk + h.
El método de Runge-Kutta tiene orden cuatro, O(h^4).
Ejemplo 5.1. Consideremos el problema de valor inicial
- y' = 1/2· (x^2 - y), con x en [0, 1],
- y(0) = 1
La solución exacta, obtenida mediante separación de variables, es
- y(x) = x^2 - 4x + 8 - 7e^(-x/2)
Es claro que el método de Runge-Kutta produce resultados mucho más precisos que los métodos de
Euler o de Heun para el mismo paso de malla.
En Sage existe una implementación del método RK4, cuyo uso puede verse en el siguiente código:
sage : var(’x, y’)
sage : fun(x, y) = (x^2-y)/2 # segundo miembro de la ecuacion
sage : a = 0. # extremo inicial intervalo
sage : b = 1. # extremo final intervalo
sage : y0 = 1. # valor inicial
sage : N = 10 # particiones
sage : h = (b-a)/N # paso de malla
sage : desolve_rk4 (fun , y, ics =[a, y0], end_points =[a, b], step =h)
[[0.000000000000000 , 1.00000000000000] ,
[0.1 , 0.951394036458] ,
[0.2 , 0.906138090168] ,
[0.3 , 0.865044190327] ,
[0.4 , 0.828884763004] ,
[0.5 , 0.798394562606] ,
[0.6 , 0.774272509145] ,
[0.7 , 0.757183435905] ,
[0.8 , 0.747759751871] ,
[0.9 , 0.746603023076] ,
[1.0 , 0.754285476837]]
