Libro de métodos-5 - MET PDF

Preview

Summary

Download Libro de métodos-5 - MET PDF

Description

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

Ing. Leonardo Flores González

1

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

INDICE

Introducción

3

Capítulo 1. Soluciones de Ecuaciones de una Variable

4

Método de la Bisección

5

Método del Punto Fijo

7

Método de Newton

10

Método de la Secante

12

Método de Bairstow

15

Capítulo 2. Sistema de Ecuaciones Lineales

17

Métodos Directos de Solución de Sistemas de Ecuaciones Lineales

18

Algoritmos para la Factorización LU

22

Método de Jacobi

28

Capítulo 3. Valores y Vectores Propios

32

Valores y Vectores Propios

33

Método de Jacobi Clásico

35

Método de la Potencia

37

Otros Algoritmos para Valores y Vectores Propios

41

Vibración Forzada Amortiguada – Sistema de Múltiples Grados de Libertad

45

Capítulo 4. Interpolación y Aproximación Polinomial

48

Interpolación de Funciones con Polinomios

49

Formas de Representar el Polinomio de Interpolación

50

Error de Interpolación Polinomial

50

Diferencias Divididas

51

Diferencias Finitas Centrales

54

Estabilidad del Método de Newmark

56

Capítulo 5. Integración Numérica

59

Cuadratura Cerrada de Newton-Cotes

60

Cuadratura de Gauss Legendre

62

Bibliografía Anexos

2

65 Anexo I

66

Anexo II

68

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

Introducción

El presente libro está dedicado a todos los estudiantes de la Facultad de Ingeniería Civil, especialmente a los alumnos del curso de Métodos Numéricos, quienes a través de sus consultas e inquietudes hicieron posible que se desarrollarán cada uno de los capítulos del presente libro.

3

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

Capítulo 1 SOLUCIONES DE ECUACIONES DE UNA VARIABLE

4

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

MÉTODO DE LA BISECCIÓN Dado f ( x)  0 y un intervalo a, b  , si la función tiene sólo una raíz r en a, b, f (a ) f (b)  0 y

f es continua a, b , el método converge a la raíz de f en el intervalo mencionado. Sea a n , bn  un intervalo donde existe una raíz de

f , si cn 

a n  bn 2

, f (a n ). f (cn )  0 entonces

bn1  cn , en caso contrario an 1  cn , esto nos indica lo siguiente:

bn  1  an  1 

1 2

(bn  an ) para todo n  0.

Además se puede ver que a0  a1   a n  b0 y a 0  bn   b1  b0 , de lo anterior podemos decir que:

bn  a n 

1 (b0  a0 ) . 2n

Si cn es una aproximación de la raíz, se tiene que el error absoluto cumple con la siguiente relación: r  cn 

1 1 bn  a n  n 1 b0  a 0 , tomando límites se tiene: 2 2

lim r  c n  0, lo que indica que la sucesión de c n , converge a la raíz.

n 

Una ecuación importante en el estudio de pavimentos es la ecuación AASHTO 93, esta sirve para calcular un parámetro llamado número estructural SN , en función de este parámetro se determina el espesor de las diferentes capas de un pavimento, a continuación describimos la solución de esta ecuación con ayuda del método de la bisección.

La ecuación AASHTO 93 se Indica a continuación:

 PSI  Log10   4.2  1.5 2.32. Log ( M ) 8.07 Log10 W18   Z R. S O  9.36.Log10( SN 1)  0.20  R 10 1094 0.40  ( SN  1) 5.19 Para la solución de esta ecuación Se emplean los siguientes datos:

5

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

1

2

3

4

5

6

7

8

9

10

11

CALCULO DEL NUMERO ESTRUCTURAL DE PAVIMENTOS FLEXIBLES Y ESPESOR DE LA LOSA DE PAVIMENTOS RIGIDOS Cd: coeficiente de drenaje, k: Módulo de reacción de la subrasante W 18: Tráfico estimado

W18

W 18

1.047E+06

ZR: Desviación estándar normal

ZR

ZR

-2.327

k

So: Error estándar combinado

So

So

0.35

Ec

 PSI MR

 PSI Sc

PSI: Pérdida de serviciabilidad MR: Módulo Resilente de la subrasante, Sc: Módulo de rotura del concreto Extremo izquierdo del intervalo de búsqueda del Número Estructutal o Espesor de Losa

SNi

Extremo derecho del intervalo de búsqueda del Número Estructutal o Espesor de Losa tol: Error Tolerancia Número Estructural Requerido o Espesor de Losa requerido

tol SNreq

Ejes Equivalentes

psi

2 440

Di

SNf

Cd psi

t

3.2 581

pci

3.120E+06

psi

1 2.5

5

