• Author / Uploaded
• alejandro silva

Citation preview

TALLER MÉTODOS NUMÉRICOS

Ing. Matilde Montealegre Madero, MSc.

Fecha de envío : Mayo 24 de 2018

Preparacion Parcial 3, Lunes 28 de Mayo de 2018.

Temas para el parcial 3. 1. Solucion de sistemas de ecuaciones algebraicas lineales por los metodos de: Jacobi, Gauss Seidel, Relajacion, Algoritmo de Thomas. 2. Ecuaciones Diferenciales Ordinarias. 3. Programacion en Macros en excel o el lenguaje de su preferencia. ECUACIONES DIFERENCIALES ORDINARIAS Una ecuación diferencial es aquella que involucra una función desconocida y sus derivadas1. En una ecuación diferencial ordinaria (EDO) la función desconocida depende de una sola variable; ésta contiene una variable dependiente y una variable independiente únicamente. Se clasifican por su orden (p), el cual depende del valor más alto de la derivada. Una ecuación de primer orden es de la forma: dv c =g− v dt m

donde la variable dependiente es v y la variable independiente es t. Una ecuación de segundo orden es: d2x dx m 2 + c + kx = 0 dt dt

(7.1)

Las ecuaciones de más alto orden pueden reducirse a un sistema de ecuaciones de primer orden definiendo una nueva variable. Ejemplo: si

m

y=

dx dt

entonces

dy d 2 x = dt dt 2

sustituyendo

en

la

Ec.

7.1

resulta

dy + cy + kx = 0 dt

de donde:

dy cy + kx =− dt m

Una solución de una ecuación diferencial ordinaria es una función específica de la variable independiente y los parámetros que satisfacen la ecuación diferencial original. Por ejemplo, dada la ecuación:

y = 2 x 3 + 5 x 2 − x + 12

(7.2)

diferenciando la Ec. 7.2, se tiene una EDO dada por: dy = 6 x 2 + 10 x − 1 dx

(7.3)

Luego, para determinar la ecuación original dada la EDO se integra la Ec. 7.3 y =  ( 6 x 2 + 10 x − 1)dx y = 2 x3 + 5 x 2 − x + c

(7.4)

c es la constante de integración e indica que la solución no es única. Luego, un número infinito de posibles valores de c satisfacen la ecuación diferencial. Por tanto, para especificar la solución completamente, las ecuaciones diferenciales se acompañan de condiciónes auxiliares. Para una ecuación de primer orden se requiere una condición auxiliar llamada valor inicial o condición inicial para determinar la constante c y obtener una solución única. Para el ejemplo bajo consideración, una condición inicial podría ser que para x = 0, y = 1 , que sustituido en la Ec. 7.4 permite hallar el valor de c = 1, que es el valor de la Ec. 7.2. Cuando se trata de ecuaciones diferenciales de n-ésimo orden, se requieren n condiciónes para obtener una solución única. Los métodos numéricos consideran la ecuación diferencial: y'=

dy = f ( x, y ) dx

que satisfacen la condición inicial y ( x0 ) = y0 . Estos métodos aproximan la función solución y (t ) en un número discreto de valores de las variables, es decir generan una secuencia de valores de y (t ) : y1 , y2 , y3 ..., en los puntos x1 , x2 , x3 .... Algunos métodos para evaluar numéricamente ecuaciones diferenciales se presentan a continuación.

7.1. MÉTODO DE EULER1-7 Este método se utiliza para resolver EDOs de la forma dy = f ( x, y ) dx para el cual se obtiene un nuevo estimado de yi +1 en función del anterior yi . Por tanto, yi +1 = yi + φ h

(7.5)

El método de Euler usa una ecuación diferencial para estimar la pendiente en la forma de la primera derivada en el punto xi . Luego la pendiente en el punto xi es φ = f ( xi , yi ) , donde f ( xi , yi ) es la ecuación diferencial calculada en el punto xi , yi . Substituyendo en la Ec. 7.5 se tiene: yi +1 = yi + f ( xi , yi )h

(7.6)

La cual se conoce como el método de Euler (ó método de Euler-Cauchy).

