Clase No. 23:
Diferenciación numérica:Método de Euler explícito
Métodos basados en la serie de TaylorMAT–251 Dr. Alonso Ramírez Manzanares
CIMAT A.C.e-mail: [email protected]: http://www.cimat.mx/salram/met_num/
Dr. Joaquín Peña AcevedoCIMAT A.C.e-mail: [email protected]
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 1 / 27
Introducción (I)
• Una ecuación diferencial es una ecuación que relaciona a una funciónde una o varias variables con sus derivadas.
• Si la solución es una función de una sola variable, entonces decimosque la ecuación es ordinaria.
• De manera general, una ecuación ordinaria es de la forma
dny
dxn= f (x, y(x), y′(x), ..., y(n−1)(x) ).
• La derivada de mayor orden determina el orden de la ecuacióndiferencial.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 2 / 27
Introducción (II)• También podemos tener un sistema de ecuaciones diferenciales
ordinarias. En este caso tenemos n funciones y1(x),y2(x), ...,yn(x) talesque
dy1
dx= f1(x, y1(x), y2(x), ... ,yn(x) )
dy2
dx= f2(x, y1(x), y2(x), ... ,yn(x) )
......
...dyn
dx= fn(x, y1(x), y2(x), ... ,yn(x) )
Escrito de manera más compacta,
dy
dx= f(x,y(x)), con y(x) = (y1(x), y2(x), ... ,yn(x)),
Así, y : R→ Rn, y f : Rn+1 → Rn.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 3 / 27
Introducción (III)
• Hay que especificar el intervalo en que queremos encontrar la solución.
• También necesitamos una condición auxiliar para encontrar la solucióndel problema.
• Un problema de valor inicial (PVI) para una ecuación diferencial deprimer orden:
dy
dx= f(x,y(x)), x ∈ [x0,x1]
y(x0) = y0.
(1)
• Una ecuación de orden n la podemos transformar en un sistema deecuaciones de primer orden mediante
y1(x) = y(x), yi+1(x) = y(i)(x), i = 1, ...,n− 1.
y1(a) = y(a), · · · ,yn−1(a) = y(n−1)(a)
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 4 / 27
Existencia y unicidad
• Tenemos el teorema de Peano, que nos dice que si f(x,y) es continua enun dominio Ω, que contiene al punto inicial (x0,y0), entonces la soluciónexiste en un intervalo alrededor de x0.
• La unicidad puede estar garantizada si la función f es Lipschitz continuaen Ω respecto a y. Esto es, existe L tal que
‖f(x,y1)− f(x,y2)‖ < L‖y1 − y2‖
para todo (x,y1), (x,y2) ∈ Ω.
• El cálculo de la solución del problema de valor inicial puede requerir eluso de algún método numérico.
Teorema de existencia y unicidad
Si existe un dominio Ω que contiene a la condición inicial (x0,y0) en el quela función f(x,y) es continua y es Lipschitz en y, entonces existe una únicasolución al PVI (1) en el intervalo [x0 − δ,x0 + δ] con δ > 0.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 5 / 27
Métodos numéricos de solución
• En general, se define el intervalo [a,b] en el que se quiere calcular lasolución.
• Se discretiza el intervalo
a = x0 < x1 < x2 < · · · < xn = b
para tratar de calcular el valor de la solución en puntos xi, yi ≈ y(xi).
• Frecuentemente la distribución de los nodos es uniforme, por lo que elvalor h = xi+1 − xi es constante.
• Hay métodos los que se trata de calcular todos los valores de lasvariables yi al mismo tiempo, lo cual conduce a resolver sistemas deecuaciones.
• Otros métodos se basan en estimar el valor yi a partir de valores previosyi−1,yi−2, ..., ya determinados. Como esto permite calcular el valor enxi = xi−1 + h, a h se le llama el tamaño de paso.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 6 / 27
Relación con integración
Resolver una ecuación diferencial y el cálculo de integrales estánrelacionados. Consideremos el problema de valor inicial:
dy
dx= f (x,y) x ∈ [a,b]
y(a) = y0
Integrando de t a t + h ambos miembros obtenemos
∫ t+h
t
dy
dxdx =
∫ t+h
tf (x,y) dx =⇒ y(t + h) = y(t) +
∫ t+h
tf (x,y(x)) dx
Si reemplazamos la integral por alguna de las reglas de cuadratura vistasanteriormente podemos obtener una fórmula para resolver la ecuacióndiferencial.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 7 / 27
Relación con integración
Resolver una ecuación diferencial y el cálculo de integrales estánrelacionados. Consideremos el problema de valor inicial:
dy
dx= f (x,y) x ∈ [a,b]
y(a) = y0
Integrando de t a t + h ambos miembros obtenemos
∫ t+h
t
dy
dxdx =
∫ t+h
tf (x,y) dx =⇒ y(t + h) = y(t) +
∫ t+h
tf (x,y(x)) dx
Si reemplazamos la integral por alguna de las reglas de cuadratura vistasanteriormente podemos obtener una fórmula para resolver la ecuacióndiferencial.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 7 / 27
Relación con series de Taylor
Suponemos que y tiene una expansión en series de Taylor:
y(x + h) = y(x) + hy′(x) +h2
2y′′(x) + ...+
hm
m!y(m)(x) + ...
• Lo usual es truncar la serie después de cierta cantidad de términos.
• Si truncamos después de m + 1 términos, decimos que obtenemos unmétodo basado en series de Taylor de orden m.
• Lo que resta es aproximar las derivadas que permanecen en laexpresión o calcularla una a partir de la anterior.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 8 / 27
Método de Euler (I)
Queremos resolver el problema de valor inicial
dy
dx= f (x,y) x ∈ [a,b]
y(a) = y0
Al usar un método basado en series de Taylor de orden 1, obtenemos elmétodo de Euler:
y(x + h) ≈ y(x) + hy′(x) = y(x) + hf (x,y(x))
Dividimos el intervalo en n partes iguales:
a = x0 < x1 < · · · < xn = b, xi = xi−1 + h, h =b− a
n.
Si denotamos por yi al valor que se obtiene por el el método de Euler en xi,entonces
yi = yi−1 + hf (xi−1,yi−1), i = 1,2, ...,n.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 9 / 27
Método de Euler (II)
Otra forma de obtener el método de Euler es mediante integración,aproximando la integral por un rectángulo cuya altura es igual al integrandoevaluado en el límite inferior de la integral:
y(x + h) = y(x) +
∫ x+h
xf (t,y(t)) dt ≈ y(x) + hf (x,y(x))
Ejemplo 1. Considere el problema de valor inicial
dy
dx= y x ∈ [0,5]
y(0) = 1
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 10 / 27
Método de Euler (III)
0 1 2 3 4 5
050
100
150
c(0, 5)
c(0,
max
(vs)
)
0 1 2 3 4 5
050
100
150
c(0, 5)c(
0, m
ax(v
s))
n = 50 n = 150
Comparación entre la solución analítica y la que se obtiene con el métodode Euler.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 11 / 27
Método de Euler (IV)
Ejemplo 2. Considere el problema de valor inicial
dy
dx=p
y x ∈ [0,2]
y(0) = 0
Tenemos dos soluciones:
y(x) =x2
4,
y
y(x) = 0.
Si aplicamos el método de Euler con n = 20, se tiene el siguiente resultado
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 12 / 27
Método de Euler (V)
0.0 0.5 1.0 1.5 2.0
0.0
0.2
0.4
0.6
0.8
1.0
c(0, b)
c(0,
max
(vs)
)
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 13 / 27
Método de Euler (VI)
Ejemplo 3. Considere el problema de valor inicial
dy
dx= ycosx x ∈ [0,8π]
y(0) = 1
0 5 10 15 20 25
0.0
0.5
1.0
1.5
2.0
2.5
3.0
X
Y
h = 0.09856
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 14 / 27
Método de Euler (VII)
0 5 10 15 20 25
0.0
0.5
1.0
1.5
2.0
2.5
3.0
X
Y
h = 0.00614
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 15 / 27
Interpretación geométrica (I)
En el método de Euler explícito se da un paso de tamaño h en dirección dela tangente:
y′ = − sinx, y(0) = 1
0 2 4 6 8
−1.
0−
0.5
0.0
0.5
1.0
X
Y
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 16 / 27
Estabilidad de la solución (I)
De esta forma, la solución numérica del PVI requiere una discretización delintervalo [a,b]:
a = x0 < x1 < · · · < xn = b
Tenemos que h(i) = xi+1 − xi, para i = 0, ...,n− 1.
• Si h(i) es constante, h(i) = (b− a)/n.
• Si h(i) es variable, decimos que el método de solución es adaptativo.
Vamos a denotar por yi a la aproximación numérica de la solución del PVI enxi, mientras que y(xi) denota a la solución exacta.
La solución numérica está formada por el conjunto de aproximacionesyin
i=1.
Estabilidad
Decimos que la solución numérica yini=1 para el PVI con condición inicial
y(a) = y0 es estable si para cualquier ε > 0 existe δ > 0 tal que yini=1 es
la solución numérica con condición inicial y(a) = y0 y ‖y0 − y0‖ < δ, se debetener que ‖yi − yi‖ < ε para todo i = 1, ...,n.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 17 / 27
Estabilidad de la solución (II)
Se dice que la solución numérica es asintóticamente estable cuando
‖y0 − y0‖ < δ =⇒ limi→∞‖yi − yi‖ = 0.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 18 / 27
Estabilidad del método de Euler
Consideremos el PVI:
dy
dx= λy
y(0) = 0
Entonces
yk+1 = yk + hλyk
De aquí que
yk+1 = (1 + hλ)k+1y0
De la definición de estabilidad, para que el método de Euler sea estable sedebe cumplir
|1 + hλ| ≤ 1.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 19 / 27
Estabilidad del método de Euler
Consideremos el PVI:
dy
dx= λy
y(0) = 0
Entonces
yk+1 = yk + hλyk
De aquí que
yk+1 = (1 + hλ)k+1y0
De la definición de estabilidad, para que el método de Euler sea estable sedebe cumplir
|1 + hλ| ≤ 1.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 19 / 27
Error del método de Euler
También se puede definir el error local de truncamiento como
L(x,y; h) = y(x + h)− [y + hf (x,y)].
Puesto que y(x + h) = y(x) + hy′(h) + h2
2 y′′(ξ), entonces
L(x,y; h) =h2
2y′′(ξ) = O(h2)
Este error determina como afecta la aproximación a los pasos subsecuentes,es decir, la propagación del error.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 20 / 27
Fuentes de error
Las más importantes son
• el error de redondeo,rn = |fl(yi)− yi|.
• el error de truncamiento,ei = |yi − y(xi)|
El error total está dado por
|fl(yi)− y(xi)| = |fl(yi)− yi + yi − y(xi)| ≤ |fl(yi)− yi|+ |yi − y(xi)| = ri + ei.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 21 / 27
Métodos basados en series de Taylor (I)
Para obtener un método basado en series de Taylor de orden superiorpodemos simplemente derivar la función f .
Ejemplo.
dy
dx=
x
y + 1x ∈ [0,4]
y(0) = 0
Puesto que y′ = xy+1 , entonces
y′′ =1
y + 1−
xy′
(y + 1)2
y′′′ = −y′
(y + 1)2−
(y + 1)(y′ + xy′′)− 2x(y′)2
(y + 1)4
Así, podemos obtener un método de tercer orden:
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 22 / 27
Métodos basados en series de Taylor (II)
1 Tenemos el valor yi en xi.
2 Calculamos y′i
=xi
yi + 1.
3 Calculamos y′′i
=1
yi + 1−
xiy′i
(yi + 1)2
4 Calculamos y′′′i
= −y′i
(yi + 1)2−
(yi + 1)(y′i + xiy′′i )− 2xi(y
′i )2
(yi + 1)3
5 Calculamos yi+1 = yi + hy′i
+h2
2y′′
i+
h3
6y′′′
i.
Comparamos este método con Euler y la solución analítica, la cual esy =p
x2 + 1− 1.
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 23 / 27
Métodos basados en series de Taylor (III)
Euler Taylor de orden 3
0 1 2 3 4
0.0
0.5
1.0
1.5
2.0
2.5
3.0
c(0, b)
c(0,
max
(vs)
)
0 1 2 3 40.
00.
51.
01.
52.
02.
53.
0
c(0, b)
c(0,
max
(vs)
)
n = 10 n = 10
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 24 / 27
Método trapezoidal (I)
Regresando a la relación entre resolver EDOs e integrales:
y(t + h) = y(t) +
∫ t+h
tf (x,y(x)) dx
Si aproximamos la integral por la regla del trapecio, obtenemos un métodoimplícito:
yi+1 = yi +1
2h[f (xi,yi) + f (xi+1,yi+1)]
En este caso es un método implícito de un paso. En general, esto nosconduce a un sistema de ecuaciones.
Ejemplo. Para el problema de valor inicial planteado anteriormente
dy
dx= ycosx x ∈ [0,8π]
y(0) = 1
Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 09.11.2015 25 / 27