En esta publicación les presento un material diseñado por mi persona para el curso de calculo numerico de la Universidad de las Fuerzas Armadas (UNEFANE)
RESUMEN
En esta publicación se hizo una introducción a las ecuaciones diferenciales y su terminología, se definieron sus características: tipo, orden. Grado y linealidad. Luego se pasó a estudiar los métodos numéricos de Euler y Runge Kutta para resolver ecuaciones lineales de primer orden con valores iniciales Los conceptos básicos de ecuaciones diferenciales los pueden ver de los siguientes videos:
INTRODUCCION A LAS ECUACIONES DIFERENCIALES, CONCEPTOS BASICOS
[youtube 
]
EJEMPLO DE SOLUCION DE UNA ECUACION DIFERENCIAL
%
En este ejercicio se demuestra que la función y=3x+x2 es solución de la ed xy’ – y = x2
En este curso solamente se estudian los métodos numéricos para resolver ecuaciones lineales de primer orden con valores iniciales. Par estos problemas existen básicamente dos métodos: el de Euler y Runge Kutta. La teoría de estos métodos están bien explicadas en la bibliográfica anexa, de las cuales le recomiendo leer al menos una de ellas para tener una mejor comprensión de estos métodos. Otra buena exposición de estos métodos la pueden encontrar en:
https://www2.caminos.upm.es/Departamentos/matematicas/Fdistancia/PIE/Analisis%20matematico/Temas/C13_Metodo_Numerico.pdf
METODO DE EULER
Definiendo un problema de valor inicial como:
dy/dx=f(x,y), y(0)=a
El método de Euler resuelve este problema usando la formula iterativa:
yi+1 = yi + hf(xi,yi)
Ejemplo:
Sea y’= 2xy. Calular la función y para x=1.5 dados los valores iniciales: y(1)=1 y h=0,1
Solución. En este problema a=1, b=1.5, h=0.1. Notase que a es el valor inicial de x tomado de la condición inicial y(1)=1 y b es el valor en donde queremos calcular y, h se denomina el paso y es la longitud de deltax
Cálculo de número de iteraciones. En este problema nos dan a, b y h=deltax. De la formula h=(b-a)/n, despejamos n, y obtenemos: n=(b-a)/h=n=(1.5-1)/.1=5
2. Cálculo de la partición a=1, h=0.1 ->x0=1; x1=1.1; x2=1.2; x3=1.3; x4=1.4 y x5=1.5
3. Proceso iterativo: Calculo de la tabla de valores de yi para cada xi
Comenzamos con los valores iniciales x0=1 y y0=1. Usando la formula
yi+1 = yi + hf(xi,yi) para i = 0, 1, 2,..n-1, obtenemos:
i=0 -> y1=y0+hf(x0,y0) = 1+0.1(2x0y0)=1+0.1(21*1)=1.2
i=1 -> y2=y1+hf(x1,y1) = 1.2+0.1(2x1y1)=1.2+0.1(21.1*1.2)=1.464
i=2 -> y3=y3+hf(x2,y2) = 1.464+0.1(2x2y2)=1.464+0.1(21.2*1.464)=1.81536
i=3 -> y4=y3+hf(x3,y3) = 1.81536+0.1(2x3y3)=1.81536+0.1(21.3*1.81536)=2.287354
i=4 -> y5=y4+hf(x4,y4) = 2.287354+0.1(2x4y4)=2.287354+0.1(21.4*2.287354)=2.927813
RESULTADO: Y(1.5)=2.927813
Tabla de resultados parciales
i Xi Yi+1
0 1 1.2
1 1.1 1.464
2 1.2 1.81536
3 1.3 2.287354
4 1.4 2.927813
Programa en R
#metodo de Euler
f= function (x,y) 2xy
a=1; b=1.5; h=0.1
#inicialización de los vectores
x=c(0); y=c(0); yp=c(0)
La variable yp = f(xi,yi) (yp es yprima)
#calculo de n
n=(b-a)/h
x[1]=a; y[1]=1
#Estos valores de x y y corresponden
#a los valores iniciales de y y(0)=1,
#porque los vectores en R comienzan desde 1 y no desde cero
cat("\n", "Particion del intervalo [a,b]")
for (i in 2:(n+1)) {
x[i]=x[i-1]+h
cat( "\n",i-1,x[i])}
#proceso iterativo
for (i in 2:(n+1)){yp[i-1]=f(x[i-1],y[i-1])
y[i]=y[i-1]+h*yp[i-1]}
y[i+1]=y[i]+h*f(x[i],y[i])
Tabla de resultados
cat("\n", "Tabla de Resultados")
cat("\n", "i x y")
for (i in 2:(n+2)) {
cat( "\n",i-1,x[i-1],y[i-1])}
corrida del programa
source('~/Euler.R')
Particion del intervalo [a,b]
1 1.1
2 1.2
3 1.3
4 1.4
5 1.5
Tabla de Resultados
i x y
1 1 1
2 1.1 1.2
3 1.2 1.464
4 1.3 1.81536
5 1.4 2.287354
6 1.5 2.927813
Método de Runge-Kutta.
Es un método numérico de resolución de ecuaciones diferenciales que surge como una mejora del método de Euler, el cual se puede considerar como un método de Runge Kutta de primer orden, éste método logra la exactitud de una serie de Taylor pero sin requerir el cálculo de derivadas superiores. El método Runge-Kutta no es sólo un método sino una importante familia de métodos iterativos tanto implícitos como explícitos para aproximar las soluciones de ecuaciones diferenciales ordinarias.
Con el método RK4, obtener una aproximación del valor de y(1,5) para el siguiente problema de valor inicial, tomando un paso h = 0,1.
f= function (x,y) 2xy
a=1; b=1.5; h=0.1
Solución:
El primer paso para resolver este problema es determinar la particion de puntos en donde se va a obtener la solución.
Como en este conocemos a,b y h=deltax, calculamos n, se tiene que
n = (b - a)/h= (1,5 - 1)/0,1 = 5.
Por lo tanto, los puntos en donde se va a determinar la solución, dados por la fórmula
xi = 1 + 0,1 i, para i =1,2,3,4,5, son: x1 = 1.1 ; x2 = 1.2; 3 = 1.3; x4 = 1.4 y x5 = 1.5
Resulta entonces, que para cada i se aplican las formulas. Se comienza desde i=2 ya que en R los vectores comienzan por 1, y no por cero, entonces y(1) en el programa es y0 en el problema
k1=f(x[i-1],y[i-1])*h
k2=hf(x[i-1]+(1/2)h,y[i-1]+(1/2)*k1);cat("\n",i,"k2=",k2)
k3=hf(x[i-1]+(1/2)h,y[i-1]+(1/2)*k2);cat("\n",i,"k3=",k3)
k4=h*f(x[i-1]+h,y[i-1]+k3);cat("\n",i,"k4=",k4)
y[i]=y[i-1]+(1/6)(k1+2k2+2*k3+k4)}
Y aplicando sucesivamente las fórmula de RK4, para i desde 1 hasta 4, se obtienen los datos que se muestran en la siguiente tabla,
i x RK
1 1 1
2 1.1 1.233
3 1.2 1.552
4 1.3 1.993
5 1.4 2.611
6 1.5 3.490
donde además se muestra el valor de la solución exacta para cada punto de la partición
Solución exacta.
La ecuación diferencial y´ =2xy se resuelve fácilmente por el método de variables separables:
y´ = 2xy -> y´/y =2x , ahora integrando cada miembro de la igualdad, tenemos:
ln(y)=x^2 +c, ahora aplicando la función exponencial, tenemos la solución general: y=exp(x^2+c)
Para calcular c, usamos las condiciones iniciales x=1, y=1 y obtenemos 1=exp(1-c), aplicando logaritmo, se obtiene: 0=1+c -> c=-1. Por lo tanto la ecuación particular es: y=exp(x^2-1)
Ahora evaluemos y para x=1.5
y=exp(1.5^2-1)=3.490343
Al analizar la tabla anterior y comparar los resultados obtenidos con el método RK4 con los valores reales, se ve por qué es tan difundido este método. En la próxima tabla se comparan los métodos de Euler y Runge Kutta de orden 4 para un mismo problema.
Programa en R
#metodo de RK
f= function (x,y) 2xy
a=1; b=1.5; h=0.1
#inicialización de los vectores
x=c(0); y=c(0); yp=c(0)
La variable yp = f(xi,yi) (yp es yprima)
#calculo de n
n=(b-a)/h
x[1]=a; y[1]=1
#Estos valores de x y y corresponden
#a los valores iniciales de y y(0)=1,
#porque los vectores en R comienzan desde 1 y no desde cero
cat("\n", "Particion del intervalo [a,b]")
for (i in 2:(n+1)) {
x[i]=x[i-1]+h
cat( "\n",i-1,x[i])}
#proceso iterativo
for (i in 2:(n+2)){k1=f(x[i-1],y[i-1])*h
cat("\n",i,"k1=",k1)
k2=hf(x[i-1]+(1/2)h,y[i-1]+(1/2)*k1);cat("\n",i,"k2=",k2)
k3=hf(x[i-1]+(1/2)h,y[i-1]+(1/2)*k2);cat("\n",i,"k3=",k3)
k4=h*f(x[i-1]+h,y[i-1]+k3);cat("\n",i,"k4=",k4)
y[i]=y[i-1]+(1/6)(k1+2k2+2*k3+k4)}
Tabla de resultados
cat("\n", "Tabla de Resultados")
cat("\n", "i x y")
for (i in 1:(n+1)) {
cat( "\n",i,x[i],y[i])}
salida
i x y RK
1 1 1
2 1.1 1.233674
3 1.2 1.552695
4 1.3 1.993687
5 1.4 2.611633
6 1.5 3.49021
EJERCICIOS
LOS EJERCICIOS A RESOLVER POR LOS METODOS DE EULER Y RK4
Ecuaciones diferenciales -ejercicio9.1
Para este y temas similares consulta mi blog calnum.wordpress.com
Bibliografia
Steven C. Chapra; Raymond P. Canale. Métodos numéricos para ingenieros. Quinta edición. México. Mac Gaw Hill. 2007.
Dennis Zill ; Michael Cullen. Ecuaciones diferenciales. Séptima edición. México. Cengage Learning Editores, S. A. de C. V.,2009