Option Explicit Sub eulercauchy() Dim i As Double, posx As Integer Dim limiteinf As Double, limitesup As Double, h As Double Dim xi As Double, yi As Double, yim1 As Double, iteracion As Integer Dim valorverdadero As Double, elerror As Double iteracion = 0 xi = Cells(4, 3).Value yi = Cells(5, 3).Value h = Cells(6, 3).Value limiteinf = Cells(7, 3).Value limitesup = Cells(8, 3).Value yim1 = yi posx = 12 For i = limiteinf To limitesup Step h Cells(posx, 2).Value = iteracion Cells(posx, 4).Value = xi Cells(posx, 5).Value = yi Cells(posx, 6).Value = fEDO(xi, yi) Cells(posx, 7) = yim1 valorverdadero = frealEDO(xi, yi) Cells(posx, 3) = valorverdadero If yim1 0 Then elerror = Abs((valorverdadero - yim1) / valorverdadero) * 100 Cells(posx, 8) = elerror yim1 = yi + fEDO(xi, yi) * h 'actualiza valores de xi,yi para siguiente iteración xi = xi + h yi = yim1 posx = posx + 1 iteracion = iteracion + 1 Next i End Sub

Function frealEDO(x, y) frealEDO = -0.5 * x ^ 4 + 4 * x ^ 3 - 10 * x ^ 2 + 8.5 * x + 1 End Function

Function fEDO(x, y) fEDO = -2 * x ^ 3 + 12 * x ^ 2 - 20 * x + 8.5 End Function EJEMPLO 1 Use el método de Euler para integrar numéricamente f ( x, y ) = −2 x3 + 12 x 2 − 20 x + 8.5 para x=0 hasta x=4 con intervalos, h, de 0.5. La condición inicial en x=0 es y=1. La solución exacta está dada por y = −0.5 x 4 + 4 x3 − 10 x 2 + 8.5 x + 1 . (Tomado de la Ref. 1) SOLUCIÓN A partir de x0 = 0 , y0 = 1 calcule xi +1 = xi + h y yi +1 = yi + f ( xi , yi )h .

Iteración

0:

i=0;

calcule

x1 = x0 + h ,

entonces,

x1 = 0 + 0.5 = 0.5 ,

y1 = y0 + f ( x0 , y0 )h

y1 = 1 + f (0,1)0.5 = 1 + [−2(0)3 + 12(0) 2 − 20(0) + 8.5]0.5 = 1 + 8.5* 0.5 = 5.25

Iteración 1: i=1; calcule x2 = x1 + h , luego, x2 = 0.5 + 0.5 = 1.0 , y2 = y1 + f ( x1 , y1 )h

y2 = 5.25 + f (0.5,5.25)0.5 = 5.25 + [−2(0.5)3 + 12(0.5) 2 − 20(0.5) + 8.5]0.5 y2 = 5.25 + 1.25*0.5 = 5.875

La siguente tabla muestra los resultados obtenidos hasta x=4. La solución verdadera para x=1.0 es 3.0, sin embargo el porcentaje de error es bastante alto (95,8%). En la Fig. 1 se puede observar que los cálculos capturan el comportamiento de la función pero con un error considerablemente alto. En la tabla 1 se muestra los cálculos disminuyendo el valor de h a 0.25, lo que demuestra según la Fig. 1 que el error disminuye. Luego, si se calcula con incrementos bastante pequeños se puede obtener mejores resultados. h=0.5

iteración 0 1 2 3 4 5 6 7 8

yreal

1 3.21875 3 2.21875 2 2.71875 4 4.71875 3

xi

0 0.5 1 1.5 2 2.5 3 3.5 4

yi

1 5.25 5.875 5.125 4.5 4.75 5.875 7.125 7

f(xi,yi)

yeuler {y(xi)}

8.5 1.25 -1.5 -1.25 0.5 2.25 2.5 -0.25 -7.5

1 5.25 5.875 5.125 4.5 4.75 5.875 7.125 7

error % 0 63.1067 95.8333 130.9859 125 74.7126 46.875 50.993 133.333

h=0.25 yeuler {y(xi)} error %

1.00 4.18 4.34 3.55 3.13 3.62 4.84 5.87 5.00

0.00 29.85 44.79 60.21 56.25 33.05 21.09 24.34 66.67

8 7 6 5 4 3 2 Curva real h=0.5

1

h=0.25 0 0

0.5

1

1.5

2

2.5

3

3.5

4

Fig. 1. Resultados de integración por el método de Euler para el ejemplo 1

