Transcript Soluciones Numéricas

Soluciones Numéricas de Ecuaciones
Diferenciales ordinarias
CÁPITULO 6
Contenidos
•
•
•
•
•
6.1 Método de Euler y Análisis de Error
6.2 Métodos de Runge-Kutta
6.3 Métodos de Varios Pasos
6.4 Ecuaciones de Orden Superior y Sistemas
6.5 problemas de Valores en al Frontera de
Segundo Orden
6.1 Método de Euler y Análisis de Error
• Introducción
Recuerde la estructura del Método de Euler
yn+1 = yn + hf(xn, yn)
(1)
• Errores en Métodos Numéricos
Una fuente de error que siempre está
presente en los cálculos es el error de rondeo.
Errores de Truncamiento para el
Método de Euler
• Este logaritmo sólo da una aproximación en línea
recta a la solución. Este error se llama error de
truncamiento local, o error de discretización.
Para obtener una fórmula para el error de
truncamiento local del método de Euler, usamos
la fórmula de Taylor con resido.
y ( x)
xa
( x  a)k
( x  a ) k 1
(k )
( k 1)
 y (a )  y(a )
   y (a)
y
(c )
1!
k!
(k  1)!
Donde c es un punto entre a y x.
• Tomando k = 1, a = xn, x = xn+1 = xn + h, tenemos
h
h2
y ( xn1 )  y ( xn )  y( xn )  y(c)
1!
2!
ó
h2
y ( xn1 )  yn  hf ( xn , yn )  y(c)

2!
yn 1
De ahí que el error de truncamiento en yn+1 es
h2
y(c)
2!
donde xn < c < xn+1
El valor de c por lo común no se conoce, pero una
cota superior es
h2
| y( x) |, xn  x  xn1
M
donde M  x max
2!
n  x xn 1
• Observación:
Se dice que e(h) es de orden hn, representado
con O(hn), si existe una constante C tal que
|e(h)|  Chn para h suficientemente pequeña.
Ejemplo 1
Determine una cota para los errores de truncamiento
local del método de Euler aplicado a
y '  2 xy, y (1)  1
Solución
x 2 1
2 x 2 1
De la solución y  e , tenemos y  (2  4 x )e ,
so
2
2
h2
h
y(c)  (2  4c 2 )e( c 1)
2
2
En particular, para h = 0.1, se puede obtener una cota
superior remplazando c por 1.1 es
2
2
(
0
.
1
)
[2  (4)(1.1)2 ]e((1.1) 1)
 0.0422
2
Ejemplo 1 (2)
Al hacer 5 pasos, remplazando c por 1.5, se
obtiene
[2  (4)(1.5)2 ]e((1.5)
2
2
(
0
.
1
)
1)
 0.1920
2
(2)
Método de Euler Mejorado
f ( xn , yn )  f ( xn1 , yn*1 )
yn1  yn  h
2
donde
yn*1  yn  hf ( xn , yn )
(3)
(4)
se conoce comunmente como el Método de Euler
Mejorado. Fig 6.1
En general, el método de Euler mejorado es un
ejemplo de método de predictor y corrector.
Fig 6.1
Ejemplo 2
Use el método de Euler mejorado para obtener el
valor aproximado de y(1.5) para la solución de
y '  2 xy, y (1)  1. Compare los resultados para h = 0.1
y h = 0.05.
Solución
Con x0 = 1, y0 = 1, f(xn, yn) = 2xnyn , h = 0.1
y1* = y0 + (0.1)(2xy) = 1.2
Usando (3) con x1 = 1 + h = 1.1
2 x0 y0  2 x1 y1*
2(1)(1)  2(1.1)(1.2)
y1  y0  (0.1)
 1  (0.1)
 1.232