Df tol D

Ejes Equivalentes J

12 0.0001 8.81

TABLA DE ITERACIONES Elija la opción adecuada pavimento flexible pavimento rígido

Iteraciones 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Di 8.5 8.5 8.5 8.5 8.71875 8.71875 8.7734375 8.80078125 8.80078125 8.807617188 8.807617188 8.807617188 8.807617188 8.808044434 8.808044434 8.808151245

Df

Dc f(Di) 12 8.5 -0.95860961 10.25 10.25 -0.08075805 9.375 9.375 -0.08075805 8.9375 8.9375 -0.08075805 8.9375 8.71875 -0.08075805 8.828125 8.828125 -0.02335606 8.828125 8.7734375 -0.02335606 8.828125 8.80078125 -0.00906982 8.814453125 8.814453125 -0.0019372 8.814453125 8.807617188 -0.0019372 8.811035156 8.811035156 -0.00015516 8.809326172 8.809326172 -0.00015516 8.80847168 8.80847168 -0.00015516 8.80847168 8.808044434 -0.00015516 8.808258057 8.808258057 -4.3793E-05 8.808258057 8.808151245 -4.3793E-05

f(Df) 0.76624159 0.76624159 0.363440518 0.146004485 0.033615534 0.033615534 0.005188251 0.005188251 0.005188251 0.001626434 0.001626434 0.000735696 0.000290284 6.75677E-05 6.75677E-05 1.18876E-05

f(Dc) Error -0.080758051 1.75 0.363440518 0.875 0.146004485 0.4375 0.033615534 0.21875 -0.02335606 0.109375 0.005188251 0.054688 -0.009069825 0.027344 -0.001937196 0.013672 0.001626434 0.006836 -0.000155156 0.003418 0.000735696 0.001709 0.000290284 0.000854 6.75677E-05 0.000427 -4.37931E-05 0.000214 1.18876E-05 0.000107 -1.59527E-05 5.34E-05

Los datos anteriores fueron procesados con visual Basic para Excel, en un programa que se entrega en el CD de programas, sin embargo un código similar en matlab es el siguiente:

function [res,i]=biseccion(a,b,tol) i=0; A = fopen('biseccion.xls','w'); while (b-a)/2>=tol x=(a+b)/2; d1=f(a); d2=f(x); i=i+1; y=[i a b d1 d2 x (b-a)/2]; if d1*d2tol it=it+1; x0=x1; x1=f(x0); E=abs(x1-x0); y=[it x0 x1 E ]; fprintf(A,'\t%d \t%6.7f \t%6.7f \t%6.7f\n',y); end fclose(A); res=x1;

9

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

Otro ejemplo común en Ingeniería Civil es encontrar la cantidad de acero necesaria para un elemento sometido a flexión, a continuación indicamos un ejemplo. Una viga soporta un momento de Mu  2596100 kgr  cm. , si la fórmula para encontrar la cantidad

a 2

de acero que debe tener la viga viene dada por Mu  0.9 As f y (d  )

y 0.85 f ' c ba  As f y ,

encontrar la cantidad de acero As , si:

d  44 cm., b  30 cm. , f 'c  350 kgr / cm 2 y f y  4200 kgr / cm 2 . Para resolver este problema notar que:

Mu

As 

0.9 f y ( d 

As f y 1.7bf ' c

)

Si iteramos considerando que: function [y]=f(x) y=2596100/(0.9*4200*(44-x*4200/(1.7*30*350))); Obtenemos: Iteraciones 1 2 3

10

As0

As1

6 16.1264945 17.0822006

16.1264945 17.0822006 17.1782798

Error 10.1264945 0.9557061 0.0960792

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

METODO DE NEWTON El método de la Newton es un método que aproxima a una función dada con un polinomio de primer grado, mediante el valor de la función y la derivada de la función en un entorno de un punto específico, es decir emplea una interpolación hermitiana. Este método es empleado para resolver

f ( x)  0 donde f :    , la sucesión generada por este método es:

xn1  xn 

f 'n ,  n  1 , donde fn representa a f ( xn ) . fn

Las condiciones de suficiencia para que el método de newton funcione son las siguientes: 1. f es continua en a ,b . 2. f ' es positiva o negativa en todo el intervalo 3. f ' ' es positiva o negativa en todo el intervalo

a ,b . a,b .

Desarrollando una función f por series de Taylor alrededor del punto xn se tiene:

f (x )  fn  f 'n ( x  xn ) 

1 ( x  xn ) 2 f ' ' ( ) 2

 es un punto del intervalo donde se busca la raíz. Supongamos por ejemplo que f y f 'son crecientes en a, b , si evaluamos a la función en la raíz se tiene que:

fn  f 'n (r  xn ) 