7.2. MÉTODOS DE RUNGE-KUTTA1,2,7 Estos están diseñados para aproximar los métodos basados en las series de Taylor. Tienen la ventaja de que no requieren una evaluación explicita de las derivadas de f(x,y). La idea básica es la de utilizar una combinación lineal de los valores de f(x,y) para aproximar y(x). La combinación lineal se toma de manera que se adapte lo mejor posible al desarrollo de la serie de Taylor para obtener métodos con un orden p elevado. Existen muchas variaciones de los método de Runge-Kutta pero todos pueden ser obtenidos de la generalización de la Ec. 7.6: yi +1 = yi + φ ( xi , yi , n)h

(7.7)

donde φ(xi, yi, n) es llamada la función incremento la cual puede interpretarse como una pendiente representativa sobre el intervalo. En forma general puede escribirse como:

φ = a1k1 + a2 k2 + ... + an kn donde las a ' s son constantes y las k‘s están dadas por: k1 = f ( xi , yi ) k2 = f ( xi + p1h, yi + q11k1h) k3 = f ( xi + p2 h, yi + q21k1h + q22 k2 h)

kn = f ( xi + pn −1h, yi + qn −1,1k1h + qn −1,2 k2 h + ... + qn −1,n −1kn −1h) Varios tipos de métodos de Runge-Kutta pueden utilizarse empleando diferentes términos dados por la función incremento especificada por n. Así, el método de Runge-Kutta de primer orden con n=1, es el método de Euler. Una vez se escoge n, se evalúan los valores de a’s, p’s y q’s igualando los términos de la Ec. 7.7 con la expansión de las series de Taylor.

7.2.1. Método de Runge-Kutta de orden 1 (Método de Euler), n=1. De la Ec. 7.7 yi +1 = yi + φ ( xi , yi , n = 1)h , φ = a1k1 Luego: yi +1 = yi + (a1k1 )h donde, k1 = f ( xi , yi ) Reemplazando yi +1 = yi + a1 f ( xi , yi )h , si a1 = 1 , se tiene el método de Euler: yi +1 = yi + f ( xi , yi )h

7.2.2. Método de Runge-Kutta de orden 2, n=2 Para este caso particular se tiene: yi +1 = yi + φ ( xi , yi , n = 2)h , φ = a1k1 + a2 k2 Luego; yi +1 = yi + (a1k1 + a2 k2 )h

(7.8)

donde: k1 = f ( xi , yi ) , y, k2 = f ( xi + p1h, yi + q11k1h) El siguiente código calcula por el método de Runge Kutta de orden 2. Los valores de las a’s, p’s y q’s deben ser inicializados de acuerdo al método que se quiera utilizar.

Sub rungekuttan2()

Dim i As Double, posx As Integer Dim limiteinf As Double, limitesup As Double, h As Double Dim xi As Double, yi As Double, yim1 As Double, iteracion As Integer Dim valorverdadero As Double, elerror As Double Dim a1 As Double, a2 As Double, p1 As Double, q11 As Double Dim k1 As Double, k2 As Double 'inicialización de paramétros Método de Runge Kutta con n=2. Aquí valores de Heun a1 = 0.5 a2 = 0.5 p1 = 1 q11 = 1 Range("B12:I100").ClearContents iteracion = 0 xi = Cells(4, 3).Value yi = Cells(5, 3).Value h = Cells(6, 3).Value limiteinf = Cells(7, 3).Value limitesup = Cells(8, 3).Value yim1 = yi posx = 12 For i = limiteinf To limitesup Step h Cells(posx, 2).Value = iteracion Cells(posx, 4).Value = xi Cells(posx, 5).Value = yi k1 = fEDO(xi, yi) k2 = fEDO(xi + p1 * h, yi + q11 * k1 * h) Cells(posx, 6).Value = k1 Cells(posx, 7).Value = k2 Cells(posx, 8) = yim1 valorverdadero = frealEDO(xi, yi) Cells(posx, 3) = valorverdadero If yim1 0 Then elerror = Abs((valorverdadero - yim1) / valorverdadero) * 100 Cells(posx, 9) = elerror yim1 = yi + (a1 * k1 + a2 * k2) * h 'actualiza valores de xi,yi para siguiente iteración xi = xi + h yi = yim1 posx = posx + 1 iteracion = iteracion + 1 Next i End Sub Ralston1, por ejemplo, determinó que a1 = 1/3 y a2 = 2/3, aunque estas variables podrían tomar otros valores. A continuación se presentan tres versiones de las variables a’s.