2
2
Los resultados se dan en la Tabla 6.3 y 6.4.
Tabla 6.3
xn
yn
Valor
real
Error
Abs.
% error
relativo
1.00
1.10
1.20
1.0000
1.2320
1.5479
1.0000
1.2337
1.5527
0.0000
0.0017
0.0048
0.00
0.14
0.31
1.30
1.40
1.50
1.9832
2.5908
3.4509
1.9937
2.6117
3.4904
0.0106
0.0209
0.0394
0.53
0.80
1.13
Tabla 6.4
xn
yn
Valor
real
Error
Abs
% error
relativo
1.00
1.0000
1.0000
0.0000
0.00
1.05
1.1077
1.1079
0.0002
0.02
1.10
1.2332
1.2337
0.0004
0.04
1.15
1.3798
1.3806
0.0008
0.06
1.20
1.5514
1.5527
0.0013
0.08
1.25
1.7531
1.7551
0.0020
0.11
1.30
1.9909
1.9937
0.0029
0.14
1.35
2.2721
2.2762
0.0041
0.18
1.40
2.6060
2.6117
0.0057
0.22
1.45
3.0038
3.0117
0.0079
0.26
1.50
3.4795
3.4904
0.0108
0.31
• Errores de Truncamiento para el Método
Mejorado de Euler
Observe que el error de truncamiento local es
O(h3).
6.2 Runge-Kutta Methods
• Métodos de Runge-Kutta
Todos los métodos de Runge-Kutta son
generalizaciones de la fórmula básica de Euler, en la
que la función pendiente f se remplaza por un
promedio ponderado de pendientes en el intervalo
xn  x  xn+1
yn1  yn  h( w1k1  w2k2    wm km )
(1)
donde las ponderaciones wi, i = 1, 2, …, m son
constantes que satisfacen w1 + w2 + … + wm = 0, y ki
es la función evaluada en un punto seleccionado (x, y)
para el cual xn  x  xn+1.
El número m se llama el orden. Si tomamos
m = 1, w1 = 1, k1 = f(x, yn), llegamos al método
de Euler. Por consiguiente, se dice que el
método de Euler es un método de RungeKutta de primer orden.
Método de Runge-Kutta de Segundo Orden
• Tratamos de hallar unas constantes de modo que la
fórmula
yn1  yn  ak1  bk2
(2)
donde k1= f(xn, yn),
k2= f(xn+h, yn+hk1)
concuerde con un polinomio de Taylor de grado 2.
Las constantes deben satisfacer
(3)
1
1
w1  w2  1, w2  , y w2  
2
2
luego
1
1
(4)
w1  1  w2 ,  
, y 
donde w2  0.
2w2
2w2
Ejemplo: escogemos w2 = ½ , de donde w1 = ½ ,
 = 1,  = 1, y (2) se transforma en
yn+1= yn+(k1+ k2)h/2
donde k1= f(xn, yn), k2= f(xn+h, yn+hk1).
Puesto que xn + h = xn+1, yn + hk1 = yn + hf(xn, yn),
es idéntica al método de Euler mejorado.
Método de Runge-Kutta de Cuarto Orden
• Tratamos de hallar parámetros de modo que la
fórmula
yn1  yn  ak1  bk2  ck3  dk4
(5)
donde k  hf ( xn , yn )
k2  hf ( xn  1h, yn  1k1 )
k3  hf ( xn   2 h, yn   2 k1  3k2 )
k4  hf ( xn h, yn  k3 )
concuerde con un polinomio de Taylor de orden 4.
El conjunto de valores usado con más
frecuencia para los parámetros produce el
siguiente resultado
1
yn1  yn  (k1  2k2  2k3  k4 )
6
k1  hf ( xn , yn )
k2  hf ( xn  1/2 h, yn  1/2 k1 )
k3  hf ( xn  1/2 h, yn  1/2 k2 )
(6)
Ejemplo 1
Use el método RK4 con h = 0.1 para obtener y(1.5) para
la solución de y’ = 2xy, y(1) = 1.
Solución
Primero se calcula el caso n = 0.
k1  (0.1) f ( x0 , y0 )  (0.1)( 2 x0 y0 )  0.2
k2  (0.1) f ( x0  1/2(0.1), y0  1/2(0.2))
 (0.1)2( x0  1/2(0.1))( y0  1/2(0.2))  0.231
k3  (0.1) f ( x0  1/2(0.1), y0  1/2(0.231))
 (0.1)2( x0  1/2(0.1))( y0  1/2(0.231))  0.234255
k4  (0.1) f ( x0  0.1, y0  0.234255 )
 (0.1)2( x0  0.1)( y0  0.234255 )  0.2715361
