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 65
Anexos Anexo I 66
Anexo II 68

2 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 (c n )  0 entonces

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

1
bn 1  an 1  (bn  an ) para todo n  0.
2
Además se puede ver que a 0  a1   a n  b0 y a 0  bn   b1  b0 , de lo anterior podemos decir

que:
1
bn  a n  (b0  a0 ) .
2n
Si c n es una aproximación de la raíz, se tiene que el error absoluto cumple con la siguiente

relación:
1 1
r  cn  bn  an  n1 b0  a0 , 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 
Log10 W18   Z R .SO  9.36.Log10 (SN  1)  0.20   4.2  1.5   [Link] (M )  8.07
10 R
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 W 18 Ejes Equivalentes W 18 1.047E+06 Ejes Equivalentes J 3.2
ZR: Desviación estándar normal ZR ZR -2.327 k 581 pci
So: Error estándar combinado So So 0.35 Ec 3.120E+06 psi
PSI: Pérdida de serviciabilidad  PSI  PSI 2 Cd 1
MR: Módulo Resilente de la subrasante, Sc: Módulo de rotura del concreto MR psi Sc 440 psi t 2.5
Extremo izquierdo del intervalo de búsqueda del Número Estructutal o Espesor de Losa SNi Di 5
Extremo derecho del intervalo de búsqueda del Número Estructutal o Espesor de Losa SNf Df 12
tol: Error Tolerancia tol tol 0.0001
Número Estructural Requerido o Espesor de Losa requerido SNreq D 8.81

TABLA DE ITERACIONES
Elija la opción adecuada

pavimento flexible Iteraciones Di Df Dc f(Di) f(Df) f(Dc) Error

1 8.5 12 8.5 -0.95860961 0.76624159 -0.080758051 1.75
pavimento rígido 2 8.5 10.25 10.25 -0.08075805 0.76624159 0.363440518 0.875
3 8.5 9.375 9.375 -0.08075805 0.363440518 0.146004485 0.4375
4 8.5 8.9375 8.9375 -0.08075805 0.146004485 0.033615534 0.21875
5 8.71875 8.9375 8.71875 -0.08075805 0.033615534 -0.02335606 0.109375
6 8.71875 8.828125 8.828125 -0.02335606 0.033615534 0.005188251 0.054688
7 8.7734375 8.828125 8.7734375 -0.02335606 0.005188251 -0.009069825 0.027344
8 8.80078125 8.828125 8.80078125 -0.00906982 0.005188251 -0.001937196 0.013672
9 8.80078125 8.814453125 8.814453125 -0.0019372 0.005188251 0.001626434 0.006836
10 8.807617188 8.814453125 8.807617188 -0.0019372 0.001626434 -0.000155156 0.003418
11 8.807617188 8.811035156 8.811035156 -0.00015516 0.001626434 0.000735696 0.001709
12 8.807617188 8.809326172 8.809326172 -0.00015516 0.000735696 0.000290284 0.000854
13 8.807617188 8.80847168 8.80847168 -0.00015516 0.000290284 6.75677E-05 0.000427
14 8.808044434 8.80847168 8.808044434 -0.00015516 6.75677E-05 -4.37931E-05 0.000214
15 8.808044434 8.808258057 8.808258057 -4.3793E-05 6.75677E-05 1.18876E-05 0.000107
16 8.808151245 8.808258057 8.808151245 -4.3793E-05 1.18876E-05 -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('[Link]','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*d2<0
b=x;
else
a=x;
end
fprintf(A,'\t%d\t%6.7f\t%6.7f\t%6.7f\t%6.7f\t%6.7f\t %6.7f\n',y);
end
res=(a+b)/2;
fclose(A);

f es la función cuya raíz se desea encontrar.

6 Leonardo Flores González

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

Encontrar la raíz cuadrada de cinco con una tolerancia de 0.0001

En matlab la función f es:
function [y]=f(x)
y=x^2-5;

Iteraciones a b f(a) f(c) c Error

1 1.000000 4.000000 -4.000000 1.250000 2.500000 1.500000
2 1.000000 2.500000 -4.000000 -1.937500 1.750000 0.750000
3 1.750000 2.500000 -1.937500 -0.484375 2.125000 0.375000
4 2.125000 2.500000 -0.484375 0.347656 2.312500 0.187500
5 2.125000 2.312500 -0.484375 -0.077148 2.218750 0.093750
6 2.218750 2.312500 -0.077148 0.133057 2.265625 0.046875
7 2.218750 2.265625 -0.077148 0.027405 2.242188 0.023438
8 2.218750 2.242188 -0.077148 -0.025009 2.230469 0.011719
9 2.230469 2.242188 -0.025009 0.001164 2.236328 0.005859
10 2.230469 2.236328 -0.025009 -0.011931 2.233398 0.002930
11 2.233398 2.236328 -0.011931 -0.005386 2.234863 0.001465
12 2.234863 2.236328 -0.005386 -0.002112 2.235596 0.000732
13 2.235596 2.236328 -0.002112 -0.000474 2.235962 0.000366
14 2.235962 2.236328 -0.000474 0.000345 2.236145 0.000183

7 Leonardo Flores González

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

METODO DEL PUNTO FIJO

Un punto fijo de una función f , es un p Domf tal que p  f ( p ). Este método es útil para
calcular los ceros de algunas funciones.

Teorema
1. Sea f : a, b  c, d  / c, d   a, b , entonces f tiene un punto fijo en a,b .

2. Si f ' existe en a, b y existe una constante positiva k  1 / f ' ( x)  k  x  a, b , entonces el

punto fijo es único.

3. Si se cumple 1 y 2, entonces para cualquier p 0 a, b , la sucesión generada por

p n  f ( p n 1 ),  n  1 , converge al único punto fijo p a, b  .

Se puede considerar que si f (a )  a o f (b)  b , entonces el punto fijo existe, para los demás

casos, sea h( x)  f ( x)  x , entonces h( a )  0 y h(b)  0 , como h(x) es continua en a, b  , por el

teorema del valor intermedio existe un p  a, b / h( p )  0.

Para probar la unicidad, supongamos que p y q son puntos fijos de f , entonces por el teorema
del valor medio existe un  entre p y q tal que:

f ( p)  f (q)  p  q  f ' ( ) p  q  k p  q  p  q .

Lo que indica que 1  k , esto es una contradicción, por tal razón el punto fijo es único.
Para probar (3), por el teorema del valor medio tenemos que:

p n  p  f ( p n 1 )  f ( p)  f ' ( ) p n 1  p  k p n 1  p ,  es un punto del domino de f . De

esta última ecuación se tiene que p n  p  k n p 0  p . Si tomamos límites cuando n   ,

entonces se aprecia que p n converge a p.

Existen algunas expresiones para acotar el error como:

kn
pn  p  p1  p0
1 k

8 Leonardo Flores González

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

Un ejemplo que ilustra lo mencionado anteriormente es probar que la sucesión generada por
x n 1 A
xn    n  1 , converge a A.
2 2 x n 1

x A
Se sabe que si 0  x  A , entonces ( x  A ) 2  0 , esto indica que   A , lo que indica
2 2x
que f ( x)  A , si escogemos un x  A , entonces todos f ( x)  A ya que ( x  A)2  0;

sea f ( x) 
x A

2 2x
una función con dominio  
A , x 0 tal que x 0  A , se puede ver que

cualquier punto del domino cumple con x2  A esta última expresión equivale a
x A
x   f ( x) , como x0  x  f ( x)  A , entonces existe punto fijo; además
2 2x
1 A 1
f ' ( x)  1  2  , se puede ver que f ' ( x)  , esto indica que el punto fijo es único, finalmente
2 x  2

 x  
A , x0 , la sucesión generada por x n 
x n 1
2

A
2 x n 1
 n  1 converge al punto fijo, además

se verifica fácilmente que f ( A )  A.

El algoritmo en matlab del punto fijo es el siguiente:

function [res,it]=pfijo(x0,tol)
it=0;
A = fopen('[Link]','w');
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);
while abs(x0-x1)>tol
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
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 ,
2
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 
As f y
0.9 f y (d  )
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 As0 As1 Error

1 6 16.1264945 10.1264945
2 16.1264945 17.0822006 0.9557061
3 17.0822006 17.1782798 0.0960792

10 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:
f 'n
xn 1  xn  ,  n  1 , donde f n 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 a,b .

3. f ' ' es positiva o negativa en todo el intervalo a,b .

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

1
f ( x)  f n  f 'n ( x  xn )  ( 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:
1
f n  f 'n (r  xn )  (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:

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

f 'n
De xn 1  xn  ,  n  1 podemos ver que como f y f ' son positivas entonces x n 1  x n , lo
fn
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 r.

11 Leonardo Flores González

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