1) Método de Heun Toma los valores de a2 = 1/ 2 , a1 = 1/ 2 y p1 = q11 = 1 , luego al reemplazar en la Ec. 7.8, se obtiene la ecuación de Runge-Kutta de segundo orden (RK2):

1  1 yi +1 = yi +  k1 + k2  h 2  2 Donde: k1 = f ( xi , yi )

y

k2 = f ( xi + h, yi + k1h)

Se puede observar en los resultados obtenidos para el ejemplo 1, que utilizando el método de Heun, el error disminuye considerablemente. h=0.5 yreal

iteración 0 1 2 3 4 5 6 7 8

1.00 3.22 3.00 2.22 2.00 2.72 4.00 4.72 3.00

xi

0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00

yi

1.00 3.44 3.38 2.69 2.50 3.19 4.38 4.94 3.00

k1

8.50 1.25 -1.50 -1.25 0.50 2.25 2.50 -0.25 -7.50

k2

1.25 -1.50 -1.25 0.50 2.25 2.50 -0.25 -7.50 -20.75

yheun {y(xi)}

1.00 3.44 3.38 2.69 2.50 3.19 4.38 4.94 3.00

error % 0.00 6.80 12.50 21.13 25.00 17.24 9.38 4.64 0.00

h=0.25 yheun {y(xi)} error

% 1.00 3.27 3.09 2.34 2.13 2.84 4.09 4.77 3.00

0.00 1.70 3.13 5.28 6.25 4.31 2.34 1.16 0.00

2) Método de Euler Modificado Este asume que a2 = 1 , a1 = 0

y

p1 = q11 = 1/ 2 , luego sustituyendo en la Ec. 7.8, se tiene:

yi +1 = yi + k2 h donde; 1 1   k2 = f  xi + h, yi + k1h  2 2  

3. Método de Ralston Ralston toma a2 = 2 / 3 , a1 = 1/ 3 y 2  1 yi +1 = yi +  k1 + k2  h 3  3

p1 = q11 = 3 / 4 . Luego sustituyendo en la Ec. 7.8;

donde; k1 = f ( xi , yi )

y

3 3   k2 = f  xi + h, yi + k1h  4 4  

7.2.3. Métodos de Runge-Kutta de Mayor Orden La versión común de tercer orden (RK3) es: 1  yi +1 = yi +  (k1 + 4k2 + k3 )  h 6  donde, k1 = f ( xi , yi ) 1 1   k2 = f  xi + h, yi + k1h  2 2   k3 = f ( xi + h, yi − k1h + 2k2 h)

Sub rungekuttan3() Dim i As Double, posx As Integer Dim limiteinf As Double, limitesup As Double, h As Double Dim xi As Double, yi As Double, yim1 As Double, iteracion As Integer Dim valorverdadero As Double, elerror As Double Dim k1 As Double, k2 As Double, k3 As Double Range("B12:j100").ClearContents iteracion = 0 xi = Cells(4, 3).Value yi = Cells(5, 3).Value h = Cells(6, 3).Value limiteinf = Cells(7, 3).Value limitesup = Cells(8, 3).Value yim1 = yi posx = 12 For i = limiteinf To limitesup Step h Cells(posx, 2).Value = iteracion Cells(posx, 4).Value = xi Cells(posx, 5).Value = yi k1 = fEDO(xi, yi) k2 = fEDO(xi + (1 / 2) * h, yi + (1 / 2) * k1 * h) k3 = fEDO(xi + h, yi - k1 * h + 2 * k2 * h) Cells(posx, 6).Value = k1 Cells(posx, 7).Value = k2

Cells(posx, 8).Value = k3 Cells(posx, 9) = yim1 valorverdadero = frealEDO(xi, yi) Cells(posx, 3) = valorverdadero If yim1 0 Then elerror = Abs((valorverdadero - yim1) / valorverdadero) * 100 Cells(posx, 10) = elerror yim1 = yi + (1 / 6 * (k1 + 4 * k2 + k3)) * h 'actualiza valores de xi,yi para siguiente iteración xi = xi + h yi = yim1 posx = posx + 1 iteracion = iteracion + 1 Next i End Sub