Ejemplo 1 (2)
Por lo tanto,
1
y1  y0  (k1  2k2  2k3  k4 )
6
1
 1  (0.2  2(0.231)  2(0.234255 )  0.2715361)
6
 1.23367435
Véase la Tabla 6.5.
Tabla 6.5
h=0.1
xn
yn
Valor
real
Error
Abs.
% error
relativo
1.00
1.10
1.20
1.30
1.0000
1.2337
1.5527
1.9937
1.0000
1.2337
1.5527
1.9937
0.0000
0.0000
0.0000
0.0000
0.00
0.00
0.00
0.00
1.40
1.50
2.6116
3.4902
2.6117
3.4904
0.0001
0.0001
0.00
0.00
En la Tabla 6.6 comparan algunos resultados.
h = 0.1
h = 0.05
xn
Euler
Euler
mejorado
RK4
Valor
real
xn
Euler
Euler
mejorado
RK4
Valor
real
1.00
1.0000
1.0000
1.0000
1.0000
1.00
1.0000
1.0000
1.0000
1.0000
1.10
1.2000
1.2320
1.2337
1.2337
1.05
1.1000
1.1077
1.1079
1.1079
1.20
1.4640
1.5479
1.5527
1.5527
1.10
1.2155
1.2332
1.2337
1.2337
1.30
1.8154
1.9832
1.9937
1.9937
1.15
1.3492
1.3798
1.3806
1.3806
1.40
2.2874
2.5908
2.6116
2.6117
1.20
1.5044
1.5514
1.5527
1.5527
1.50
2.9278
3.4509
3.4902
3.4904
1.25
1.6849
1.7531
1.7551
1.7551
1.30
1.8955
1.9909
1.9937
1.9937
1.35
2.1419
2.2721
2.2762
2.2762
1.40
2.4311
2.6060
2.6117
2.6117
1.45
2.7714
3.0038
3.0117
3.0117
1.50
3.1733
3.4795
3.4903
3.4904
Errores de Truncamiento para el Método RK4
• Como es de grado 4, el error de truncamiento
local es O(h5) y el error de truncamiento
global es O(h4). Sin embargo, esto no se
abarca en este texto.
Ejemplo 2
Determine una cota para los errores de truncamiento
local del método RK4 aplicado a y '  2 xy, y (1)  1
Solución
Al calcular
la quinta derivada de la solución conocida
2
x 1
y( x)  e ,se obtiene
5
5
2
h
( 5)
3
5 c 1 h
(7)
y (c)  (120c  160c  32c )e
5!
5!
Así con c = 1.5, entonces (7) = 0.00028.
La Tabla 6.7 proporciona aproximaciones a la solución del
problema de valor inicial en x = 1.5 por el método RK4.
Tabla 6.7
h
0.1
Aproximación
3.49021064
Error
0.05
3.49033382
9.13776090  10-6
1.32321089  10-4
6.3 Métodos de Varios Pasos
• Método de Adams-Bashforth-Moulton
El predictor es la fórmula de Adams-Bashforth
h
 yn  (55 yn  59 yn 1  37n2  9 yn 3 ) ,
24
yn  f ( xn , yn )
yn 1  f ( xn1 , yn1 )
yn*1
yn 2  f ( xn2 , yn2 )
yn 3  f ( xn3 , yn3 )
donde n  3.
(1)
El valor de yn+1* se sustituye en el corrector de
Adams-Moulton
h
yn1  yn  (9 yn 1  19 yn  5 yn 1  yn 2 )
24
yn 1  f ( xn1 , yn*1 )
(2)
Ejemplo 1
Use el método anterior con h = 0.2 para obtener
y(0.8) para la solución de
y  x  y  1 ,
y (0)  1
Solución
Con h = 0.2, y(0.8) se aproxima mediante y4. En
principio s emplea el método RK4 con x0 = 0, y0 = 1,
h = 0.2 para obtener
y1 = 1.02140000, y2 = 1.09181796,
y3 = 1.22210646
Ejemplo 1 (2)
Ahora con x0 = 0, x1 = 0.2, x3 = 0.4, x4 = 0.6, y
f(x, y) = x + y – 1, hallamos
y0  f ( x0 , y0 )  (0)  (1)  1  0
y1  f ( x1 , y1 )  (0.2)  (1.02140000 )  1  0.22140000
y2  f ( x2 , y2 )  (0.4)  (1.09181796 )  1  0.49181796
y3  f ( x3 , y3 )  (0.6)  (1.22210646 )  1  0.82210646
El predictor (1) da
y4*
0.2
 y3 