La existencia del punto x 0 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('[Link]','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:

xn  xn 1
xn 1  xn  f n ,  n  1 / f n  f n 1 , donde f n representa a f ( xn ) .
f n  f n 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 n  f n 1 1
f ( x)  f n  ( x  xn )  ( x  xn )( x  xn 1 ) f ' ' ( )(1) donde  es un punto del intervalo
xn  xn 1 2
donde se busca la raíz.

Del método de la secante tenemos:

f n  f n 1
0  f n  ( xn 1  xn ) (2) . Si reemplazamos la raíz r de f en (1) y luego restamos (2)
xn  xn 1
obtenemos:
f n  f n 1 1
(r  xn 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 n  f n 1 f ' ' ( )
f ' ( )  ,  xn , xn 1  , de lo anterior si En  r  xn se tiene que En 1  En En 1 si
xn  xn 1 2 f ' ( )

max f ' ' ( x)

M   x I , donde I es el intervalo donde se busca la raíz, entonces
2 min f ' ( x)

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

Teorema: Sea f ( x)  0 una función, con raíz única r en a, b . Si

f C 2 a, b f ' ( x)  0  x a, b, entonces si x0  x1  V 1 (r ) , la sucesión de puntos

M

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

13 Leonardo Flores González

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

Sea   maxME0 , ME1 entonces si se cumple MEk  2    MEk 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 ki  1

5


   
1 5
2
i 1
1 5
2
i 1
 , como 0  lim ME  lim  k i  0 , se tiene
 i 
i
i 

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 En 1  CE n En 1 , donde

f ' ' ( )
C ya que cuando n     r   r , podemos determinar el orden de
2 f ' ( )

convergencia del método con la siguiente suposición En 1  KE np , de donde En  KE np1 , que

1 1
1 1
reemplazando en En 1  CE n En 1 resulta K E  CEn p
n
p p
, de donde K  C p
 p (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.

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

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

x1  x0
Se tiene x 0  1 y x1  1.05 entonces x1  x0  0.05  tol , x2  x1  f1  1.097064 . Luego
f1  f 0
se procede en forma similar para las demás iteraciones.

iteración xn 1 f ( xn  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
n
Sea P( x)  a x k
k
un polinomio de grado n , si ak  R , k  0,1,..., n , el método de Bairstow es
k 0

ú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 .

n
Teorema. Si el polinomio P( x)  a x
k 0
k
k
se divide entre el factor cuadrático x 2  ux  v , entonces

n
el cociente Q( x)  b x
k 2
k
k 2
y el residuo R ( x)  b1 ( x  u )  b0 , pueden calcularse de manera

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

Se tiene que:

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

 a k x k   bk x k  2 ( x 2  ux  v)  b1 ( x  u )  b0
k 0 k 2

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

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

bk b
R (x) sea cero bk (u , v)  0 para k  0,1, de esta forma definimos ck   d k  k 1 donde
u v
k  0,1,2,..., n , si derivamos bk  ak  ubk 1  vbk  2 con respecto a u se tiene que

c k  bk 1  uc k 1  vc k  2 , definimos de esta manera la relación de recursividad

ck  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  d n  0 , de estas relaciones de recursividad se puede ver

bk b
que ck  d k . Linealizando bk (u  u , v  v)  0 se tiene bk (u, v)  u  k v  0 para
u v
k  0,1 resultando:

c0 c1 u b c c1
  0 donde J  0 .
c1 c2 v b1 c1 c2

16 Leonardo Flores González

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

Teorema.: Sean u0 ,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  ak  ubk 1  vbk  2 con bn 1  bn  2  0 .

3. Para k  n, n  1,...0 hacer ck  bk  uck 1  vck  2 con cn 1  cn  2  0 .

c0 c1 u b
4. Resolver  0 .
c1 c2 v b1

u u u
5. Hacer   .
v v v

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

CODIFICACION EN MATLAB

function [u,v]=bairstow(a,u,v,it,tol)
n=length(a);delta=[tol+10 tol+10]';
for j=1:it
if((abs(delta(1))+abs(delta(2)))>tol)
[b]=formula(a,u,v,n);
[c]=formula(b,u,v,n-1);
delta=[c(n-2) c(n-3);c(n-1) c(n-2)]\[-b(n-1) -b(n)]';
u=u+delta(1);v=v+delta(2);
end
end

function [m]=formula(f,u,v,nn)
m(1)=f(1);m(2)=f(2)+u*f(1);
for k=3:nn
m(k)=f(k)+u*m(k-1)+v*m(k-2);
end

17 Leonardo Flores González

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

Capítulo 2

SISTEMA DE ECUACIONES LINEALES

18 Leonardo Flores González

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

METODOS DIRECTOS DE SOLUCION DE SISTEMAS DE ECUACIONES LINEALES

Los sistemas de ecuaciones lineales son importantes en la solución de problemas relacionados al

campo de la ingeniería, frecuentemente se presentan al resolver armaduras, sistemas
aporticados, etc. Estos problemas involucran sistemas de ecuaciones lineales con un gran número
de variables, dependiendo del orden del sistema que se quiera solucionar, existen algoritmos
recursivos como los métodos de Jacobi, Gauss-Seidel y métodos directos como la eliminación
gaussiana.

Por ejemplo en una armadura espacial, se resuelve un sistema de ecuaciones lineales para
encontrar los desplazamientos verticales y horizontales en cada uno de sus elementos.

Antes de continuar se hace mención a algunas propiedades de matrices y sistemas lineales, útiles
para la comprensión de la solución de los sistemas lineales por métodos directos.

Propiedades:
1. La inversa de una matriz triangular inferior o superior es otra matriz triangular inferior o superior.
2. Si un sistema de ecuaciones se obtiene a partir de otro por medio de una sucesión finita de
operaciones elementales, los dos sistemas son equivalentes.
3. Una matriz cuadrada puede tener a lo sumo una matriz inversa por la derecha.

19 Leonardo Flores González

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

Definición: Se llamará matriz elemental triangular inferior a una matriz de la forma:

Ek  I  mk ekT donde mk  0,0,,0, mk 1,1 ,, mn ,1  , sea ek  0,0,, ek ,1 ,,0 donde e k ,1  1, I

T T

es una matriz identidad de orden n . Esta matriz presenta las siguientes propiedades:

1. Ek1  I  mk ekT .
k
2. Ek Ek 1  E2 E1  I   m j eTj .
j 1

Sea Ax  b un sistema de n ecuaciones lineales, donde A es una matriz cuadrada de orden n ,

b y x vectores de orden nx1 ; para resolver el sistema de ecuaciones por eliminación gaussiana
procederemos a descomponer la matriz A  LU . A continuación se indican condiciones
suficientes para que la factorización A  LU sea posible.

a. Si los n menores principales de la matriz A son no singulares, entonces la matriz A tiene

descomposición LU .

b. Si A  A , es matriz real simétrica y definida positiva, entonces tiene factorización única

T

A  LLT , donde L , es una matriz triangular inferior con los elementos de la diagonal positivos.
c. Se consideran mij , k   a ij , k / aki , k , j  k  1 , para factorizar A  LU .

Teorema: Si A es una matriz que cumple con la propiedad (a), esta es una condición suficiente
para que exista una matriz L triangular inferior y una matriz U triangular superior, ambas de
orden n tal que A  LU .

 U k 1, k 1 M k 1, n  k 1 
Si A( k 1)  Ek 1Ek  2  E1 A    , luego A( k )  Ek A( k 1) , 1  k  n  1 , sea
 0n  k 1, k 1 An  k 1, n  k 1 
*

 I k 1, k 1 0n  k 1, n  k 1 
Ek   *
 , donde Ek* es una matriz elemental triangular inferior de orden

 0n  k 1, n  k 1 Ek 
n  k  1 , que hace cero los elementos debajo de la diagonal principal de la primera columna de
An* k 1,n  k 1 , luego se tiene:

20 Leonardo Flores González

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

 I k 1, k 1 0 n  k 1, n  k 1  U k 1, k 1 M k 1, n  k 1 

A( k )   
 0

 0n  k 1, n  k 1 Ek*  n  k 1, k 1 An* k 1, n  k 1 
 U k ,k M nk ,nk 
A( k )   *


 0n  k , k An  k , n  k 

luego A( n 1)  En 1En  2  E1 A donde A( n 1)  U , de la propiedad (2) se tiene que:

L  E11  En11 es una matriz triangular inferior con unos en la diagonal. Se puede probar que esta
descomposición con unos en la diagonal de L es única. Este algoritmo toma el nombre de
factorización LU de Doolittle.

Para resolver un sistema de ecuaciones se procede como se indica a continuación:

1. Se factoriza A  LU de modo que LUx  b , se hace Ux  z , de modo que Lz  b.

i 1
2. Resolver Lz  b por sustitución directa: zi  (bi  l
k 1
z ) / li ,i , 2  i  n, z1  b1 / l1,1.
i,k k

n
3. Resolver Ux  z por sustitución inversa: xi  (bi  u
k  i 1
z ) / ui ,i , n  i  2, xn  bn / an, n .
i,k k

Ejemplo: Se tiene una Armadura como la que se muestra en la figura (a) la cual soporta fuerzas
de 20000 y 25000 libras como se indica en la figura (b). El area de las secciones transversales de
la armadura es de 1 in 2 y E  29 .5  10 6 psi , se desea encontrar los desplazamientos verticales y
horizontales en los puntos en los puntos 2 y 3.

(a)

21 Leonardo Flores González

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

(b)

El sistema de ecuaciones lineales que resulta luego de formular la matriz de rigidez de la

estructura es:

15 0 0 U H 2   20000 

29.5  106     
 0 22.68 5.76  U H 3    0 
600     
 0 5.76 24.35  UV 3    25000 

Donde U H 2 ,U H 3 y U V 3 representan los desplazamientos horizontal en el punto 2, horizontal en el

punto 3 y vertical en el punto 3 respectivamente, descomponiendo la matriz cuadrada en su forma

LU se tiene:

1 0 0  0.7375 0 0 
   
L  0 1 0 , U  106  0 1.1151 0.2832 
 0 0.2540 1   0 1.1238 
   0

Resolviendo Lz  b por sustitución directa se tiene z  20000 ,0,25000  .

T

Resolviendo Ux  z por sustitución inversa se tiene x  10 3 (27 .12,5.65,22 .25)T todos los
desplazamientos están en pulgadas.

22 Leonardo Flores González

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

ALGORITMOS PARA LA FACTORIZACION LU

La factorización LU consiste en representar una matriz A como el producto de una matriz

triangular inferior L y otra matriz triangular superior U es decir A=LU.

l11 0 0  0 u11 u12 u13  u1n 

l  0  0 u u23  u2 n 
 21 l22 0  22

L  l31 l32 l33  0 

U  0 0 u33  u3n 
   
          
ln1 ln 2 ln 3  lnn   0 0 0  unn 

min( i , j )
Además el elemento aij  l u
s 1
is sj ya que L y U son matrices triangulares.

k k 1 k 1
Si i  j se tiene: akk   lksusk  lksusk  lkkukk de donde ukk  (akk   lksusk ) / lkk .
s 1 s 1 s 1

k 1 k 1
Si k  1  j  n se tiene akj   lksusj  lkkukj de donde ukj  (akj   lksusj ) / lkk .
s 1 s 1

k 1 k 1
Si k  1  i  n se tiene aik   lisusk  lik ukk de donde lik  (aik   lisusk ) / ukk .
s 1 s 1

Cuando lkk  1 se llama el algoritmo de Doolittle que es el que desarrollaremos:

Algoritmo 1

Paso 1. Multiplicar la primera fila de L por todas las columnas de U para obtener la primera fila de
U.

Paso 2. Multiplicar todas las filas de L por la primera columna de U para obtener la primera
columna de L.

Paso 3. Repetir el Paso 1 es decir multiplicar la segunda fila de L por las Columnas de U a partir
de la segunda columna para obtener los elementos de la segunda fila de U. (Si se desea
obtener la k-ésima fila de U se multiplicara la k-ésima fila de L por las columnas de U a
partir de la k-ésima columna).

Paso 4. Repetir el Paso 2 es decir multiplicar las filas de L a partir de la segunda fila por la
segunda columnas de U para obtener los elementos de la segunda de L. (Si se desea
obtener la k-ésima fila de L se multiplicaran las filas de L a partir de la k-ésima fila por la
k-ésima columna de U).

23 Leonardo Flores González

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

function [L,U]=lablu(a) %Algoritmo de Doolittle

n=length(a);
for k=1:n
L(k,k)=1;
for j=k:n
suma=0;
for s=1:k-1
suma=suma+L(k,s)*U(s,j);
end
U(k,j)=a(k,j)-suma;
end
for i=k:n
suma=0;
for s=1:k-1
suma=suma+L(i,s)*U(s,k);
end
L(i,k)=(a(i,k)-suma)/U(k,k);
end
end

Algoritmo 2
Como segundo algoritmo desarrollaremos en matlab el criterio expuesto en la parte teórica del
presente capítulo.

function [L,U]=lu_2(A) %Algoritmo de Doolittle

n=length(A);
for i=1:n-1
A(i+1:n,i)=A(i+1:n,i)/A(i,i);
A(i+1:n,i+1:n)=A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n);
end
L=tril(A,-1)+eye(n);
U=triu(A);

24 Leonardo Flores González

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

Algoritmo 3
El método PA  LU , es una variante del anterior algoritmo, se permutan las filas en caso de
encontrar un cero en la diagonal principal. Pij es la matriz que resulta de intercambiar las fila i y

j en la matriz I nxn . A continuación indicamos el algoritmo en matlab:

function [L,U,P]=palu_2(A)

n=length(A); P=eye(n);

for i=1:n-1

[aux,pos] = max(A(i:n,i:n)); j=pos(1)+i-1;

c=A(i,:);A(i,:)=A(j,:);A(j,:)=c;

c=P(i,:);P(i,:)=P(j,:);P(j,:)=c;

A(i+1:n,i)=A(i+1:n,i)/A(i,i);

A(i+1:n,i+1:n)=A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n);

end

L=tril(A,-1)+eye(n);

U=triu(A);

Algoritmo 4
El algoritmo de cholesky será una variante del LU cuando A  A , a cada fila de la matriz U se
t

le dividirá por uii ; será suficiente que los valores de uii sean mayores que cero, para que todos

lo elementos de la matriz L sean números reales.

Si consideramos que la descomposición de Doolittle es única y además que A  A , entonces se

t

cumple que A  LU  U t Lt , lo que indica que:

 1  0   u1,1  u1, n   u1,1  0   1  ln ,1 

     
A     0         0    0  ,
l     
 n ,1  1   0  un , n   u1, n  un, n   0  1 

si multiplicamos a la columna i-ésima de la matriz triangular inferior por 1 / ui , i y a la i-ésima fila de

la matriz triangular superior por ui ,i , el producto no se altera y se tiene:

25 Leonardo Flores González

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

 
 1  0   u1,1  ln ,1u1,1 
  
A   0   0 
 u1, n  
u  1   0  un , n 
 1,1 

por la unicidad de la factorización de Doolittle se tiene que:

   
 1  0  1  0  1  0   u1,1  u1, n 
      
   0      0  entonces A     0       , si multiplicamos a la
 u1, n     u1, n  
u  1   ln ,1  1  u  1   0  un , n 
 1,1   1,1 

columna i-ésima de la matriz triangular inferior por ui ,i y a la i-ésima fila de la matriz triangular

superior por 1 / ui , i , se tiene:

   u 
u1, n 

  0   1,1

u1,1
 u1,1 
A   0      , lo que indica que A  LLt
 u1, n  
  un , n   0  un , n 
u1,1  
  

function [L, LT]=cholesky(A)

n=length(A);

for i=1:n-1

A(i+1:n,i)=A(i+1:n,i)/A(i,i);

A(i+1:n,i+1:n)=A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n);

A(i,i:n)=A(i,i:n)/A(i,i)^0.5;

end

A(n,n)=A(n,n)^0.5;LT=triu(A);L=LT';

26 Leonardo Flores González

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

A continuación se explica un ejemplo ilustrativo:

Dada una matriz de orden 40x40 como se indica a continuación:

A 0 A 0  0
 
0 B 0 B  B
0 0 A 0  0  1 2
  x  c donde c  0 1 0 0 1 0 1 0 0 1 0 t A  
0 0 B  B ,  2 1  ,
0  
     

0 0 0 0 B 
 0

 1 2 3
 
B   2 1 0  . Resuelva el sistema usando factorización LU . Dar como respuesta x t x.
 2 1 1
 

Se puede hacer la siguiente descomposición:

A 0 A 0  0   LA 0 0 0  0  U A 0 U A 0  0 
    
0 B 0 B  B  0 LB 0 0  0   0 UB 0 UB  UB 
0 0 A 0  0  0 0 LA 0  0  0 0 UA 0  0 
   
0 0 0 B  B  0 0 0 LB  0  0 0 0 UB  UB 
      
                  
0 0 0 0 0 B   0   U B 
 0 0 0 0 LB   0 0 0 0 0

Para resolver el sistema se tiene:

U A 0 U A 0  0 
 
 0 UB 0 UB  UB 
 0 0 UA 0  0 
  x  z , entonces el sistema para la sustitución directa es:
 0 0 0 UB  UB 
       

 0 U B 
 0 0 0 0

27 Leonardo Flores González

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

 LA 0 0 0  0
 
0 LB 0 0  0
0 0 LA 0  0
  z  c , se puede ver que c  c1 c2 c3 c4  , donde
t

0 0 0 LB  0
       

0 LB 
 0 0 0 0

c2  c4  c6  c8  0 1 0  , c1  c3  c5  c7  0 1 . Descomponiendo a las matrices A y B se

tiene:

1 0 01 2 3
 1 01 2   
A      y B   0 1 00 1 0  podemos realizar los siguientes cálculos por
 2 1   0  3  2  3 1   0 0  5
  
sustitución directa:

LA z2i 1  c2i 1 , LB z2i  c2i1 donde i  1,2,3,4 de esta manera se obtiene:

z1  z3  z5  z7  0 1 t , z2  z4  z6  z8  0 1 3 t .

Por sustitución inversa obtenemos:

1 3
U B x8  z8 , luego x8   1   t , luego se tiene que U B x6  U B x8  z6  z8 luego
5 5

x6  0 0 0 , de igual forma U B x4  U B x6  U B x8  z4  z8 , luego x4  0 0 0  , finalmente

t t

U B x2  U B x4  U B x6  U B x8  z2  z8 , luego x2  0 0 0  .
t

De la misma forma para los valores de luego x impares se tiene:

2 1
U A x7  z7 , luego x7     t , luego se tiene que U A x5  U A x7  z5  z7 luego
3 3

x5  0 0 0 , de igual forma U A x3  U A x5  U A x7  z3  z7 , luego x3  0 0 0 , finalmente

t t

U A x1  U A x3  U A x5  U A x7  z1  z7 , luego x1  0 0 0  .
t

88
La respuesta sería: xt x  x7 2  x8
2 2
2
 .
45

28 Leonardo Flores González

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

METODO DE JACOBI

Sea Ax  b un sistema de n ecuaciones lineales, si aii  0 donde i  1,..., n , el sistema se puede

bi n a aij
escribir de la siguiente forma xi    x j .sea M una matriz tal que mij  
ij
si
aii j 1 aii aii
j i

bi
j  i  mii  0 , sea C un vector tal que ci  . El sistema se puede escribir de la forma
aii
x  Mx  C .

Definición: El vector error en la k-ésima iteración se denota por E ( k )  x  x ( k ) tal que

 xi  xi , i  1,..., n .
(k ) (k )
ei

Definición: Se llamará error absoluto en la k-ésima iteración a  ( k )  E ( k ) .



Definición: Se llamará factor de convergencia de M al número   M 

.

Definición: Sea el vector diferencia D ( k )  x ( k )  x ( k 1) se llamará k-ésima diferencia absoluta a

 (k )  D(k ) .


Teorema 1 : Una condición suficiente para que el método de Jacobi converja a la solución de
x  Mx  C para cualquier aproximación inicial x ( 0) es que   1 .

Sea x ( k )  Mx ( k 1)  C  E ( k )  ME ( k 1) , k  1,..., n , tomando valores absolutos se tiene:

n
  mij ei donde i  1,..., n  ei ( k )    ( k 1) luego  ( k )    ( k 1)   ( k )   k  ( 0) de donde
(k ) ( k 1)
ei
j 1

lim  k  0 por tanto lim x ( k )  x .

k  k 

Teorema 2 : Sea el sistema Ax  b . Una condición suficiente para que el método de Jacobi sea
convergente a la solución del sistema es que A sea diagonalmente dominante.
Es una consecuencia del teorema 1.

29 Leonardo Flores González

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

 (k )
Teorema 3 : Si el factor de convergencia de M es   1   ( k )   .
1
E ( k 1)  E ( k )  D ( k )  E ( k 1)  E (k )  D(k ) ,   ( k 1)   ( k )   ( k ) , además  ( k )    ( k 1) de donde
  

 (k )
se obtiene  ( k )   .
1
 (k )
Se puede escoger como criterio de parada   Tol . Si consideramos que el sistema es de
1
n ecuaciones con n incógnitas, que tiene solución única y que A es una matriz diagonalmente
dominante entonces podemos considerar el algoritmo siguiente:

Algoritmo

1. Para k  0 escoger una aproximación inicial x ( 0) , asignar a  ( 0) un valor grande, indicar tol .

2. Hacer k  k  1 , x ( k )  Mx ( k 1)  C ,  ( k )  D ( k ) .


 (k )
3. Si   tol repetir para k  1,2,...
1
4. x (k ) es una aproximación con un error menor que tol .

Programa en matlab
function [x1,it]=jacobilab(A,b,xo,tol)
x1=xo;xo=xo+10*tol;it=0;n=length(b);
while (norm(x1-xo,inf)>tol)
xo=x1;
for i=1:n
x1(i)=(b(i)-A(i,:)*xo+A(i,i)*xo(i))/A(i,i);
end
it=it+1;
end

Se puede notar que si Ax  b  Qx  (Q  A) x  b de donde x  ( I  Q 1 A) x  Q 1b , esta última

fórmula sugiere x ( k )  ( I  Q 1 A) x ( k 1)  Q 1b , en el método de Jacobi Q es una matriz diagonal tal

que qii  aii , i  1,..., n .

El algoritmo de Gauss-Seidel es similar al algoritmo de Jacobi con la particularidad de que la
matriz Q es la matriz triangular inferior de A incluyendo su diagonal. Las condiciones de
suficiencia son las mismas es decir A debe ser diagonalmente dominante.

30 Leonardo Flores González

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

Programa en matlab

function [xo,it]=gaussseidellab(A,b,xo,tol)
u=xo+10*tol;it=0;n=length(b);
while (norm(u-xo,inf)>tol)
u=xo;
for i=1:n
xo(i)=(b(i)-A(i,:)*xo+A(i,i)*xo(i))/A(i,i);
end
it=it+1;
end

Un ejemplo donde se aplican sistemas de ecuaciones y luego estos son resueltos con ayuda de
algoritmos LU o similares, es el caso de sistemas aporticados, por ejemplo dado el pórtico que
indicamos en la figura, deseamos conocer los desplazamientos de cada uno de los elementos.

La Longitud de las vigas es de 800 cm., la altura de las columnas es de 450 cm., las cargas
verticales son de 50 toneladas y la carga horizontal es de 10 toneladas, con ayuda de estos datos
más las respectivas inercias, áreas y módulos de elasticidad de cada elemento, se construye la
matriz de rigidez del pórtico. Par esto, en los anexos se entrega un programa para construir la
matriz de rigidez de una estructura bidimensional, este programa está basado íntegramente en los
apuntes de clase del Dr. Hugo Scaletti Farina (1996).

31 Leonardo Flores González

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

El sistema que se tiene que resolver es de 18 incógnitas, que se reduce a 12 incógnitas,

considerando las condiciones de frontera del mismo.

El sistema tiene la forma: Ku  f , donde K es la matriz de rigidez de la estructura y f es el vector

de fuerzas que en nuestro caso sería:

f  0 0 0 0  50 0 10  50 0 0 0 0 0  50 0 0  50 0
T

La matriz de rigidez sería:

6584.36214 0 -1481481.48 -6584.36214 0 -1481481.48 0 0 0 0 0 0 0 0 0 0 0 0

0 666666.667 0 0 -666666.667 0 0 0 0 0 0 0 0 0 0 0 0 0
-1481481.48 0 444444444 1481481.48 0 222222222 0 0 0 0 0 0 0 0 0 0 0 0
-6584.36214 0 1481481.48 263168.724 0 0 -6584.36214 0 -1481481.48 0 0 0 -250000 0 0 0 0 0
0 -666666.667 0 0 1333919.27 234375 0 -666666.667 0 0 0 0 0 -585.9375 234375 0 0 0
-1481481.48 0 222222222 0 234375 1013888889 1481481.48 0 222222222 0 0 0 0 -234375 62500000 0 0 0
0 0 0 -6584.36214 0 1481481.48 256584.362 0 1481481.48 0 0 0 0 0 0 -250000 0 0
0 0 0 0 -666666.667 0 0 667252.604 234375 0 0 0 0 0 0 0 -585.9375 234375
K= 0 0 0 -1481481.48 0 222222222 1481481.48 234375 569444444 0 0 0 0 0 0 0 -234375 62500000
0 0 0 0 0 0 0 0 0 6584.36214 0 -1481481.48 -6584.36214 0 -1481481.48 0 0 0
0 0 0 0 0 0 0 0 0 0 666666.667 0 0 -666666.667 0 0 0 0
0 0 0 0 0 0 0 0 0 -1481481.48 0 444444444 1481481.48 0 222222222 0 0 0
0 0 0 -250000 0 0 0 0 0 -6584.36214 0 1481481.48 263168.724 0 0 -6584.36214 0 -1481481.48
0 0 0 0 -585.9375 -234375 0 0 0 0 -666666.667 0 0 1333919.27 -234375 0 -666666.667 0
0 0 0 0 234375 62500000 0 0 0 -1481481.48 0 222222222 0 -234375 1013888889 1481481.48 0 222222222
0 0 0 0 0 0 -250000 0 0 0 0 0 -6584.36214 0 1481481.48 256584.362 0 1481481.48
0 0 0 0 0 0 0 -585.9375 -234375 0 0 0 0 -666666.667 0 0 667252.604 -234375
0 0 0 0 0 0 0 234375 62500000 0 0 0 -1481481.48 0 222222222 1481481.48 -234375 569444444

Si empleamos el algoritmo LU Propuesto obtenemos como solución:

u= 0 0 0 2.3936 -0.1404 -0.0073 6.2579 -0.2108 -0.0065 0 0 0 2.3935 -0.1596 -0.0073 6.238 -0.2392 -0.0065 .

32 Leonardo Flores González

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

Capítulo 3

VALORES Y VECTORES PROPIOS

33 Leonardo Flores González

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

VALORES Y VECTORES PROPIOS

Se estudia el problema clásico A   , antes de enunciar algún algoritmo para la solución del
problema clásico, se indican algunos teoremas importantes.

Definición: Una matriz cuadrada de orden n , que pertenece al conjunto de matrices C nxn es

unitaria si A A  I donde A*  ( A)T

*

Definición: Se llamará a una matriz hermitiana si A  A .

*

Definición: Una matriz es definida positiva si x* Ax  0  x C n  0.

Definición: Se dice que A  B son matrices cuadradas semejantes, si P / B  P 1 AP .

Se debe notar que si A  R nxn  aij  R  A*  AT .

Teorema 1 Todas las matrices semejantes tiene los mismos valores propios.

Teorema 2 (Schur) Toda matriz cuadrada es unitariamente semejante a una matriz triangular.

Sea x1  C n / Ax1  1 x1 un vector tal que x 1* x1  1, se puede encontrar un vector

h 1
xh  C n / xh   h (eh   mkh xk ) tal que xi* xh   ih , donde xh   hVh , mih   xi*eh , xi  C n ,
k 1

i  N / 1  i  h y  h  1 / Vh 2 .

Empleando la construcción anterior se construye la matriz T1  x1 x2 ...xn 1 , la cual es una matriz

unitaria, supongamos que toda matriz de orden n es unitariamente semejante a una matriz
triangular, entonces se puede ver que:

 v1xn   1 01xn 
T1* AT1   1  si T2* BT2 es una matriz triangular superior y T3    , donde T2  T3
 0nx1 Bnxn   0nx1 T2 
son matrices unitarias, se puede ver también que T  T3T1 es una matriz unitaria, comprobándose

finalmente que:

34 Leonardo Flores González

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

 c1xn 
T * AT   1 *
 es una matriz triangular superior.
 0nx1 T2 BT2 

Teorema 3 Una matriz Anxn es semejante a una matriz diagonal Dnxn si y solamente si Anxn tiene

n vectores propios linealmente independientes.

1
Si X AX  D entonces det( X )  0 , luego los vectores columna de X son L.I., sea xk la k-

ésima columna de X se puede ver que Ax k  d kk xk , luego xk es un vector propio de A .

Si xk son vectores L.I. de A donde k  1,2,...n , luego AX  Ax1    xn   1 x1    n xn  de

donde:

 1  0 
 
AX  x1    xn      XD , como xk son L.I., se tiene X 1 AX  D .
0   
 n

Consecuencias de los teoremas 2 y 3

1. Toda matriz Hermitiana tiene valores propios reales.

2. Si A  A , existe una matriz unitaria X / X * AX  D , D es matriz diagonal.
*

3. Si A  AT  aij  R , existe una matriz ortogonal X / X T AX  D , D es matriz diagonal.

4. Si A es una matriz real y simétrica tiene n vectores linealmente independientes.

5. Una matriz Hermitiana es definida positiva si y solamente si sus valores propios son positivos.

Propiedades de las matrices ortogonales

T
1. Si A es una matriz real y simétrica, y P una matriz ortogonal entonces P AP es simétrica.
n
2. Si Pj representa una matriz ortogonal, entonces  P es una matriz ortogonal.
k 1
k

35 Leonardo Flores González

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

Método de Jacobi Clásico

Este método es utilizado para la solución de A   , donde A es una matriz cuadrada real y
simétrica. La idea de jacobi es la siguiente:
n n

lim  J A J k  D , donde J k es una matriz ortogonal y D una matriz diagonal.

T
nk
n   k 1 k 1

1
 n n 2
   aij
2
La norma de Frobenius para A  R nxn se define como A F  . Dada la matriz:

 j 1 i 1 

1  0  0  0
 
      
0  c  s  0  p
 
J ( p , q, )        
0   s  c  0   q

      
 
0  0  0  1
 
p q

Una matriz ortogonal donde c  cos   s  sen . El proceso consiste en:

A( k 1)  J ((pk ,)q , ) A( k ) J (( pk ,)q , )  k  0, donde A( 0)  A,

T
1. p, q indican la posición del

máximo aij / i  j .

2. Luego de (1) se tiene aij( k 1)  aij( k ) si i, j son diferentes de p, q , en caso contrario se tiene
k 1)
a (pm ( k 1)
 amp  camp
(k )
 samq
(k )
donde m  1,..n  p, q , del mismo modo se tiene que
( k 1)
aqm ( k 1)
 amq  samp
(k )
 camq
(k )
donde m  1,..n  p, q.

c 2  s 2 aqq  a pp
(k ) (k )
s
3. Para que a ( k 1)
pq a ( k 1)
qp  0 , se tiene se tiene que    (k )
, si t   tan
2cs 2a pq c

signo( ) 1
luego t 2  2 t  1  0 , de donde t  ,c y s  ct.
  1 2 1 t2

4. a (ppk 1)  a (ppk ) c 2  aqq

(k ) 2 ( k 1)
s  2a (pqk ) cs y aqq  a (ppk ) s 2  aqq c  2a (pqk )cs.
(k ) 2

( k 1) ( k 1)
De lo anterior se tiene que amp  amq  amp  amq
2 2 2 2
(k ) (k )
5. , sea W (k ) la suma de los cuadrados
( k 1)
de todos los elementos fuera de la diagonal de A
(k )
, luego W  W ( k )  2a 2pq , sea N la

36 Leonardo Flores González

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

mitad del número de elementos fuera de la diagonal, entonces W ( k )  2 Na 2pq , se tiene

k
( k 1)  1  1
finalmente que W  1  W ( k ) de donde W ( k )  1   W ( 0) , si k    W ( k )  0.
 N  N
Algoritmo
1. Para k  0 ingresar A y tol .
(0)

2. V  I n , eps  tol A( 0 ) .
F

3. Mientras W ( k )  eps

3a. Escoger p, q / a (pqk )  max i  j aij( k ) .

3b. Calcular  , J y t .

3c. A( k 1)  J (( pk ,)q , ) A( k ) J (( pk ,)q , )  k  0.

T

3d. V  VJ (( pk ,)q , ) .

Codificación en Matlab
function [d,O]=jacobiclasico(A,tol)
d=A; V=eye(length(A));
while norm(d-diag(diag(d)),'fro')>tol*norm(A,'fro')
[maximos,filas]=max(abs(triu(d,1)));
[maximo,q]=max(maximos); p=s(q);
ang=atan(2*d(p,q)/(d(q,q)-d(p,p)))/2;
J=eye(length(d));
J([p,q],[p,q])=[cos(ang) sin(ang);-sin(ang) cos(ang)];
d=J'*d*J;V=V*J;
end

37 Leonardo Flores González

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

METODO DE LA POTENCIA

Es un método empleado para encontrar el valor propio dominante con su vector propio asociado.
Para poder emplear este método es necesario que se cumpla los siguiente:

1. Exista un valor propio dominante es decir: 1  2  3   n .

2. La matriz tiene n vectores propios linealmente independientes.

El problema clásico es el siguiente A  . Sea x0 un vector tal que x 0 2

 1 , sea u1, u2 ,un

los vectores propios de A que generan R n , de lo anterior podemos decir Au i  i u i , donde

n
i  1,2...n , además x0   uk . Se debe considerar que u i  (u i1 , u i 2 ,....u in ) T .
k 1
Para deducir la convergencia del método al valor propio dominante se tienen que seguir los
siguientes pasos:
n n n

  k uk , de esta manera x1 
y1 1
y1  Ax 0  y1  Au k 
y1 2

y1
 u k k , la primera
k 1 k 1 2 k 1

aproximación para el valor propio será 1  x1T Ax1 ; si se repite lo anterior tenemos
n n n
1 1 y2 1
y2  Ax1  y2 
y1
 k Auk  y1
 2k uk , de esta manera x2  
y 2 2 s2
 u 2
k k , donde
2 k 1 2 k 1 k 1

n
s2   u
k 1
2
k k , luego 2  x2T Ax 2 .

n n
1 ym 1
De esta manera podemos decir: ym 
sm1
 mkuk  xm 
k 1 ym 2
, luego xm 
sm
 u
k 1
m
k k donde

 j
m
n n

sm   v además vm   (u1i   
m
 u ji ) 2 .
j  2  1
1 m
i 1 

1m n
 m  1 n
 m 
De estas últimas expresiones se puede decir xm  (u1    mk uk )  (u1    mk uk ) ,
sm k 2  1  vm k 2  1 

luego m  xmT Ax m , de esta forma tenemos:

1 T n
 mk  T 1 n
 mk  1 T n  mk  T n
 mk1 
m  (u1    m uk )( ( Au1    m  Au k )  2 (u1    m uk )(u1    m1 uk )
vm k 2  1  vm k 2  1  vm k 2  1  k 2  1 
 im  T n   j 
m1
1  mk  T n n
Donde m  2 (u u    m u1 uk    m ui  m1 u j ) , si m   resulta:
T 
vm k 2  1 
1 1
i 2  1 

j 1  1


38 Leonardo Flores González
 UNIVERSIDAD NACIONAL DE INGENIERÍA
FACULTAD DE INGENIERÍA CIVIL

1
lim m  ( u1 2 )  1 .
2

m  2
u1 2

De lo anterior podemos deducir el siguiente algoritmo:

Algoritmo
1. Ingresar A, x0 , tol.
2. xk  Ax k 1  k  1.

xk
3. xk   k  1.
xk 2

4. k  xkT Ax k .

5. Mientras k  k 1  tol repetir los pasos 2,3,4,5.

6. La respuesta es  k y xk .

El algoritmo propuesto en matlab es:

function [x,L,it]=potencia(A,x0,tol)
x0=x0/norm(x0,2);it=1;
L0=x0'*A*x0;x1=A*x0;x1=x1/norm(x1,2);L1=x1'*A*x1;
y=[it x1' L1 abs(L1-L0)];n=length(x1);
[s1]=impresion(y,n);
while abs(L1-L0)>tol
x0=x1;L0=L1;it=it+1;
x1=A*x0;
x1=x1/norm(x1,2);
L1=x1'*A*x1;
y=[it x1' L1 abs(L1-L0)];
[s1]=impresion(y,n);
end
x=x1;L=L1;

39 Leonardo Flores González

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

Esta rutina es empleada para la impresión de los resultados:

function [s1]=impresion(y,n)
CC = fopen('[Link]','a');
s1='%d\t';
for i=1:n
s1=strcat(s1,'%8.5f\t');
end
s1=strcat(s1,'%8.5f\t');
s1=strcat(s1,'%8.7f\n');
fprintf(CC,s1,y);
fclose(CC);

Se debe tener en cuenta que el método de la potencia inversa es similar al método de la directa
solo que la matriz es A 1 .

De esta manera el algoritmo de la potencia inversa resulta:

function [x,L,it]=potencia_inversa(A,x0,tol)
x0=x0/norm(x0,2);it=1;A=inv(A);
L0=x0'*A*x0;x1=A*x0;x1=x1/norm(x1,2);L1=x1'*A*x1;
y=[it x1' L1 abs(L1-L0)];n=length(x1);
[s1]=impresion(y,n);
while abs(L1-L0)>tol
x0=x1;L0=L1;it=it+1;
x1=A*x0;
x1=x1/norm(x1,2);
L1=x1'*A*x1;
y=[it x1' L1 abs(L1-L0)];
[s1]=impresion(y,n);
end
x=x1;L=1/L1;

40 Leonardo Flores González

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

Como ejemplo Ilustrativo encontraremos el valor propio dominante de:

 2 1 0
 
A   1 2 1  con x0   1
3
, 1
3
, 
1 T
3
y tol  0.00001 .
 0 1 2
 

 0.5145 
x1  
0  x Ax 0  3.33333 , x1  Ax 0 y x1 
T
0   0.6859  , 1  x1T Ax1  3.41176
x1 2  
 0.5145 
1  0  0.0784  tol. Se repiten los mismos pasos hasta que k  k 1  tol. De acuerdo ha
esto tenemos:

Iteraciones xi,1 xi,2 xi,3 Lambda Error

x1 0.51450 0.68599 0.51450 3.411760 0.078431
x2 0.50252 0.70353 0.50252 3.414140 0.002377
x3 0.50043 0.70649 0.50043 3.414210 0.000070
x4 0.50007 0.70700 0.50007 3.414210 0.000002

La respuesta es  4  3.414210 y x4  0.50007 ,0.70700 ,0.50007  .

T

41 Leonardo Flores González

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

Capítulo 4

INTERPOLACION Y APROXIMACION POLINOMIAL

42 Leonardo Flores González

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

INTERPOLACIÓN DE FUNCIONES CON POLINOMIOS

Los polinomios de interpolación se utilizan para representar funciones cuando no es posible

conocer la función en forma explícita o se conoce un número finito de puntos de una función, sea
f :    , la función se aproxima: f ( x)  p ( x)  E ( x) , donde p es el polinomio de
interpolación y E el error cometido al evaluar f en un punto de su dominio.

Figura 1

En general se dice que el polinomio p de grado n interpola f en los puntos x i si

p ( xi )  f ( xi )  0  i  n .

Teorema: Si x0 , x1 ,..., x n son números reales distintos, entonces para valores arbitrarios

y 0 , y1 ,..., y n existe un polinomio único p n de a lo sumo grado n tal que:

p n ( xi )  y i  0  i  n .

43 Leonardo Flores González

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

FORMAS DE REPRESENTAR EL POLINOMIO DE INTERPOLACIÓN

El polinomio de interpolación puede ser representado por la forma de Newton o por la forma de
Lagrange.

n k 1
Polinomio de Interpolación de Newton: p n ( xi )  c0   ck  ( x  x j ), c k representan las
k 1 j 0

diferencias divididas como veremos a más adelante.

n n (x  x j )
Polinomio de Interpolación de Lagrange: p n ( x )  y l
k 0
k k ( x) , l k ( x)  
j 0 ( xk  x j )
jk

ERROR DE INTERPOLACIÓN POLINOMIAL

Teorema: Si f  C n 1 a, b y p un polinomio de grado  n que interpola f en n  1 puntos

diferentes x0 , x1 ,..., x n en a, b  . Para cada x en a, b  corresponde un punto  x en a, b tal que:

n
1
f ( x)  p ( x)  f ( n 1) ( x ) ( x  xk )
(n  1)! k 0

44 Leonardo Flores González

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

DIFERENCIAS DIVIDIDAS

n k 1
Sea p n ( xi )  c0   c  (x  x
k 1
k
j 0
j ), el polinomio que interpola f en n  1 puntos diferentes

x0 , x1 ,..., x n , se define la diferencia dividida cn  f x0 , x1 , x2 ,...xn  como el coeficiente de x n del

polinomio p n de a lo mas grado n que interpola f en x0 , x1 ,..., x n .

Teorema: Las diferencias divididas satisfacen la ecuación.

f x1.x2 ,...xn   f x0 .x1 ,...xn 1 

f x0 .x1 ,...xn  
xn  x0

A partir de este teorema se construye la tabla de diferencias divididas.

Propiedades:

1. La diferencia dividida es una función simétrica de sus argumentos.

n
2. Sea t un punto  a los nodos entonces f (t )  p(t )  f x0 , x1 ,..., xn , t   (t  x ) . k
k 0

3. Si f  C n a, b y x0 , x1 ,..., xn son puntos distintos de a, b  , entonces hay un punto  en a, b  ,

tal que: f x0 , x1 ,..., xn  

1 (n)
f ( ) .
n!

45 Leonardo Flores González

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

APLICACIONES

Ejemplo 1: Se desea construir el techo de un grifo con las siguientes características:

 El diámetro de la columna central es de 16.5 metros.
 Losa de la superficie cúbica de 0.1 metros de espesor.
 Cuatro vigas, la longitud de cada viga tiene la forma un polinomio cúbico de sección variable
entre 0.25 m. x 0.45 m. y 0.25m x 0.30 m., la viga situada en el plano XZ debe pasar por los
siguientes puntos:

x (metros) f(x) (metros) f'(x)

0 0 -
0.5 0.35 -
3 0.85 0

Tabla de diferencias
Tabla de diferencias

0.0000 0.0000 0.7000 -0.1667 0.0289

0.5000 0.3500 0.2000 -0.0800
3.0000 0.8500 0.0000
3.0000 0.8500

De acuerdo ha esto el polinomio de interpolación es:

P( x)  0.7( x  0.5)  0.1667 ( x  0.5)( x  3)  0.0289 ( x  5)(1  3) 2

luego los puntos calculados dibujan (las unidades se encuentran en metros):

X Z Zfinal Larco
0.0000 0.0000 4.0000 0.2529
0.2000 0.1549 4.1549 0.2412
0.4000 0.2897 4.2897 0.2313
0.6000 0.4058 4.4058 0.2231
0.8000 0.5047 4.5047 0.2166
1.0000 0.5878 4.5878 0.2114
1.2000 0.6563 4.6563 0.2075
1.4000 0.7118 4.7118 0.2047
1.6000 0.7555 4.7555 0.2028
1.8000 0.7889 4.7889 0.2015
2.0000 0.8133 4.8133 0.2007
2.2000 0.8302 4.8302 0.2003
2.4000 0.8410 4.8410 0.2001
2.6000 0.8469 4.8469 0.2000
2.8000 0.8495 4.8495 0.2000
3.0000 0.8500 4.8500
Longitud Total 3.1941

46 Leonardo Flores González

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

Polinomio Cúbico

1.0
0.8
0.6
f(x) 0.4
0.2
0.0
0.0 1.0 2.0 3.0 4.0
x

47 Leonardo Flores González

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

DIFERENCIAS FINITAS CENTRALES

La diferencia finita central es una función  :RR , definida por

f ( x)  f ( x  h / 2)  f ( x  h / 2), los nodos se encuentran igualmente espaciados de manera que
xk  x0  kh , para abreviar nomenclatura tomar en cuenta que f ( xi )  f i , de manera que

f i  f i 1 / 2  f i 1 / 2 .

Como f es una función lineal se tiene que  2 f i  f i 1  2 f i  f i 1 y

 3 f i  f i  3 / 2  3 f i 1 / 2  3 f i 1 / 2  f i  3 / 2 , se puede probar sin mayor dificultad que

n
 n fi   (1) k Ckn f n .
i  k
k 0 2

Los operadores f i  f i 1  f i y f i  f i  f i 1 son conocidos como diferencias finitas hacia

n
delante y diferencias finitas hacia atrás, ambos son generalizados como n f i   (1) C
k 0
k n
f
k i nk y

n
 n f i   (1) k Ckn f i  k respectivamente.
k 0

APLICACIONES

Las diferencias finitas centrales se utilizan para aproximar derivadas de funciones y también para
resolver ecuaciones diferenciales.

Por ejemplo la primera derivada puede ser expresada a partir de las series de Taylor:
h 2 ( 2 ) h 3 ( 3) h 2 ( 2 ) h 3 ( 3)
fi 1  fi  hfi (1)  fi  f (1 ) y fi 1  fi  hfi (1)  fi  f ( 2 ) de ambas fórmulas
2 6 2 6
resulta:
f i 1  f i 1 h 2 (3)
f i (1)   f (3 ) , análogamente se puede obtener una expresión para la segunda
2h 6
fi 1  2 f i  f i 1 h 2 ( 4)
derivada f i ( 2)   f ( ) .
h2 12
Ambas expresiones pueden emplearse para resolver: mx ( 2 )  cx (1)  kx   mu ( 2 ) (t ) , para valores

dados de m, c, k , u ( 2 ) (t ) .

48 Leonardo Flores González

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

La ecuación anterior puede ser escrita de la siguiente manera:

x ( 2)  2wx (1)  w 2 x  u ( 2) (t ) , con las siguientes condiciones iniciales x (0)  x' (0)  0.
Esta ecuación se puede escribir como una ecuación en diferencias de la siguiente manera:
xi 1  2 xi  xi 1 x  xi 1
 2w i 1  w 2 xi  ui( 2)
(t ) 2
2(t )
De la ecuación anterior obtenemos la siguiente sucesión:
1 2 1 w
xi 1  ( xi (  w 2 )  xi 1 (  )  ui )
 1 w  (t ) 2
(t ) 2
t
  
 (t )
2
t 

Si consideramos que x0  x1  0, t  0.02 y las aceleraciones las del registro sísmico de Lima de

1966, el cual será normalizado a 400 cm/seg 2.

Registro de Aceleraciones
1.5

0.5

0
0 10 20 30 40 50 60 70
-0.5

-1

-1.5

Se puede emplear el primer programa del Anexo I para obtener la respuesta del sistema y simular
la misma, queda como ejercicio para el lector elaborar otros algoritmos de integración directa para
obtener respuestas de sistemas de un grado de libertad.

49 Leonardo Flores González

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

Estabilidad del método de Newmark

Matriz de Amplificación: Se tiene el problema de valores y vectores propios A  B ,

donde A, B son matrices cuadradas de orden n sea la sucesión de vectores d definida por
n
Bd n 1  Ad n , cualquier vector de  n se puede expresar como d 0   c k  k , donde c k y  k son
k 1

constantes y vectores linealmente independientes. Empleando inducción matemática podemos

n
considerar d n  c  
k 1
k
n
k k de donde:

n n n
Ad n   ck nk A k  Bd n 1  B( c k nk 1 k )  d n 1   c k nk 1 k , de esta forma se confirma la
k 1 k 1 k 1

n
hipótesis inductiva d n  c   k
n
k k si además  k  1, k  1,2...n entonces lim d n  0. .
n 
k 1

A B
Teorema: Sean Q    una matriz cuadrada de orden 2n , A, B, C , D matrices diagonales
C D 
de orden n , si a j , b j , c j , d j son los elementos que ocupan la j-ésima posición en las diagonales

de A, B, C , D , entonces:

 H1  02x2 
a j bj 
det Q  det       det H 1  det H 2    det H n donde H j  
d j 
.
0 2 x 2  H n  c j

Por inducción supongamos que el teorema anterior se cumple para matrices de orden
2( n  1) , el determinante de Q lo podemos encontrar con el desarrollo de lo menores
complementarios de su n-ésima columna:

 An 1xn 1 Bn 1xn 1 0n 1x1   An 1xn 1 Bn 1xn 1 0n 1x1 

det Q  an det Cn 1xn 1 Dn 1xn 1 0n 1x1   (1) n cn det  0Tn 1x1 0Tn 1x1 bnn 
 0Tn 1x1 0Tn 1x1 d nn  Cn 1xn 1 Dn 1xn 1 0n 1x1 

donde An 1x1n 1 , Bn 1x1n 1 , C n 1x1n 1 , Dn 1x1n 1 son la matrices diagonales de orden n  1 con los

elementos en las diagonales a j , b j , c j , d j respectivamente, j  1...n  1, tomando nuevamente

menores complementarios tenemos:

50 Leonardo Flores González

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

A Bn 1xn 1  A Bn 1xn 1 
det Q  a n d n det  n 1xn 1  (1) 2 n 1 c n bn det  n 1xn 1
C n 1xn 1 Dn 1xn 1  C n 1xn 1 Dn 1xn 1 

agrupando términos:

A Bn 1xn 1 
det Q  det H n det  n 1xn 1  det H 1  det H 2    det H n
C n 1xn 1 Dn 1xn 1 

probando finalmente:

 H1  02x2 
det Q  det     
0 2 x 2  H n 

Estabilidad del método de Newmark: Para encontrar la solución de la ecuación diferencial lineal
homogénea:
Ma  Cv  f int  0

donde a, v, f int son el vector aceleración, el vector velocidad y el vector de fuerzas internas del

pórtico, se linealiza la ecuación anterior en cada paso de integración, la ecuación en diferencias

asociada a cada paso de integración es Ma n  Cv n  Kd n  0 , donde C  1M   2 K es la

matriz de amortiguamiento, los integradores consistentes son:

d n 1  d n  Δtvn  Δt 2  β a n  βa n 1 

v n 1  v n  Δt γ a n  γa n 1  .

donde   1   y   (1  2  ) / 2 . Para probar la estabilidad numérica del método de

Newmark, se multiplican los integradores anteriores por la matriz M, y se reemplaza
Ma n  Cv n  Kd n :

Md n 1  Md n  ΔtMv n  Δt 2 β Cv n  Kd n   β Cv n  Kd n 

Mv n 1  Mv n  γ Δt Cv n  Kd n   γΔt Cv n 1  Kd n 1 

expresando en forma matricial:
M  βΔt 2 K βΔt 2 C  d n 1  M  βΔt 2 K ΔtM  β Δt 2 C  d n 
    
 γΔtK M  γΔtC   v n 1    γΔtK M  γ ΔtC   v n 

considerando la matriz de amplificación:

M  β Δt 2 K ΔtM  βΔt 2 C  M  βΔt 2 K βΔt 2 C 
  z   μ  z 
  γ ΔtK M  γ ΔtC   γΔtK M  γΔtC 

51 Leonardo Flores González

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

   
donde z   d T v T  z1T z 2T , como el problema de valores y vectores propios K  M tiene n
vectores propios linealmente independientes, sea  la matriz de vectores propios, entonces los
vectores z1  a, z 2  b se pueden expresar como una combinación lineal de las columnas de

 ( a, b son vectores de constantes de orden n ), multiplicando la matriz de amplificación por la

matriz  se obtiene:

 0nxn   I  β Δt 2 L ΔtI  β Δt 2G  a   I  βΔt 2 L βΔt 2G  a 

 ,     μ  
0nxn     γ ΔtL I  γ ΔtG  b   γΔtL I  γΔtG  b 

ω12  0  2ξ1ω1  0 
  
donde: L      , y además G 
    
 0  ωn2   0  2ξ n ωn 
 
la penúltima expresión se puede expresar así: Hz  0 , donde:

 I  β Δt 2 L   ( I  βΔt 2 L) ΔtI  β Δt 2 G  μ( βΔt 2 G )  z1   A B   z1 

Hz           0    (c.25).
  γ ΔtL  μ(γΔtL) I  γ ΔtG  μ( I  γΔtG)   z 2  C D  z 2 

Se puede notar que las matrices A, B, C , D son matrices diagonales de orden n , los

elementos de sus diagonales son a j , b j , c j , d j , donde j  1, 2, ...n , de acuerdo al teorema

anterior, el determinante de la matriz H se puede expresar:

H1  0 
a j bj 
det H  det      , donde H j  
c j d j 
 0  H 2 

a j  1  β Δt 2 w 2j   (1  βΔr 2 ω jj )   (c.26)
b j  Δt  2 Δt 2 β ξω   (2 βΔt 2 ξω j )    (c.27)
c j   Δtγ w 2   ( Δtγω 2j )    (c.28)
d j  1  2γ Δtξω   (1  2γΔtξω j )    (c.29)

si los valores propios de H son menores que uno el algoritmo es estable, desarrollando el det H j
1
   (  0.25     2  2 )    (  0.25     2 2 )
se tiene:     si   1  t 
 ( 0.5   ) tw  (0.5   ) w j
 j 
para   0.5 y   0.25 la expresión anterior es incondicionalmente estable.

52 Leonardo Flores González

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

Capítulo 5

INTEGRACIÓN NUMÉRICA

53 Leonardo Flores González

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

CUADRATURA CERRADA DE NEWTON-COTES

b
Para Integrar I ( f )   f ( x)dx , la aproximación para este tipo de integrales se obtiene a partir de:
a

n n (x  x j ) f ( n 1) ( x ) n
f ( x)   f k    ( x  xk ) , integrando esta expresión:
k 0 j 0 ( xk  x j ) (n  1)! k  0
jk

b n b n (x  x j ) b
f ( n 1) ( x ) n
 f ( x)dx   f k   dx    ( x  xk )dx.
a k 0 a j 0
( xk  x j ) a
(n  1)! k  0
jk

La segunda integral expresa el error cometido al integrar, mientras que la primera integral expresa
el valor aproximado de la integral:
b n b n (x  x j )
 f ( x)dx   f k   dx sea x  a  ht , además los puntos de interpolación están
a k 0 a j 0
( xk  x j )
jk

ba
igualmente espaciados h  , sea los nodos xi  a  ih donde i  0,1,2..., n se tiene que:
n
b n (x  x j ) n n
(t  j )
n

  (x
a j 0 k  xj)
dx  h  
0 j 0
(k  j )
dt  h   k (t )dt.
0
jk jk

b b n


Sea  k   k (t)dt de esta expresión se tiene:  f ( x)dx  h f k k , para diferentes valores de n ,
a a k 0

n
por ejemplo 1 y 2 se tienen las reglas del trapecio y la de Simpson, se debe notar 
k 0
k  n que y

que además los  k son simétricos, existen también fórmulas de integración compuestas, de estas

reglas.
Una forma de representar el error es mediante la forma de Peano:
n b
Sea: R( f )  h  f k k   f ( x)dx .
k 0 a

54 Leonardo Flores González

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

Teorema: Sea R( P)  0 , para todo polinomio P  n . Luego para toda función f  C n 1 a, b :

b
R( f )   f ( n 1) (t ) K (t )dt
a

donde: K (t ) 
1
 n
 n ( x  t ) n si x  t
Rx x  t  , x  t   

.
n!  0 si x t

55 Leonardo Flores González

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

CUADRATURA DE GAUSS LEGENDRE

Método para integrar funciones con puntos que no están igualmente espaciados, este tipo de
cuadratura tiene mejor precisión que la cuadratura de Newton, el problema trata de resolver la
integral:
1 b
I1   h( x)dx , pero generalmente el problema es I 2   f ( x)dx , esta última integral se puede
1 a

ba ba
transformar a la primera con el siguiente cambio de variable , x   t    , de acuerdo
 2   2 

ba
1
ha este cambio de variable resulta I 2     g (t )dt donde g (t )  f ( b 2 a t  b 2 a ) .
 2 1

El polinomio de Lagrange que interpola a f en ( x1 , x2 ,..., xn ) se puede escribir así:

n n
p ( x)   f k   ( x k  xj j ) , se debe notar que los nodos no están igualmente espaciados, en general
xx

k 1 j 1
jk

se puede escribir a la función f como se indica a continuación:

n n n
f ( x)   f k   ( x k  xj j )  (
x x f (n)
( x ) x x j
n! xk  x j
) , si integramos esta expresión tendremos:
k 1 j 1 j 1
jk

1 n 1 n 1 n
f ( x)dx   f k   ( xk  x j )dx   (
x x j x x j

f (n)
( x )
n! xk  x j
)dx
1 k 1 1 j 1 1 j 1
jk

1 n 1 n

( f ( x)dx   f k Ak ,para este método los

x x j
donde: Ak 
1 j 1
xk  x j
)dx , luego podemos decir que 
1 k 1
jk

nodos son las raíces del polinomio de legendre de grado n.

Polinomio de Legendre

El Polinomio de Legendre se puede obtener de dos maneras, una a través fórmula de recurrencia:
2n  1 n 1
 n ( x)  x n 1 ( x)   n  2 ( x)
n n
 0 ( x)  1
1 ( x)  x

56 Leonardo Flores González

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

La otra forma es a partir de la ecuación diferencial ( x 2  1) y ' '2nxy ' n(n  1) y  0 , esta da como
solución la forma de Rodríguez del polinomio de Legendre:

1 d n ( x 2  1)n
 n ( x)  .
2n n! dx n

Estos polinomios cumplen con las siguientes propiedades:

1. Sus raíces son simétricas con respecto al origen y se encuentran en el intervalo  1,1  .
1
2. Son ortogonales es decir:   ( x )
1
i j ( x)dx  0 si i  j . Para esto ver (Earl A. Coddington

“Introducción a las ecuaciones diferenciales”, Editorial C.E.C.S.A, 1940).

3. De la propiedad 2 Se puede ver que los polinomios de Legendre  0 ( x),  0 ( x),...,  n ( x) son

linealmente independientes, luego cualquier polinomio de grado n se puede escribir así:

n
p( x)   ck  k ( x) .
k 0

1
4. Si q (x) es un polinomio de grado menor o igual a n-1 se cumple que  q ( x )
1
n ( x)dx  0 .

TEOREMA. Sean ( x1 ,..., xn ) las raíces del polinomio de Legendre de grado n, entonces se puede

integrar de forma exacta polinomios de grado menor o igual a 2n-1 con la siguiente fórmula de
cuadratura:

1 n 1 n

 p( x)dx   p A (
x x j
k k donde Ak  xk  x j
)dx .
1 k 1 1 j 1
jk

Si p (x) es de grado menor o igual a n-1, entonces la fórmula es exacta ya que el error de
interpolación viene dado por:
n

(
p ( n ) ( x ) x x j
E ( x)  n! xk  x j
)  0.
j 1

1 n
Lo que indica que la integración es exacta, luego  p( x)dx   pk Ak .
1 k 1

Si p (x) es de grado menor o igual a 2n-1 y de grado mayor o igual a n, se puede decir
que p ( x)   n ( x)q ( x)  r ( x) , integrando esta última expresión se tiene:

57 Leonardo Flores González

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

1 1 1 1 n

 p( x)dx    n ( x)q( x)dx   r ( x)dx   r ( x)dx   rk Ak .

1 1 1 1 k 1

1 n
Además p ( xk )  r ( xk ) , luego  p( x)dx   pk Ak .
1 k 1

Para el cálculo de los Ak se emplea el teorema anterior:

Si m es menor o igual que n-1 se tiene que:

1  (1) m  2
1

 x dx  A1 x1 A2 x2  A3 x3      An xn 
m m m m m
, esto determina el sistema de ecuaciones:
1
m 1

 1 1  1  A1   2 
    
 x1 x2    xn  A2   0 
 La solución de este sistema nos da los Ak
          
 n 1    1 ( 1) n1 
x n 1
 xnn 1  An   n 
 1 x 2

De lo anterior podemos deducir el siguiente algoritmo:

Algoritmo
b
1. Ingresar f , a, b, n , para integrar I   f ( x)dx.
a

1
2. Convertir I  b 2 a  g (t )dt donde g (t )  f ( b 2 a t  b 2 a ) .
1

3. Encontrar las raíces de  n (x) , estas son ( x1 ,..., xn ) .

4. Resolver el sistema:
 1 1  1  A1   2 
    
 x1 x2    xn  A2   0 

          
 n 1    1 ( 1) n1 
x n 1
 xnn 1  An   n 
 1 x
2

n
5. I  ba
2 g
k 1
k Ak .

58 Leonardo Flores González

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

BIBLIOGRAFÍA

Bathe Klaus-Jurgen, “Finite Element Procedure in Engineering Análisis”, Prentice Hall, Inc.,
Englewood Cliffs, New Yersey, Estados Unidos de Norte América,1982.

Burden Richard, Faires Douglas, “Análisis Numérico” Internacional Thomson Editores, sexta
edición, México, 1988.

Coddington Earl, “Introducción a las Ecuaciones Diferenciales” C.E.C.S.A., segunda edición,

México, 1940.

Kincaid D., Cheney W., “Análisis Numérico” Addison Wesley Iberoamericana, sexta edición,
México, 1988.

Stoer J., Bulirsch R., “Introduction to Numérical Analysis”, Springer Verlag, segunda edición,
Estados Unidos de Norte América,1980.

59 Leonardo Flores González

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

ANEXO

Como anexo se entregan los siguientes programas:

Anexo I
El primer Programa está referido a la solución de un sistema vibratorio de un grado de libertad, el
cual es útil para el ingreso de registros sísmicos y simula la respuesta del sistema, el archivo
[Link].

load previo
caras1=[1 2 3 4;1 2 6 5;2 3 7 6;4 3 7 8;1 4 8 5;5 6 7 8];
ver1=[-1 -1 0;-1 1 0;1 1 0;1 -1 0;-1 -1 10;-1 1 10;1 1 10;1 -1 10];
ver2=[-3 -3 10;-3 3 10;3 3 10;3 -3 10;-3 -3 13;-3 3 13;3 3 13;3 -3 13];
acel=sacel';
u=t';
uu=t';
uuu=t';
figure(4)
clf
subplot(2,2,1)
hdl1=patch('faces',caras1,'vertices',[ver1(:,1) ver1(:,2)
ver1(:,3)],'FaceColor','r');
hdl2=patch('faces',caras1,'vertices',[ver2(:,1) ver2(:,2)
ver2(:,3)],'FaceColor','g');
hold on
set(gcf,'renderer','zbuffer');
grid on
view(30,30)
title('Sistema Masa-Resorte');
xlabel('X')
ylabel('Y')
zlabel('Z')
axis('equal')
axis('square')
axis([-10 10 -10 10 -1 15])
subplot(2,2,2)
grid on
ef=plot(0,0);
axis([min(u) max(u) min(acel) max(abs(acel))]);
grid on

60 Leonardo Flores González

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

title('Registro Sísmico');
xlabel('tiempo');
ylabel('aceleraciones del registro');
subplot(2,2,3)
ff=plot(0,0);
axis([min(uu) max(uu) min(x) max(abs(x))]);
grid on
title('desplazamiento vs. tiempo');
xlabel('tiempo');
ylabel('desplazamientos');
subplot(2,2,4)
gf=plot(0,0); %set(gca, 'Color', 'r');
axis([min(uuu) max(uuu) min(xx) max(abs(xx))]);
grid on
title('velocidades vs. tiempo');
xlabel('tiempo');
ylabel('velocidades');
for i=[Link]length(acel);
a=x(:,i);
aa=2*a;
ttt=[-3+aa -3 10;-3+aa 3 10;3+aa 3 10;3+aa -3 10;-3+aa -3 13;-3+aa 3 13;3+aa 3
13;3+aa -3 13];
ppp=[-1 -1 0;-1 1 0;1 1 0;1 -1 0;-1+aa -1 10;-1+aa 1 10;1+aa 1 10;1+aa -1 10];
set(hdl2,'vertices',[ttt(:,1) ttt(:,2) ttt(:,3)]);
set(hdl1,'vertices',[ppp(:,1) ppp(:,2) ppp(:,3)]);
set(ef,'xdata',u(:,1:i),'ydata',acel(:,1:i));
set(ff,'xdata',uu(:,1:i),'ydata',x(:,1:i));
set(gf,'xdata',uuu(:,1:i),'ydata',xx(:,1:i));
drawnow
end

61 Leonardo Flores González