UNIVERSIDAD NACIONAL DE INGENIERÍA

FACULTAD DE INGENIERÍA CIVIL

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   2.32.Log (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('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*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('pfijo.xls','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('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:

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