(55 y3  59 y2  37 y1  9 y0 )  1.42535975
24
Ejemplo 1 (3)
Para usar el corrector (2), se necesita
y4  f ( x4 , y4* )  0.8  1.42535975  1  1.22535975
0.2
y4  y3 
(9 y4  19 y3  5 y2  y1 )  1.42552788
24
Estabilidad de Métodos Numéricos
• Decimos que un método numérico es estable,
si cambios pequeños en la condición inicial
dan como resultado sólo cambios pequeños
en la solución calculada.
6.4 Ecuaciones de Orden Superor y Sistemas
• PVI de Segundo Orden
Una PVI
y  f ( x, y, y) ,
y( x0 )  y0 ,
y( x0 )  u0
(1)
puede expresarse como
y  u
u   f ( x, y , u )
(2)
Como y’(x0) = u0, entonces y(x0) = y0, u(x0) = u0.
Aplicando el método de Euler (2)
yn1  yn  hun
un1  un  hf ( xn , yn , un )
(3)
Mientras que al aplicar el método RK4:
donde
1
yn1  yn  (m1  2m2  2m3  m4 )
6
1
un1  un  (k1  2k2  2k3  k4 )
6
(4)
m1  hun  hun
k1  hf ( xn , yn , un )
m2  h(un  1/2 k1 )
k2  hf ( xn  1/2 h, yn  1/2 m1 , un  1/2 k1 )
m3  h(un  1/2 k2 ) k3  hf ( xn  1/2 h, yn  1/2 m2 , un  1/2 k2 )
m4  h(un  k3 )
k4  hf ( xn  h, yn  m3 , un  k3 )
(n)
( n1)
y

f
(
x
,
y
,
y
'
,

,
y
)
En general,
Ejemplo 1
Use el método de Euler para obtener y(0.2),
donde
y  xy  y  0 ,
y (0)  1 ,
y(0)  2
Solución
Sea y’ = u, entonces (5) se transforma en
y  u
u   xu  y
De (3)
yn1  yn  hun
un1  un  h[ xnun  yn ]
(5)
Ejemplo 1 (2)
Usando h = 0.1, y0 = 1, u0 = 2, determinamos
y1  y0  (0.1)u0  1  (0.1)2  1.2
u1  u0  (0.1)[ x0u0  y0 ]  2  (0.1)[ (0)( 2)  1]  1.9
y2  y1  (0.1)u1  1.2  (0.1)(1.9)  1.3
u2  u1  (0.1)[  x1u1  y1 ]  1.9  (0.1)[ (0.1)(1.9)  1.2]  1.761
Fig 6.2
• En la Fig 6.2 se compara la curva solución generada
mediante el método de Euler con la curva solución
generada mediante el método RK4.
Ejemplo 2
Escribir
x  x  5 x  2 y  et
 2 x  y  2 y  3t 2
como un sistema de ecuaciones diferenciales de primer
orden.
Solución
Escribimos
x  2 y  et  5 x  x
y  3t 2  2 x  2 y
Al simplificar:
x  9 x  4 y  x  et  6t 2
Ejemplo 2 (2)
Sea
x'  u, y '  v
u  x  9 x  4 y  u  et  6t 2
v  y  2 x  2 y  3t 2
El sistema original se puede escribir en la forma
x  u
y  v
u  9 x  4 y  u  et  6t 2
v  2 x  2 y  3t 2
Solución Numérica de un Sistema
• La solución de un sistema de la forma
dx1
 f (1 t , x1 , x2 , , xn )
dt
dx2
 f 2 (t , x1 , x2 , , xn )
dt


dxn
 f n (t , x1 , x2 , , xn )