1 (r  xn )2 f ' ' ( )  0 . 2

Como la segunda derivada es positiva, se tiene que:

f n  f ' n ( r  x n )  0 , de esta última

expresión se puede ver fácilmente que:

r  xn 

fn  x n1 . f 'n

Esta última expresión indica que xn  r . Entonces una cota inferior de la sucesión xn es r . De xn1  xn 

f 'n ,  n  1 podemos ver que como f y f ' son positivas entonces fn

x n 1  x n , lo

que indica que tenemos una sucesión decreciente y acotada inferiormente, esto quiere decir que la sucesión converge a r . Todo lo anterior se puede resumir en el siguiente teorema: Teorema: Sea f (x) una función, con raíz única r en a , b . Si se satisfacen las tres condiciones de suficiencia indicadas anteriormente y se escoge un x 0 que satisface f 0 f ' ' 0  0, entonces la sucesión de puntos x0 , x1 ,  generada por el método de Newton, converge a

11

r. Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

La existencia del punto x0 es fácil de probar, ya que la función tendrá raíz única en el intervalo

a, b , por lo tanto existe una infinidad de valores positivos a la derecha o la izquierda de la raíz, según sea el caso. Programa en matlab del método de Newton function [res,it]=newton(x0,tol) it=1; A = fopen('newton.xls','w'); x1=x0-f(x0)/ f’(x0); E=abs(x1-x0); y=[it x0 f(x0) f’(x0) x1 E ]; fprintf(A,'\t%d \t%6.7f \t%6.7f \t%6.7f \t%6.7f

\t%6.7f\n',y);

while E>tol it=it+1; x0=x1; x1=x0-f(x0)/ f’(x0); E=abs(x1-x0); y=[it x0 f(x0) f’(x1) x1 E ]; fprintf(A,'\t%d \t%6.7f \t%6.7f \t%6.7f \t%6.7f

\t%6.7f\n',y);

end fclose(A); res=x1;

12

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

METODO DE LA SECANTE

El método de la secante es un método lagrangiano, se utiliza para encontrar los ceros de algunas funciones. La aproximación de la raíz de la ecuación

f ( x)  0 donde f :    , se calcula

conociendo dos aproximaciones x0 y x1 a partir de la sucesión:

x n1  x n  f n

x n  x n 1 ,  n  1 / f n  f n 1 , donde f n representa a f ( xn ) . fn  fn 1

Al aproximar una raíz de f con el método de la secante, se comete un error, para obtener una expresión del error cometido se utiliza el polinomio de interpolación que pasa por los puntos

xn y xn 1 , se representa a f (x) por medio del polinomio de interpolación más su expresión de error.

f ( x)  fn  ( x  xn )

f n  f n 1 1  (x  xn )(x  xn 1 ) f ' ' ( ) (1) donde  es un punto del intervalo x n  x n 1 2

donde se busca la raíz.

Del método de la secante tenemos:

0  fn  ( xn 1  xn )

fn  fn 1 (2) . Si reemplazamos la raíz r de f en (1) y luego restamos (2) xn  x n 1

obtenemos:

( r  xn 1 )

f n  f n 1 1  ( r  xn )( r  xn1 ) f ' ' ( )  0 (3) xn  xn  1 2

Como f tiene segunda derivada continua, en el intervalo en estudio, por el T.V.M. entonces existe

f ' ()  M 

f ' ' ( ) f n  f n 1 En En 1 si ,  xn , xn 1  , de lo anterior si En  r  xn se tiene que E n 1  xn  xn 1 2 f ' ( )

max f ' ' (x )  x I , donde I 2 min f ' ( x)

es el intervalo donde se busca la raíz, entonces

En 1  MEn En 1  ( 4) .

Teorema:

Sea

f ( x)  0

f C 2 a , b  f ' (x )  0  x  a , b ,

una

función,

entonces

si

con

raíz

x0  x1 V 1 (r ) ,

única la

r en a , b  sucesión

de

.

Si

puntos

M

x0 , x1 ,generada por el método de la secante, converge a

13

r.

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

Sea   max ME 0 , ME1  entonces si se cumple ME k  2    ME k 1   de (4) se cumple que

ME k   , finalmente se tiene que MEi   k i donde ki  ki 1  ki  2 , k0  k1  1 es la sucesión de Fibbonacci. Se debe tener presente que k i 

1 5

 

    1 5 2

i 1

 , como 0  lim ME  lim  k i  0 , se tiene i i  i  

i 1 5  1 2

que lim Ei  0 lo que indica que la sucesión converge r . i 

Suponga que el método de la secante converge, además se sabe que

C

f ' ' ( ) 2 f ' ()

ya que cuando n  

En 1  CE n En 1 , donde

  r    r , podemos determinar el orden de En1  KE np , de donde En  KE np1 , que

convergencia del método con la siguiente suposición 1

reemplazando en En 1  CE n En 1 resulta K E  CEn p

p n

1

1 p

, de donde K  C p  p 

1 (1  5) . 2

Nota: Sea f ( x )  0 y d  f ' ( x)  D,  x I  D  2d , entonces se puede obtener la raíz con un error absoluto menor que tol en el paso j si se cumple que x j  x j 1  tol.

3 2 Ejemplo: se sabe que la ecuación x  3x  x  6  0 tiene una sola raíz en el intervalo,

encontrar dicha raíz 1,1 .05 , considerando una tolerancia de 0.003:

Se tiene x 0  1 y x1  1.05 entonces x1  x0  0.05  tol , x2  x1  f1

x1  x0  1.097064 . Luego f1  f 0

se procede en forma similar para las demás iteraciones.

iteración

xn 1

f (x n 1 )

xn 1

xn

Error

1

1.097064

0.028081

1

1.05

0.05

2

1.094487

-0.000715

1.05

1.097064

0.047064

3

1.094551

-0.000001

1.097064

1.094487

0.002576

4

1.094551

-0.000001

1.094487

1.094551

0.0

Rpta: 1.094551

14

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

Programa en matlab del método de la secante. function [altura,teta,it]=fsecante(func,x0,x1,tol) it=0; d=feval(func,x1)*(x1-x0)/(feval(func,x1)-feval(func,x0)); while abs(d)>tol x2=x1-d; it=it+1; x0=x1; x1=x2; d=feval(func,x1)*(x1-x0)/(feval(func,x1)-feval(func,x0)); end respuesta=x1-d;it=it+1;

15

Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

METODO DE BAIRSTOW Sea P (x ) 

n

a x

k

k

k 0

un polinomio de grado n , si ak  R , k  0,1,..., n  , el método de Bairstow es

útil para encontrar los ceros complejos de P (x) . Teorema. Si P es un polinomio donde ak  R , k  0,1  , n y w es un cero complejo de P , entonces w es también un cero de P .

Teorema. Si el polinomio P (x ) 

n

a x k

k

se divide entre el factor cuadrático x 2  ux  v , entonces

k 0

el cociente Q ( x ) 

n

b x

k 2

k

y el residuo R ( x )  b1 ( x  u )  b0 , pueden calcularse de manera

k 2

recursiva, bk  a k  ubk 1  vbk  2 donde k  n, n  1 ..., 0 y bn 1  bn  2  0 .

Se tiene que:

P (x )  Q ( x )(x2  ux  v)  R( x), n

n

k 0

k 2

 ak x k   bk x k  2 (x 2  ux  v ) b1 (x  u )  b0

de donde bk  a k  ubk 1  vbk  2 , k  n, n  1 ..., 0 y bn 1  bn  2  0 .

2 De acuerdo a bk  a k  ubk 1  vbk  2 se tiene que bk : R  R es una función de (u , v ) , para que

R (x) sea cero bk (u , v )  0 para k  0,1 , de esta forma definimos c k 

bk b  dk  k 1 donde u v

k  0,1,2 ,..., n , si derivamos bk  a k  ubk 1  vbk  2 con respecto a c k  bk 1  uck 1  vc k  2 ,

definimos

de

esta

manera

la

relación

u de

se tiene que recursividad

c k  bk 1  uck 1  vck  2 donde cn 1  cn  0 , análogamente derivando con respecto v se tiene que d k  bk 1  ud k 1  vd k  2 donde d n 1  dn  0 , de estas relaciones de recursividad se puede ver que ck  d k . Linealizando bk (u  u , v  v)  0 se tiene bk (u , v ) 

bk b  u  k v  0 para u v

k  0,1 resultando:

c0 c1 16

c1  u b0 c 0 c1 donde J  .  c2  v b1 c1 c 2 Leonardo Flores González

UNIVERSIDAD NACIONAL DE INGENIERÍA FACULTAD DE INGENIERÍA CIVIL

Teorema.: Sean u 0 ,v0 valores reales tales que los ceros de x 2  u0 x  v0 son ceros simples, entonces J  0 . Algoritmo: 1. Indicar valores iniciales de (u0 , v0 ) y   0 . 2. Para k  n , n  1,... 0 hacer bk  a k  ubk 1  vbk  2 con bn 1  bn  2  0 . 3. Para k  n , n  1,... 0 hacer ck  bk  uck 1  vc k  2 con cn1  cn  2  0 . 4. Resolver

5. Hacer

c0 c1

c1 u b0  . c2 v b1

u u u .   v v v

6. Repetir (2,3,4,5,6) hasta u  ...