En el ejemplo 1, dado que la ecuación diferencial es una cúbica, el método de RK3 se ajusta perfectamente y por tanto el error es cero. h=0.5

iteración 0 1 2 3 4 5 6 7 8

yreal

xi

yi

k1

k2

k3

1.00 3.22 3.00 2.22 2.00 2.72 4.00 4.72 3.00

0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00

1.00 3.22 3.00 2.22 2.00 2.72 4.00 4.72 3.00

8.50 1.25 -1.50 -1.25 0.50 2.25 2.50 -0.25 -7.50

4.22 -0.59 -1.66 -0.47 1.47 2.66 1.59 -3.22 -13.28

1.25 -1.50 -1.25 0.50 2.25 2.50 -0.25 -7.50 -20.75

La versión clásica de cuarto orden (RK4) es: 1  yi +1 = yi +  (k1 + 2k2 + 2k3 + k4 )  h 6  donde: k1 = f ( xi , yi )

yRK3 {y(xi)}

1.00 3.22 3.00 2.22 2.00 2.72 4.00 4.72 3.00

error % 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1 1   k2 = f  xi + h, yi + k1h  2 2   1 1   k3 = f  xi + h, yi + k2 h  2 2   k4 = f ( xi + h, yi + k3 h) Cuando se requiere mayor precisión se recomienda el método de Runge-Kutta de orden 5 (RK5 o método de Butcher): 1  yi +1 = yi +  (7 k1 + 32k3 + 12k4 + 32k5 + 7 k6 )  h  90  Siendo; k1 = f ( xi , yi ) 1 1   k2 = f  xi + h, yi + k1h  4 4   1 1 1   k3 = f  xi + h, yi + k1h + k2 h  4 8 8   1 1   k4 = f  xi + h, yi − k2 h + k3 h  2 2   3 3 9   k5 = f  xi + h, yi + k1h + k4 h  4 16 16   3 2 12 12 8   k6 = f  xi + h, yi − k1h + k2 h + k3 h − k4 h + k5 h  7 7 7 7 7  

EJERCICIOS PROPUESTOS 1. Use el método de Euler para integrar numéricamente f ( x, y ) = x y para x=0 hasta x=1 con intervalos (h) de 0.25 y 0.5. La condición inicial en x=0 es y=1. La solución exacta está dada por ( x 2 + 4) 2 y= . Grafique y compare los resultados con la solución exacta. 16 a) Use el método de Heun con h=0.5 para el resolver el problema 1. b) Use el método de Runge-Kutta de orden 3 con h=0.5 para resolver el problema 1. c) Use el método de Runge-Kutta de orden 4 con h=0.5 para resolver el problema 1. d) Use el método de Runge-Kutta de orden 5 con h=0.5 para resolver el problema 1.

2

2. Utilice los métodos de Runge-Kutta para integrar y ' = y + 1 , para x = 0 hasta x = 2 con diferentes intervalos de h. Las condiciónes iniciales para x = 0 es y = 0 . Modifique el valor de h para mejores resultados. Grafique la solución. Compare los resultados con el resultado analítico

y = tan x

y los obtenidos con el método de Euler. 2

3. Utilice Runge-Kutta de orden 5 integrar y ' = yx − y , para x = 0 hasta x = 2 con diferentes intervalos de h. Las condiciónes iniciales para x = 0 es y = 1 . Modifique el valor de h para mejores resultados. Grafique la solución. Compare los resultados con el resultado analítico obtenidos con el método de Euler.

y=e

x3 −x 3

y los

0.8 x

− 0.5 y , para x = 0 hasta x = 4 4. Utilice Runge-Kutta de orden 5 para integrar y ' = 4e con intervalos de h=1. Las condiciones iniciales para x = 0 es y = 2 . Grafique. Compare los resultados con el resultado analítico

y=

4 0.8 x −0.5 x e −e + 2e −0.5 x y los obtenidos con el ( ) 1.3

método de Euler.

5. Utilice los métodos de Runge-Kutta para integrar

x2 + y2 y'= xy

, para x = 0 hasta x = 4 con

intervalos de h=1. Las condiciones iniciales para x = 1 es y = 3 . Grafique. Compare los resultados con el resultado analítico

y = x 9 + ln x 2

y los obtenidos con el método de Euler.