dt
se puede aproximar mediante métodos
numéricos.
Por ejemplo, mediante el método RK4:
x  f (t , x, y )
y  g (t , x, y )
x(t0 )  x0 ,
y(t0 )  y0
(6)
se parece a esto:
1
xn1  xn  xn  (m1  2m2  2m3  m4 )
6
1
yn1  yn  (k1  2k2  2k3  k4 )
6
(7)
donde
m1  hf (tn , xn , yn )
k1  hg (tn , xn , yn )
m2  hf (tn  1/2 h, xn  1/2 m1 , yn  1/2 k1 )
k2  hg (tn1/2 h, xn  1/2 m1 , yn  1/2 k1 )
m3  hf (tn  1/2 h, xn  1/2 m2 , yn  1/2 k2 ) k3  hg (tn  1/2 h, xn  1/2 m2 , yn  1/2 k2 )
m4  hf (tn  h, xn  m3 , yn  k3 )
k4  hg (tn  h, xn  m3 , yn  k3 )
(8)
Ejemplo 3
Considere
x  2 x  4 y
y   x  6 y
x(0)  1 ,
y (0)  6
Use el método RK4 para aproximar x(0.6) y y(0.6)
con h = 0.2 y h = 0.1.
Solución
Con h = 0.2 y los datos proporcionados, de (8)
Ejemplo 3 (2)
m1  hf (t0 , x0 , y0 )  0.2 f (0,  1, 6)  0.2[2(1)  4(6)]  4.4000
k1  hg (t0 , x0 , y0 )  0.2 g (0,  1, 6)  0.2[1(1)  6(6)]  7.4000
m2  hf (t0  1/2 h, x0  1/2 m1 , y0  1/2 k1 )  0.2 f (0.1, 1.2, 9.7)  8.2400
k2  hg (t0  1/2 h, x0  1/2 m1 , y0  1/2 k1 )  0.2 g (0.1, 1.2, 9.7)  11.4000
m3  hf (t0  1/2 h, x0  1/2 m2 , y0  1/2 k2 )  0.2 f (0.1, 3.12, 11.7)  10.6080
k3  hg (t0  1/2 h, x0  1/2 m2 , y0  1/2 k2 )  0.2 g (0.1, 3.12, 11.7)  13.4160
m4  hf (t0  h, x0  m3 , y0  k3 )  0.2 f (0.2, 8, 20.216)  19.3760
k4  hg (t0  h, x0  m3 , y0  k3 )  0.2 g (0.2, 8, 20.216)  21.3776
Ejemplo 3 (3)
Por lo tanto, de (7) obteenmos
1
x1  x0  (m1  2m2  2m3  m4 )
6
1
 1  (4.4  2(8.24)  2(10.608)  19.3760 )  9.2453
6
Observe Fig 6.3 y Tabla 6.8, 6.9.
Fig 6.3
Tabla 6.8
tn
xn
yn
0.00
0.20
-1.0000
9.2453
6.0000
19.0683
0.40
0.60
46.0327
158.9430
55.1203
150.8192
Tabla 6.9
tn
0.00
0.10
0.20
0.30
0.40
0.50
0.60
xn
-1.0000
2.3840
9.3379
22.5541
46.5103
88.5729
160.7563
yn
6.0000
10.8883
19.1332
32.8539
55.4420
93.3006
152.0025
6.5 Problemas de Valores en la Frontera de
Segundo Orden
• Aproximaciones por Diferencias Finitas
El desarrollo en serie de Taylor en a de y(x) es
xa
( x  a) 2
( x  a )3
y ( x)  y (a)  y(a)
 y(a)
 y(a)

1!
2!
3!
Si ponemos h = x – a, entonces
h
h2
h3
y ( x)  y (a)  y(a)  y(a)  y(a)  
1!
2!
3!
Escribiendo la última expresión como
h2
h3
y ( x  h)  y ( x)  y( x)h  y( x)  y( x)  
2
6
y
h2
h3
y ( x  h)  y ( x)  y( x)h  y( x)  y( x)  
2
6
(1)
(2)
Si h es pequeña, podemso despreciar y” y términos de
orden mayor, luego 1
y( x)  [ y ( x  h)  y ( x)
(3)
h
1

y ( x)  [ y ( x)  y ( x  h)]
(4)
h
Al restar (1) de (2) se obtiene también
1
(5)
y( x)  [ y ( x  h)  y ( x  h)]
2h
Si despreciamos los términos con h3 y superores,
entonces al sumar (1) y (2)
1
(6)
y( x)  2 [ y ( x  h)  2 y ( x)  y ( x  h)]
h
Los lados derechos de (3), (4), (5), (6) se
denominan cocientes de diferencias, y estas
diferencias se llaman diferencias finitas.
y(x + h) – y(x) : diferencia hacia delante
y(x) – y(x – h) : diferencia hacia atrás
y(x + h) – y(x – h): diferencia central
y(x + h) – 2y(x) + y(x – h): diferencia
central
Método de Diferencias Finitas
• Considere el PVF:
y  P( x) y  Q( x) y  f ( x) ,
y (a)   ,
y (b)   (7)
Supongase que a = x0 < x1 < … < xn < b representa
una partición regular del intervalo [a, b], esto es,
xi = a + ih, donde i = 0, 1, 2, ..., n, y h = (b – a)/n.
Estos puntos se llaman puntos de malla
interiores.
x1  a  h , x2  a  2h ,  , xn1  a  (n  1)h
Si permitimos que sea
yi = y(xi), Pi = P(xi), Qi = Q(xi), fi = f(xi),
y si y’ e y” en (7) se remplazan por (5) y (6),
entonces tenemos
yi 1  2 yi  yi 1
yi 1  yi 1
 Pi
 Qi yi  fi
2
2h
h
ó
1  h P  y  (2  h2Q ) y  1  h P  y  h2 f


i  i 1
i
i
i  i 1
i
2
2




Esto se conoce como ecuación de diferencias
finitas.
(8)
Ejemplo 1
Use (8) con n = 4 para aproximar la solución del
PVF
y  4 y  0 ,
y (0)  0 ,
y (1)  5
Solución
Tenemos P = 0, Q = –4, f(x) = 0, h = (1 – 0)/4 = ¼ .
De ahí que
yi 1  2.25 yi  yi 1  0
(9)
Los puntos interiores son x1= 0 + 1/4, x2= 0 + 2/4,
x3= 0 + 3/4, entonces (9) genera
Ejemplo (2)
y2  2.25 y1  y0  0
y3  2.25 y2  y1  0
y4  2.25 y3  y2  0
Junto con y0 = 0, y4 = 5, luego
 2.25 y1 
y2
y1  2.25 y2 
0
y3  0
y2  2.25 y3  5
Obtenemos y1 = 0.7256, y2 = 1.6327, y3 = 2.9479.
Ejemplo 2
Use (8) con n = 10 para aproximar la solución del
PVF
2
y  3 y  2 y  4 x ,
y (1)  1 ,
y (2)  6
Solución
Tenemos P = 3, Q = 2, f(x) = 4x2,
h = (2 – 1)/10 = 0.1, de ahí que (8) se transforma
1.15 yi 1  1.98 yi  0.85 yi 1 
2
0.04 xi
(10)
Los puntos interiores x1= 1.1, x2= 1.2, …, x9= 1.9,
y0 = 1, y10 = 6, luego (10) genera
Ejemplo 2 (2)
1.15 y2  1.98 y1  0.8016
1.15 y3  1.98 y2  0.85 y1  0.0576
1.15 y4  1.98 y3  0.85 y2  0.0676
1.15 y5  1.98 y4  0.85 y3  0.0784
1.15 y6  1.98 y5  0.85 y4  0.0900
1.15 y7  1.98 y6  0.85 y5  0.1024
1.15 y8  1.98 y7  0.85 y6  0.1156
1.15 y9  1.98 y8  0.85 y7  0.1296
1.98 y9  0.85 y8  6.7556
Podemos resolver este sistema de ecuaciones para
obtener y1, y2, …, y9.
Método de Disparos
• Otra manera de aproximar una solución se denomina
método de disparos. El punto de partida de este
método es remplazar el PVF por un PVI:
y  f ( x, y, y) ,
y (a)   ,
y(a)  m1
(11)
donde m1 es simplemente una suposición.
Esto se deja como ejercicio. Mirese el problema 14.