Iniciacion a la Investigaci on en Matematica Aplicada

DEPARTAMENTO DE MATEMÁTICAS UNIVERSIDAD DE EXTREMADURA 06071-BADA JOZ (Spain)

Universidad de ExtremaduraDepartamento de Matematicas

Iniciacion a la Investigacion enMatematica Aplicada

Pedro Martın Jimenez

Badajoz, noviembre de 2014

Indice

1. Errores, redondeo, estabilidad, condicionamiento. 71.1. Cifras significativas. Exactitud y precision. Errores. . . . . . . . 71.2. Calculos estables e inestables. Condicionamiento. . . . . . . . . 10

1.2.1. Inestabilidad . . . . . . . . . . . . . . . . . . . . . . . . 101.2.2. Condicionamiento . . . . . . . . . . . . . . . . . . . . . 10

1.3. Aritmetica de la computadora . . . . . . . . . . . . . . . . . . . 111.3.1. Aritmetica de punto flotante. . . . . . . . . . . . . . . . 111.3.2. Operaciones con computadoras . . . . . . . . . . . . . . 131.3.3. Epsilon de la maquina . . . . . . . . . . . . . . . . . . . 14

2. Resolucion de ecuaciones no lineales 152.1. Metodo de la biseccion . . . . . . . . . . . . . . . . . . . . . . . 17

2.1.1. Descripcion del metodo . . . . . . . . . . . . . . . . . . 172.1.2. Convergencia del metodo . . . . . . . . . . . . . . . . . 172.1.3. Aproximacion y error . . . . . . . . . . . . . . . . . . . 182.1.4. Variaciones del metodo: Regula Falsi . . . . . . . . . . . 19

2.2. Metodo de Newton-Raphson . . . . . . . . . . . . . . . . . . . . 202.2.1. Descripcion del metodo . . . . . . . . . . . . . . . . . . 202.2.2. Convergencia del metodo . . . . . . . . . . . . . . . . . 212.2.3. Aproximacion y error . . . . . . . . . . . . . . . . . . . 222.2.4. Variaciones del metodo: Metodo de la secante . . . . . . 25

2.3. Metodo iterativo de punto fijo . . . . . . . . . . . . . . . . . . . 262.3.1. Descripcion del metodo . . . . . . . . . . . . . . . . . . 262.3.2. Convergencia del metodo . . . . . . . . . . . . . . . . . 272.3.3. Aproximacion y error . . . . . . . . . . . . . . . . . . . 28

2.4. Raıces de polinomios . . . . . . . . . . . . . . . . . . . . . . . . 302.4.1. Separacion de raıces. Sucesion de Sturm. . . . . . . . . . 322.4.2. Acotacion de raıces . . . . . . . . . . . . . . . . . . . . . 352.4.3. Raıces de polinomios con el algoritmo de Horner . . . . 36

3

4 INDICE

2.4.4. Raıces multiples . . . . . . . . . . . . . . . . . . . . . . 37

3. Sistemas Lineales 39

3.1. Algebra de matrices . . . . . . . . . . . . . . . . . . . . . . . . 39

3.1.1. Valores propios y vectores propios . . . . . . . . . . . . 433.1.2. Matriz definida . . . . . . . . . . . . . . . . . . . . . . . 43

3.2. Resolucion de sistemas de ecuaciones lineales: metodo de Gauss 44

3.2.1. Metodo de Gauss . . . . . . . . . . . . . . . . . . . . . . 45

3.2.2. Metodo de Gauss con pivoteo . . . . . . . . . . . . . . . 47

3.3. Factorizacion LU. Factorizacion de Cholesky . . . . . . . . . . . 51

3.3.1. Metodo de Crout . . . . . . . . . . . . . . . . . . . . . . 533.3.2. Metodo de Cholesky . . . . . . . . . . . . . . . . . . . . 53

3.3.3. Sistemas triangulares . . . . . . . . . . . . . . . . . . . . 53

3.4. Normas y analisis del error . . . . . . . . . . . . . . . . . . . . 55

3.4.1. Numero condicion de una matriz . . . . . . . . . . . . . 58

3.5. Mejora de soluciones . . . . . . . . . . . . . . . . . . . . . . . . 603.5.1. Refinamiento iterativo . . . . . . . . . . . . . . . . . . . 60

3.5.2. Escalamiento . . . . . . . . . . . . . . . . . . . . . . . . 61

3.6. Metodos iterativos . . . . . . . . . . . . . . . . . . . . . . . . . 62

3.6.1. Metodo de Jacobi . . . . . . . . . . . . . . . . . . . . . 64

3.6.2. Metodo de Gauss-Seidel . . . . . . . . . . . . . . . . . . 66

3.6.3. Metodos de relajacion . . . . . . . . . . . . . . . . . . . 67

4. Aproximacion de funciones 69

4.1. Aproximacion alrededor de un punto . . . . . . . . . . . . . . . 69

4.1.1. Polinomio de Taylor . . . . . . . . . . . . . . . . . . . . 69

4.1.2. Teorema de Taylor. Formula de Taylor . . . . . . . . . . 72

4.2. Aproximacion simultanea en varios puntos . . . . . . . . . . . . 744.2.1. Polinomio interpolador . . . . . . . . . . . . . . . . . . . 74

4.3. Error del polinomio interpolador . . . . . . . . . . . . . . . . . 77

4.3.1. Eleccion de nodos. Polinomios de Chebyshev . . . . . . 78

4.4. Interpolacion a trozos y con condiciones sobre la derivada . . . 81

4.4.1. Interpolacion a trozos . . . . . . . . . . . . . . . . . . . 81

4.4.2. Interpolacion con condiciones sobre la derivada . . . . . 824.5. Aproximacion por mınimos cuadrados . . . . . . . . . . . . . . 90

4.6. Aproximacion de funciones periodicas. Series de Fourier . . . . 92

4.6.1. Funciones periodicas. Funciones pares e impares . . . . 92

4.6.2. Desarrollo de una funcion en serie de Fourier . . . . . . 94

INDICE 5

4.7. Practicas con SAGE . . . . . . . . . . . . . . . . . . . . . . . . 96

5. Diferenciacion e integracion numerica 995.1. Diferenciacion numerica y extrapolacion de Richarson . . . . . 99

5.1.1. Diferenciacion mediante interpolacion . . . . . . . . . . 1025.1.2. Extrapolacion de Richarson . . . . . . . . . . . . . . . . 104

5.2. Integracion numerica mediante interpolacion . . . . . . . . . . . 1085.2.1. Regla del trapecio . . . . . . . . . . . . . . . . . . . . . 1085.2.2. Regla de Simpson . . . . . . . . . . . . . . . . . . . . . 1105.2.3. Reglas compuestas . . . . . . . . . . . . . . . . . . . . . 1125.2.4. Metodo de los coeficientes indeterminados . . . . . . . . 114

5.3. Cuadratura gaussiana . . . . . . . . . . . . . . . . . . . . . . . 1155.4. Integracion de Romberg . . . . . . . . . . . . . . . . . . . . . . 1165.5. Cuadratura adaptativa . . . . . . . . . . . . . . . . . . . . . . . 1195.6. Practicas con SAGE . . . . . . . . . . . . . . . . . . . . . . . . 122

6. Resolucion de ecuaciones diferenciales ordinarias de primerorden 1256.1. Existencia y unicidad de soluciones . . . . . . . . . . . . . . . . 1256.2. Resolucion numerica de EDO de primer orden . . . . . . . . . . 128

6.2.1. Metodo de la serie de Taylor . . . . . . . . . . . . . . . 1286.2.2. Errores . . . . . . . . . . . . . . . . . . . . . . . . . . . 1326.2.3. Metodos de Runge-Kutta . . . . . . . . . . . . . . . . . 1336.2.4. Metodos multipaso . . . . . . . . . . . . . . . . . . . . . 139

6.3. Resolucion exacta de EDO de primer orden . . . . . . . . . . . 1476.3.1. Ecuaciones en variables separadas . . . . . . . . . . . . 1496.3.2. Ecuaciones diferenciales homogeneas . . . . . . . . . . . 1506.3.3. Ecuaciones reducibles a homogeneas . . . . . . . . . . . 1516.3.4. Ecuaciones diferenciales lineales . . . . . . . . . . . . . . 1526.3.5. Ecuaciones de Bernoulli . . . . . . . . . . . . . . . . . . 154

6.4. Sistemas de primer orden y ecuaciones diferenciales de ordensuperior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1546.4.1. Sistemas de ecuaciones de primer orden . . . . . . . . . 1546.4.2. Ecuaciones de orden superior . . . . . . . . . . . . . . . 156

6.5. Practicas con SAGE . . . . . . . . . . . . . . . . . . . . . . . . 157

Problemas 159

Capıtulo 1

Errores, redondeo,estabilidad, condicionamiento.

1.1. Cifras significativas. Exactitud y precision. Erro-res.

Cualquier numero real x puede representarse en forma decimal con unnumero finito o infinito de dıgitos:

x = αmαm−1 . . . α1α0 . β1β2 . . . βn . . . con αi, βj ∈ 0, 1, 2, . . . , 9

Con la expresion anterior queremos representar que

x = αm10m+αm−110m−1+· · ·+α110+α0100+β110−1+β210−2+· · ·+βn10−n+. . .

El Calculo Numerico consiste en desarrollar metodos para encontrar solucioneslo mas aproximadas posibles a la verdadera solucion de un cierto problema.Se trabaja, por tanto, con numeros aproximados en lugar de numeros exactos,por lo que se utilizan aproximaciones con un numero finito de decimales.

Toda expresion aproximada tiene un numero determinado de cifras signi-ficativas. El concepto de cifra significativa intenta trasmitir la idea de cuandouna cifra que aparece en la expresion decimal trasmite informacion .esencial 2

”de fiar”sobre el numero que se intenta aproximar. La significacion de unacifra depende del contexto y no solo de la expresion decimal en la que aparece.

Los ceros que se anaden a la izquierda de una expresion decimal para situarla primera cifra distinta de cero no son significativos:

7

8 Capıtulo 1

Ejemplo: Los numeros

0,00001845 0,0001845 0,001845

tienen todos cuatro cifras significativas.

Los ceros que se anaden al final pueden ser o no significativos:Ejemplo: En la expresion ”se manifestaron 45000 personas”, los tres ulti-

mos ceros parecen ser no significativos. Sin embargo en la expresion ”mecosto 45000 euros exactos”sı son significativos.

Para evitar equıvocos acerca del numero de cifras significativas, se emplea

la notacion cientıfica (en el ejemplo anterior, 4,5× 104, 4,5000× 104).

Cualquier metodo del Calculo Numerico que se emplee en la resolucionde un problema, debe ser suficientemente exacto, es decir, los valores apro-ximados que se obtengan en diferentes intentos de solucion deben estar cercadel valor verdadero. Cuanto mayor sea el numero de cifras significativas de laaproximacion, mayor sera la exactitud. Tambien deben ser preciso, es decirlos diferentes valores obtenidos para la resolucion de un problema deben estarcercanos entre sı.

Con el uso de aproximaciones para representar operaciones y cantidadesmatematicas se generan errores numericos. Se denomina error E a la can-tidad:

E = valor verdadero− valor aproximado

A menudo se trabaja con el error absoluto (|E|). El error relativo es

e =E

valor verdadero.

Este ultimo compara la magnitud del error cometido con la magnitud del valorque se pretende estimar y puede interpretarse en terminos de %.

Las causas de los errores pueden ser:• Truncamiento: surgen al sustituir un procedimiento matematico exacto

por uno aproximado. Ejemplo: utilizar un polinomio de Taylor para encontrarel valor de una funcion en un punto.• Redondeo: surgen al sustituir un numero exacto por uno aproximado.

Ejemplo: sustituir 1/3 por 0,33333333. Para minimizar este tipo de error serecurre a las reglas de redondeo:

Errores, redondeo, estabilidad, condicionamiento 9

1. En los calculos, se conservan las cifras significativas y el resto se des-cartan. El ultimo dıgito que se conserva se aumenta en uno si el primerdıgito descartado es mayor que 5. Si es 5 o es 5 seguido de ceros, enton-ces el ultimo dıgito retenido se incrementa en uno solo si este ultimo esimpar:

numero 6 cifras significativas 8 cifras significativas

5.6170431500 5.61704 5.6170432

5.6170462500 5.61705 5.6170462

2. En la suma y en la resta, el redondeo se lleva a cabo de forma tal que elultimo dıgito retenido en la respuesta corresponda al ultimo dıgito MASsignificativo de los numeros que se estan sumando o restando. Noteseque un dıgito en la columna de las centesimas es mas significativo queuna en la columna de las milesimas.

2,2− 1,768 = 0,432 −→ 0,4

4,68× 10−7 + 8,3× 10−4 − 228× 10−6 = 6,02468× 10−4 −→ 6,0× 10−4

3. En la multiplicacion y division se conserva el numero mınimo de cifrassignificativas que tenga los numeros que intervienen:

0,0642× 4,8 = 0,30816 −→ 0,31945/0,3185 = 2967,032967 −→ 297× 10

4. Para combinaciones de las operaciones aritmeticas, existen dos casosgenerales.

(multiplicacion o division )± (multiplicacion o division)

(suma o resta )×÷ (suma o resta)

En ambos casos, se ejecutan primero las operaciones entre parentesis yse redondea el resultado antes de proceder con otra operacion.

• Otras causas de error: equivocacion de usuario, mala formulacion delmodelo, incertidumbre en los datos fısicos recogidos, etc.

10 Capıtulo 1

1.2. Calculos estables e inestables. Condicionamien-to.

1.2.1. Inestabilidad

Decimos que un proceso numerico es inestable cuando los pequenos erro-res que se producen en alguna de sus etapas se agrandan en etapas posterioresy degradan seriamente la exactitud del calculo en su conjunto.

Ejemplo: Consideremos la sucesion

x0 = 1, x1 =1

3, xn+1 =

13

3xn −

4

3xn−1 ∀n ≥ 1

La sucesion anterior verifica que xn = (13)n. Si calculamos el termino x11 conel programa Maxima, obtenemos −1,086162170293413×10−6 con la definicioninductiva y 2,867971990792441 × 10−10 con la otra definicion. El calculo delos valores de la sucesion por el metodo inductivo es inestable, pues cualquiererror que se presente en xn se multiplica por 13/3 al calcular xn+1.

1.2.2. Condicionamiento

Un problema esta mal condicionado si pequenos cambios en los datos pue-den dar lugar a grandes cambios en las respuestas. En ciertos problemas sepuede calcular un numero de condicion. Si este es grande significa que el pro-blema esta mal condicionado.

Ejemplo: Dado un sistema de ecuaciones lineales en forma matricial Ax =b, se puede calcular el numero de condicion ası:

k(A) = ∥A∥ · ∥A−1∥,

donde ∥A∥ = maxi,j |aij |. Dado el sistema

x + 1,01y = 10,99x + y = 1

el numero de condicion de la matriz del sistema serıa k(A) = ∥A∥ · ∥A−1∥ =(1,01)(10100) = 10201. Esto significa que pequenos cambios en b producirangrandes cambios en la solucion. De hecho:

x + 1,01y = 10,99x + y = 1

⇒ [x = −100, y = 100]


yx + 1,01y = 1,010,99x + y = 1

⇒ [x = 0, y = 1]

1.3. Aritmetica de la computadora

1.3.1. Aritmetica de punto flotante.

Las computadoras utilizan la llamada aritmetica del punto flotante, es de-cir, almacenan cada numero desplazando sus dıgitos de modo que aparezca dela forma

0.β1β2...βn × 10e.

Por ejemplo:

13,524→ 0,13524× 102

−0,0442→ −0,442× 10−1

Como todos los numeros se traducen segun la aritmetica del punto flotante,todos empezaran con 0. · · · y terminaran con 10e. El ordenador, para almace-narlos ahorrando memoria, prescinde de lo que es comun a todos los numeros,y ası el 13,524 se almacena como 13524E2 y el −0,0442 como −442E − 1.

Por otra parte, la mayor parte de las computadoras trabajan con numerosreales en sistema binario, en contraste con el sistema decimal que normalmentese utiliza. Ası, el numero 9.90625 en sistema decimal, se convierte en binarioen:

9 = 2 ∗ (2 ∗ (2 ∗ 1 + 0) + 0) + 1 = 1 ∗ 23 + 0 ∗ 22 + 0 ∗ 21 + 1 ∗ 20 = (1001)2

0,90625 ∗ 2 = 1,81250 1 ∗ 2−1

0,81250 ∗ 2 = 1,625 1 ∗ 2−2

0,625 ∗ 2 = 1,25 1 ∗ 2−3

,25 ∗ 2 = 0,5 1 ∗ 2−4

,5 ∗ 2 = 1 1 ∗ 2−5

⇒0,90625 = 1 ∗ 2−1 + 1 ∗ 2−2 + 1 ∗ 2−3 + 0 ∗ 2−4 + 1 ∗ 2−5 = (0,11101)2

Con lo que

1,90625 = 1 ∗ 23 + 0 ∗ 22 + 0 ∗ 21 + 1 ∗ 20 + 1 ∗ 2−1 +

+1 ∗ 2−2 + 1 ∗ 2−3 + 0 ∗ 2−4 + 1 ∗ 2−5 = (1001,11101)2

12 Capıtulo 1

En general, la representacion de un numero en una computadora sera dela forma:

±′d1d2...dp ×Be

y se almacenara de la forma

±d1d2...dpEe

donde B es la base, di ∈ 0, 1, 2, ..., B − 1, d1 > 0, p es el numero maximode bits significativos y e es el exponente de la potencia de B. El 0 es un casoespecial que se almacena como 00..,0E0. La parte d1d2...dp se llama partesignificativa o mantisa. p es el numero maximo de bits que puede tener esaparte significativa del numero que se almacena y es, por tanto, finita. Cuantomas grande sea p, mayor sera la precision de la maquina. Como consecuencia detodo lo anterior, en cualquier maquina solo se pueden almacenar una cantidadfinita de numeros, llamados numeros maquina. En los ordenadores de base 2,d1 siempre es 1, con lo que se puede suprimir, aumentando la precision. Segana con ello un bit que se denomina bit escondido.

Ejemplo: Si en una computadora se tiene que B = 2, p = 2 y −3 ≤ e ≤ 3,entonces solo se podran almacenar numeros de la forma

±,102 × 2e ± ,112 × 2e con − 3 ≤ e ≤ 3.

Como,102 = 1× 2−1 + 0× 2−2 = 1/2

los numeros ±,102 × 2e son

±4,±2,±1,±1/2,±1/4,±1/8,

y como,112 = 1× 2−1 + 1× 2−2 = 1/2 + 1/4 = 3/4

los numeros ±,112 × 2e son

±6,±3,±3/2,±3/4,±3/8,±3/16.

En este sistema el numero 2′9 se almacenarıa, segun este disenada lacomputadora, bien como 2 si se hace por truncamiento, o bien como 3 sise hace por redondeo. Tambien se almacenarıan como 2 los numeros 2′3 o 2′4.


Los parametros B, p y e de una maquina real pueden ser los siguientes:B = 2, p = 23 bits, e = 8 bits, con lo cual el numero maximo que se podrıaalmacenar serıa 1′701E38. Todo numero mayor que el anterior estarıa en eldesbordamiento positivo de la maquina. El numero mınimo en valor absolutoserıa 1′755E−38. Todo numero positivo menor estarıa en el subdesbordamientopositivo de la maquina. De forma similar se definen desbordamiento negativoy subdesbordamiento negativo.

1.3.2. Operaciones con computadoras

Supongamos una computadora cuyos parametros sean B = 10, p = 3 y−9 ≤ e ≤ 9. Veamos como trabajan y como se pueden generar errores almanejar las operaciones basicas y el redondeo:

1. Para la suma, el de menor exponente se desplaza para alinear la comadecimal

1,37 + 0,0269 = 0,137× 101 + 0,269× 10−1 =

= 0,137× 101 + 0,00269× 101 = 0,13969× 101

que se almacenarıa como 0,139×101 por truncamiento y 0,140×101 porredondeo.

2. En la resta:

4850− 4820 = 0,485× 104 − 0,482× 104 = 0,003× 104 = 0,300× 102

que se almacenarıa como 0,300×102 tanto por redondeo como por trun-camiento.

3. En el producto, se multiplican las cifras significativas y se suman losexponentes:

403000 · 0,0197 = 0,403× 106 · 0,197× 10−1 = 0,079391× 105

que se almacena como 0,793 × 104 por truncamiento y 0,794 × 104 porredondeo.

4. En la division, se dividen las cifras significativas y se restan los exponen-tes:

0,0356/1560 = 0,356× 10−1/0,156× 104 = 2,28205× 10−5

que se almacena como 0,228× 10−4.

14 Capıtulo 1

Otra fuente de error puede ser el paso de sistema decimal a sistema binario.Ası el numero 1/10 en sistema binario es (0,0001100110011001...)2.

1.3.3. Epsilon de la maquina

La diferencia mas pequena entre dos numeros que puede detectar unamaquina se denomina epsilon de la maquina. Tambien se puede definir comoel menor numero ε > 0 tal que

1,0 + ε = 1,0

y se puede calcular programando este algoritmo:

input s← 1,0for k = 1, 2, . . . , 100 do

s← 0,5st← s + 1,0

if t ≤ 1,0 thens← 2,0s

output k-1,sstop

end ifend

Capıtulo 2

Resolucion de ecuaciones nolineales

-10 -8 -6 -4 -2

-1

-0.5

0.5

1

1.5

f(x) = ex − sen(x)

Dada una funcion f : R → R, trata-remos de encontrar metodos que permitanlocalizar valores aproximados de las solu-ciones de f(x) = 0.

Ejemplo: Dada la ecuacion ex −sen(x) = 0, encuentra la solucion mas cer-cana a 0.

Para resolver ejercicios como el anterior es conveniente recordar los siguien-

tes resultados:

Teorema [Bolzano]: Sea f : R → R una funcion continua en [a, b] talque f(a)f(b) < 0. Entonces existe c ∈ (a, b) tal que f(c) = 0.

El teorema de Bolzano se utiliza en este contexto para localizar intervalosque contengan una solucion de f(x) = 0 cuando la funcion f(x) sea continua.

Ejercicio: Dada la ecuacion ex − sen(x) = 0, localiza un intervalo dondeexista una solucion.

15

16 Capıtulo 2

Teorema [Rolle]: Sea f : R→ R una funcion continua en [a, b], derivableen (a, b) y tal que f(a) = f(b). Entonces existe c ∈ (a, b) tal que f ′(c) = 0.

Del teorema de Rolle se deduce que, si f(x) cumple las hipotesis, entrecada dos soluciones de la ecuacion f(x) = 0 debe existir al menos una solucionde f ′(x) = 0. De este modo, el numero maximo de soluciones de f(x) = 0sera el numero de soluciones de f ′(x) = 0 mas uno. Por tanto, el teorema deRolle se utilizara para determinar el numero maximo de soluciones que puedetener una ecuacion f(x) = 0.

Ejercicio: Determina el numero de soluciones de la ecuacion ex − x = 0.

Teorema [Taylor]: Sea f : R→ R una funcion n veces derivable en [a, b]y n + 1 veces derivable en (a, b). Entonces, para cada x0 ∈ [a, b] existe unpolinomio de grado menor o igual que n que tiene un punto de contacto conf(x) de grado n en x0. Dicho polinomio es

Pn,x0(x) = f(x0) + f ′(x0)(x− x0) + · · ·+ fn)(x0)

n!(x− x0)

n.

Ademas, para cada x ∈ [a, b] existe un punto θx entre x y x0 tal que

f(x)− Pn,x0(x) =fn+1)(θx)

(n + 1)!(x− x0)

n+1.

Teorema [de valor medio]: Sea f : R → R una funcion continua en[a, b] y derivable en (a, b). Entonces existe c ∈ (a, b) tal que

f(b)− f(a)

b− a= f ′(c).

Resolucion de ecuaciones no lineales 17

2.1. Metodo de la biseccion

2.1.1. Descripcion del metodo

Consiste en aplicar de forma reiterada el teorema de Bolzano en intervaloscada vez mas pequenos. Supongamos que f(a) < 0 y f(b) > 0 (de formasimilar se razonarıa si f(a) > 0 y f(b) < 0).

Paso 0. Partimos del intervalo [a, b] y calculamos c0 = a+b2 .

Paso 1. Si f(c0) < 0, definimos a1 = c0 y b1 = b. Si f(c0) > 0, definimosa1 = a y b1 = c0. En ambos casos, estudiamos el intervalo [a1, b1] y calculamosc1 = a1+b1

2 .Paso 2. Si f(c1) < 0, definimos a2 = c1 y b2 = b1. Si f(c1) > 0, definimos

a2 = a1 y b2 = c1. En ambos casos, estudiamos el intervalo [a2, b2] y calculamosc2 = a2+b2

2 .Paso n. Si f(cn−1) < 0, definimos an = cn−1 y bn = bn−1. Si f(cn−1) > 0,

definimos an = an−1 y bn = cn−1. En ambos casos, estudiamos el intervalo[an, bn] y calculamos cn = an+bn

2 .

2.1.2. Convergencia del metodo

Teorema: Si f(x) es una funcion continua en [a, b] y [a, b], [a1, b1], . . . , [an, bn], . . .denotan los intervalos utilizados en el metodo de la biseccion, entonces

lımn→∞

an = lımn→∞

bn = lımn→∞

cn = r

y ademas f(r) = 0.

Demostracion. Por la construccion de las sucesiones se tiene que

a ≤ a1 ≤ a2 ≤ · · · ≤ b2 ≤ b1 ≤ b

y

(bn − an) =1

2(bn−1 − an−1) = · · · = 1

2n(b− a).

Por tanto, la sucesion (an)n∈N converge por ser creciente y acotada superior-mente y la sucesion (bn)n∈N converge por ser decreciente y acotada inferior-mente. Veamos que sus lımites son iguales

lımn→∞

bn − lımn→∞

an = lım(bn − an) = lımn→∞

1

2n(b− a) = 0.

18 Capıtulo 2

Dado lo anterior y puesto que an ≤ cn ≤ bn se concluye que

lımn→∞

an = lımn→∞

cn = lımn→∞

bn.

Llamemos a ese lımite r. Veamos ahora que r es la solucion de f(x) = 0utilizando las hipotesis de continuidad de f(x) y la condicion f(an)f(bn) ≤ 0 :

0 ≤ f(r)f(r) = f( lımn→∞

an)f( lımn→∞

bn) =

lımn→∞

f(an) lımn→∞

f(bn) = lımn→∞

f(an)f(bn) ≤ 0.

Por tanto f(r) = 0 que era lo que se tenıa que probar.

2.1.3. Aproximacion y error

Una aproximacion a la solucion de la ecuacion f(x) = 0 en el paso n es cn.Una cota del error cometido sera:

|r − cn| ≤1

2|bn − an| =

1

22|bn−1 − an−1| =

1

2n+1|b− a|.

Ejemplo: Encuentra una aproximacion de la solucion de ex = sen(x) conun error menor que 10−2.

Aplicaremos el metodo de la biseccion a la funcion f(x) = ex− sen(x). Porel teorema de Bolzano sabemos que hay una solucion de f(x) = 0 en el intervalo[−4,−3], puesto que f(−4) < 0 y f(−3) > 0. Si cn es la aproximacion y r es lasolucion exacta, como queremos que el error sea menor que 10−2, tendremosque:

|r − cn| ≤1

2n+1|b− a| = 1

2n+1| − 3− (−4)| = 1

2n+1< 10−2

por lo que n ≥ 6, es decir, tendremos que calcular c6.Los calculos serıan los siguientes:

a0 = −4 b0 = −3 c0 = −3,5 f(c0) = −0,3206a1 = −3,5 b1 = −3 c1 = −3,25 f(c1) = −0,0694a2 = −3,25 b2 = −3 c2 = −3,125 f(c2) = 0,0554a3 = −3,25 b3 = −3,125 c3 = −3,1875 f(c3) = −0,0046a4 = −3,1875 b4 = −3,125 c4 = −3,15625 f(c4) = 0,0277a5 = −3,1875 b5 = −3,15625 c5 = −3,175 f(c5) = 0,0084a6 = −3,1875 b6 = −3,175 c6 = −3,18125

La solucion aproximada propuesta serıa −3,18125.


2.1.4. Variaciones del metodo: Regula Falsi

El metodo de la Regula-falsi difiere respecto de la biseccion en el calculodel punto cn. Con este metodo se calcula ası:

cn = bn − f(bn)bn − an

f(bn)− f(an).

El paso siguiente es elegir el intervalo formado por los puntos an y cn o bienel formado por cn y bn asegurando que la funcion en los extremos sea de signocontrario.

Ejemplo: Aplicar el metodo de la regula falsi para encontrar una solucionde f(x) = 3x + sen(x)− ex = 0 a partir del intervalo [0, 1].

a0 = 0 b0 = 1 c0 = 0,470990 f(c0) = 0,265160a1 = 0 b1 = 0,470990 c1 = 0,372277 f(c1) = 0,029533a2 = 0 b2 = 0,372277 c2 = 0,361598 f(c2) = 2,94 ∗ 10−3a3 = 0 b3 = 0,361598 c3 = 0,360538 f(c3) = 2,94 ∗ 10−4a4 = 0 b4 = 0,360538 c4 = 0,360433 f(c4) = 2,94 ∗ 10−5

0.2 0.4 0.6 0.8 1

-1

-0.5

0.5

1

regula falsi para f(x) = 3x + sen(x)− ex

20 Capıtulo 2

2.2. Metodo de Newton-Raphson


Consiste en utilizar el polinomio de Taylor de grado 1 como aproximacionf(x).Paso 1. Partimos de un punto inicial x0. Calculamos el polinomio de Taylorde f(x) de grado 1 en x0

P1,x0(x) = f(x0) + f ′(x0)(x− x0).

Paso 2. Utilizamos P1,x0(x) como aproximacion de f(x) y, en vez de resolverf(x) = 0, resolvemos P1,x0 = 0, es decir

f(x0) + f ′(x0)(x− x0) = 0⇒ x = x0 −f(x0)

f ′(x0).

Paso 3. Definimos

x1 = x0 −f(x0)

f ′(x0)

y repetimos los pasos 1, 2 y 3 sustituyendo el punto xn−1 por el punto xn, esdecir construimos la sucesion

xn = xn−1 −f(xn−1)

f ′(xn−1).

Si la funcion f(x) ∈ C1 y la sucesion x0, x1, x2, . . . es convergente, su lımite res la solucion puesto que

r = lımn→

xn = lımn→

[xn−1 −

f(xn−1)

f ′(xn−1)

]= r − f(r)

f ′(r)⇒ f(r)

f ′(r)= 0⇒ f(r) = 0.

Ejemplo Aplica el metodo de Newton-Raphson para encontrar una apro-ximacion de

√2.

Puesto que queremos una aproximacion de√

2, bastara encontrar una apro-ximacion de x2 − 2 = 0. Utilizamos el metodo de Newton-Raphson aplicado ala funcion f(x) = x2 − 2, es decir, construimos la sucesion

xn = xn−1 −f(xn−1)

f ′(xn−1)= xn−1 −

x2n−1 − 2

2xn−1.


Por el teorema de Bolzano, en el intervalo [1, 2] hay una solucion. Comenzamospor el punto x0 = 1,5 :

x0 = 1,5, x1 = 1,4167, x2 = 1,4142, x3 = 1,4142, . . .

La aproximacion serıa 1,4142.


Teorema: Si en un cierto intervalo se verifica que f ∈ C2(R), es creciente,convexa (f ′′(x) ≥ 0) y tiene un cero, entonces el cero es unico y la iteracionde Newton-Raphson convergera a partir de cualquier punto inicial.

-1 1 2 3 4

-2

2

4

6

8

Newton Raphson converge

-2 -1 1 2

0.1

0.2

0.3

0.4

0.5

Newton Raphson no converge

Demostracion. Sea x0, x1, . . . , xn la sucesion construida por el metodo deNewton-Raphson, r la solucion de la ecuacion f(x) = 0 y en = xn − r el erroren el paso n. Por el teorema de Taylor sabemos que:

0 = f(r) = f(xn) + f ′(xn)(−en) +f ′′(θn)

2(−en)2.

Por tanto

enf′(xn)− f(xn) =

f ′′(θn)

2e2n.

De modo que

22 Capıtulo 2

en+1 = xn+1 − r = xn −f(xn)

f ′(xn)− r = en −

f(xn)

f ′(xn)=

=enf

′(xn)− f(xn)

f ′(xn)=

f ′′(θn)e2n2f ′(xn)

≥ 0 (2.1)

Por lo tanto xn+1 ≥ r y por ser f(x) creciente, f(xn) ≥ f(r) = 0. Con estoultimo, se deduce que en es una sucesion decreciente:

en+1 = en −f(xn)

f ′(xn)≤ en.

Puesto que xn+1 ≥ r y el error va decreciendo, la sucesion (xn)n∈N tambien esdecreciente y acotada inferiormente por r, por lo tanto tiene lımite. Ası mismo,la sucesion (en)n∈N es decreciente y acotada inferiormente por 0, por lo tantotambien tiene lımite. Se deduce que

lımn→∞

en+1 = lımn→∞

(en −

f(xn)

f ′(xn)

)=

lımn→∞

en −f(lımn→∞ xn)

f ′(lımn→∞ xn)⇒ f( lım

n→∞xn) = 0

y puesto que f(x) es estrictamente creciente se obtiene que lımn→∞ xn = r(en caso contrario f(lımn→∞ xn) > f(r) = 0), es decir, la sucesion de Newton-Raphson converge a la solucion.


La aproximacion a la solucion en el paso n del metodo de Newton-Raphsones xn. Del apartado anterior (ver (2.1)) sabemos que existe un θn entre xn yr tal que

en+1 =f ′′(θn)

2f ′(xn)e2n.

Si encontramos una constante C tal que

C ≥ maxx∈I |f ′′(x)|2 mınx∈I |f ′(x)|


siendo I un intervalo que contenga a la sucesion (xn)n∈N y a la solucion r,tendremos que

|en+1| ≤ Ce2n

lo que significa que la convergencia de la sucesion de errores es cuadratica. Enestas condiciones

C|en+1| ≤ C2e2n ≤ C2(Ce2n−1)2 ≤ · · · ≤ (Ce0)

2n+1

es decir

C|en| ≤ (Ce0)2n ⇒ |en| ≤

1

C(Ce0)

2n

que es una cota del error que se puede obtener previamente a aplicar el metodo.Si |Ce0| < 1, o lo que es lo mismo |e0| < 1

C , siendo e0 = x0− r, tendremos que

|en| ≤1

C(Ce0)

2n → 0

con lo que tenemos otro criterio de convergencia del metodo.

Ejemplo: Calcula el numero de pasos necesarios para encontrar la raızcuadrada de 3 con un error menor que 10−6

1. Con el metodo de Newton-Raphson

2. Con el metodo de la biseccion.

Puesto que la raız cuadrada de 3 es solucion de x2 − 3 = 0, aplicaremoslos metodos a la funcion f(x) = x2 − 3. Puesto que f(x) es continua en [1, 2]y tiene distinto signo en los extremos del intervalo, por el teorema de Bolzanose sabe que en [1, 2] hay una solucion de f(x) = 0.

Apartado 1). Aplicamos el primer criterio para saber si el metodo deNewton-Raphson convergera en el intervalo [1, 2]. f(x) es un polinomio y portanto es de clase C2. Como f ′(x) = 2x > 0 en [1, 2], la funcion es estrictamentecreciente en [1, 2]. Ademas f ′′(x) = 2 > 0 en [1, 2]. Por todo ello, podemosasegurar que el metodo de Newton-Raphson convergera a partir de cualquierpunto de [1, 2].

Veamos el numero de iteraciones necesarias para que el error sea menorque 10−6. Sabemos que

|en| ≤1

C(Ce0)

2n

24 Capıtulo 2

siendo

C ≥maxx∈[1,2]|f ′′ (x)|

2 mınx∈[1,2] |f ′(x)|=

2

2 mınx∈[1,2] |2x|=

2

2 · 2=

1

2.

Tomando C = 12 y dado que e0 ≤ 1

|en| ≤1

1/2

(1

2e0

)2n

= 2

(1

2

)2n

Si n = 4, entonces |e4| ≤ 3,051 ∗ 10−5

Si n = 5, entonces |e5| ≤ 4,66 ∗ 10−10. Luego basta con 5 iteraciones.Aplicamos el metodo empezando en el punto x0 = 1

xn = xn−1 −f(xn−1)

f ′(xn−1)= xn−1 −

x2n−1 − 3

2xn−1=

x2n−1 + 3

2xn−1.

x0 = 1 x1 = 2 x2 = 1,75 x3 = 1,73414 x4 = 1,73205 x5 = 1,73205

Apartado 2). Sabemos que

|en| ≤1

2n+1(b− a) =

1

2n+1< 10−6

y si n = 19, entonces |e19| ≤ 9,54 ∗ 10−7, es decir, necesitamos 19 iteraciones.

Ejemplo: Dada la ecuacion e−x − x = 0,

1. Aplica el metodo de Newton-Raphson con 4 iteraciones para encontrarla solucion.

2. Calcula una cota del error cometido.

Sea f(x) = e−x − x. Aplicando el teorema de Bolzano al intervalo [0, 1] sepuede asegurar que en dicho intervalo hay una solucion de la ecuacion f(x) = 0.Como f ′(x) = e−x−1 < 0 en [0, 1] no podemos aplicar el criterio para asegurarque el metodo funcione empezando en cualquier punto del intervalo [0, 1]. Noobstante, lo aplicamos empezando en el punto x0 = 0

xn+1 = xn −e−xn − xn−e−xn − 1

x0 = 0 x1 = 0,5 x2 = 0,566311 x3 = 0,567143 x4 = 0,567143


Calculemos una cota del error con la formula

C|en| ≤ (Ce0)2n

siendo

C ≥ maxx∈I |f ′′(x)|2 mınx∈I |f ′(x)|

=maxx∈[0,1] |e−x|

2 mınx∈[0,1] | − e−x − 1|=

1

2(1/e + 1)=

e

2(e + 1)≈ 0,365529.

Por tanto|e4| ≤ C2n−1e2

n

0 ≤ 0,3715116 ≈ 0,33 ∗ 10−7

2.2.4. Variaciones del metodo: Metodo de la secante

Teniendo en cuenta quef ′(xn−1) = lımh→0f(xn−1+h)−f(xn−1)

h este metodose diferencia del de Newton-Raphson en que se sustituye f ′(xn−1) por el valor

aproximado f(xn−2)−f(xn−1)xn−2−xn−1

de modo que la sucesion que queda es

xn = xn−1 − f(xn−1)xn−2 − xn−1

f(xn−2)− f(xn−1).

Por tanto hay que utilizar dos puntos iniciales en lugar de uno. Es un metodomas rapido que el de la biseccion pero mas lento que el de Newton-Raphson.

Ejemplo: Calcula una aproximacion de√

2, con 4 iteraciones del metodode la secante.

Aplicamos el metodo a la solucion de la ecuacion f(x) = x2 − 2 = 0 en elintervalo [1, 2].

xn = xn−1 − (x2n−1 − 2)xn−2 − xn−1

x2n−2 − x2n−1

y tomando tomando x0 = 1 y x1 = 2 nos queda

x0 = 1 x1 = 2 x2 = 1,3333 x3 = 1,4 x4 = 1,4146

26 Capıtulo 2

0.5 1.5 2 2.5

-1

1

2

3

4

Metodo de la secante para f(x) = x2 − 2

2.3. Metodo iterativo de punto fijo


Ciertos problemas de aproximacion se pueden resolver encontrando unasolucion de F (x) = x. Una solucion r de F (x) = x se denomina un puntofijo de la funcion F (x), debido a que F (r) = r. A menudo, para resolver laecuacion F (x) = x se construye la sucesion

xn = F (xn−1).

Cuando la sucesion xn es convergente y F (x) es continua, si r = lımn→∞

xn,

entoncesF (r) = F ( lım

n→∞xn) = lım

n→∞F (xn) = lım

n→∞xn+1 = r,

es decir, el lımite de xn es un punto fijo de F (x). La busqueda de solucionesde F (x) = x empleando sucesiones de la forma xn = F (xn−1), se denominametodo iterativo de punto fijo.

Ejemplo: Encuentra una solucion real de x3 − x − 1 = 0 usando cincopasos de un metodo iterativo de punto fijo.

Como pretendemos resolver la ecuacion f(x) = x3 − x − 1 = 0, por elteorema de Bolzano sabemos que hay una solucion en el intervalo [1, 2].

Resolver f(x) = x3 − x − 1 = 0 es equivalente a resolver una de estasecuaciones

x3 − 1 = x o bien 3√

1 + x = x,


con lo que podemos aplicar el metodo iterativo de punto fijo para resolver laecuacion F (x) = x siendo F (x) = x3−1 o bien F (x) = 3

√1 + x. Construyamos

la sucesion xn = F (xn−1) para ambas funciones, comenzando por un puntodel intervalo [1, 2], por ejemplo por x0 = 1 resulta lo siguiente:

xn = F (xn−1) x0 x1 x2 x3 x4 x5F (x) = x3 − 1 1 0 −1 −2 −9 −730F (x) = 3

√1 + x 1 1,2600 1,3123 1,3224 1,3243 1,3246

Es claro que para la funcion F (x) = x3 − 1, la sucesion no resulta conver-gente. Sin embargo, para la funcion F (x) = 3

√1 + x la sucesion que sale sı es

convergente a un numero r. Dicho numero, puesto que F (x) es continua, veri-ficara que F (r) = r, es decir, f(r) = 0. La mejor aproximacion de r obtenidaserıa 1,3246.


Definicion: Sea F : R → R una funcion. Se dice que F es contractiva enun conjunto C ⊂ R si existe 0 ≤ λ < 1 tal que

|F (x)− F (y)| ≤ λ|x− y| ∀x, y ∈ C.

Ejemplo: Comprueba que la funcion F (x) = 3√

1 + x es contractiva en[0, 1].

Por el teorema del valor medio

|F (y)− F (x)| = |F ′(θ)(y − x)| =

∣∣∣∣∣ 1

3 3√

(1 + θ)2(y − x)

∣∣∣∣∣ ≤≤∣∣∣∣ 1

3 3√

1(y − x)

∣∣∣∣ =1

3|y − x|

para todo x, y ∈ [0, 1]. Por tanto la funcion es contractiva con λ = 1/3.

Teorema: Si F es una funcion contractiva que va de un cerrado C ⊂ Ren C (F (C) ⊂ C), entonces F tiene un unico punto fijo r ∈ C. Ademasr es el lımite de cualquier sucesion que se obtenga a partir de la expresionxn+1 = F (xn) siendo x0 ∈ C.

28 Capıtulo 2

Demostracion. Sea x0, x1, x2, . . . , xn, . . . la sucesion generada mediante laexpresion xn+1 = F (xn). Veamos que (xn)n∈N converge. Por ser contractivase tiene que

|xn − xn−1| = |F (xn−1)− F (xn−2)| ≤ λ|xn−1 − xn−2| ≤λ2|xn−2 − xn−3| ≤ · · · ≤ λn−1|x1 − x0|.

Por otra parte

xn = x0 + x1 − x0 + x2 − x1 + · · ·+ xn − xn−1 = x0 +

n∑k=1

(xk − xk−1)

por lo que tomando lımites

lımn→∞

xn = x0 +

∞∑k=1

(xk − xk−1)

La sucesion (xn)n∈N sera convergente si la serie∑∞

k=1(xk − xk−1) es conver-gente. Esta ultima lo es porque∣∣∣∣∣

∞∑k=1

(xk − xk−1)

∣∣∣∣∣ ≤∞∑k=1

|xk − xk−1| ≤∞∑k=1

λk−1(x1 − x0) = |x1 − x0|1

1− λ

Sea r = lımn→∞ xn. Por la construccion de la sucesion, se tiene que F (r) = r,es decir, hay un punto fijo. r ∈ C por ser C un conjunto cerrado y xn ∈ C. Faltaver que es unico. Supongamos que hubiese dos puntos fijos r1 y r2. Entonces

|r1 − r2| = |F (r1)− F (r2)| ≤ λ|r1 − r2| < |r1 − r2|

llegando a contradiccion.


Si r es solucion de f(x) = 0 obtenida por iteracion de una funcion F (x),entonces si f(x) es continua y derivable, existe θn ∈ R entre r y xn tal que

f(r)− f(xn)

r − xn= f ′(θn).


Suponiendo que xn ≤ r y puesto que f(r) = 0 y −en = r − xn, se deduce que

f(xn)

en= f ′(θn)⇒ en =

f(xn)

f ′(θn)⇒ |en| ≤

|f(xn)|mınθ∈[xn,r] |f ′(θ)|

.

La formula anterior nos proporciona un metodo para encontrar una cota deerror a posteriori y se puede aplicar a cualquier metodo iterativo de punto fijo.

Si F (x) es contractiva de constante λ entonces

|xn+p − xn| = |xn+p − xn+p−1|+ |xn+p−1 − xn+p−2|+ · · ·+ |xn+1 − xn| ≤≤ λn+p−1|x1 − x0|+ λn+p−2|x1 − x0|+ · · ·+ λn|x1 − x0| =

= (λn+p−1 + λn+p−2 + · · ·+ λn)|x1 − x0| == λn(λp−1 + λp−2 + · · ·+ 1)|x1 − x0|

de modo que cuando p tiende a infinito

|en| = |r − xn| ≤ λn∞∑k=0

λk = λn

(1

1− λ

)|x1 − x0|

que es una formula que permite calcular una cota de error a priori.

Ejemplo: Encuentra una solucion de la ecuacion 4+ 13 sin(2x)−x = 0 con

un metodo iterativo de punto fijo estimando una cota del error a priori y otraa posteriori.

Sea f(x) = 4 + 13 sin(2x)− x y F (x) = 4 + 1

3 sin(2x). Aplicando el teoremade Bolzano a la funcion f(x) en el intervalo [−3,5, 4,5] se deduce que hayuna solucion de f(x) = 0 en dicho intervalo. Veamos si F (x) es contractivaaplicando el teorema del valor medio

|F (y)− F (x)| = |F ′(θ)(y − x)| = |23

cos(2x)(y − x)| ≤ 2

3|y − x|.

Se deduce que F (x) es contractiva con constante λ = 23 .

Por otra parte, puesto que −3,5 ≤ 4 + 13 sin(2x) ≤ 4,5, se deduce que

F ([−3,5, 4,5]) ⊂ [−3,5, 4,5].Por tanto, la sucesion xn+1 = F (xn) convergera a un punto fijo a partir de

cualquier valor inicial x0. Calculamos cinco iteraciones empezando en x0 = 0 :

30 Capıtulo 2

x0 x1 x2 x3 x4 x50 4 4,329786082 4,230895147 4,273633799 4,256383406

Una cota del error para x5 a priori serıa

|e5| ≤(

2

3

)5(

1

1− 23

)|4− 0| ≈ 1,5802

que no es una cota muy buena.

Una cota del error para x5 a posteriori serıa

|e5| ≤|f(x5)|

mınθ∈[−3,5,4,5] |23 cos(2θ)− 1|≤ 0,00719542

1/3= 0,00239847.

2.4. Raıces de polinomios

No existe ninguna formula algebraica para la resolucion de raıces de poli-nomios de grado mayor que 4, por tanto, hay que aplicar metodos numericospara el calculo de dichas raıces.

Definicion: Una solucion r de una ecuacion f(x) = 0 se dice que tienemultiplicidad n si

f(r) = f ′(r) = · · · = fn−1)(r) = 0 fn)(r) = 0.

Ejemplo: r = 1 es una solucion de multiplicidad 2 para la ecuacion

x2 − 2x + 1 = 0.

Para la localizacion de raıces de polinomios con multiplicidad par (como

en el ejemplo anterior) no podemos emplear el teorema de Bolzano. Sin em-bargo podemos intentar calcular las raıces de otro polinomio con las mismassoluciones que el nuestro pero todas con multiplicidad 1.


0.5 1 1.5 2

0.2

0.4

0.6

0.8

1

f(x) = x2 − 2x + 1

Teorema [Teorema fundamental del algebra]: Un polinomio de grado ncon coeficientes reales tiene exactamente n ceros o raıces entre los numeroscomplejos C, contando cada cero tantas veces como indique su multiplicidad.

De modo que dado un polinomio Pn(x), podemos descomponerlo ası:

Pn(x) = an(x− rk)mk(x− rk−1)mk−1 . . . (x− r1)

m1

siendo mi la multiplicidad de la raız i-esima. Si derivamos, obtenemos que

P ′n(x) = an(x− rk)mk−1(x− rk−1)

mk−1−1 . . . (x− r1)m1−1M(x)

siendo M(x) un polinomio tal que M(ri) = 0. El maximo comun divisor dePn(x) y de P ′

n(x) es

D(x) = M.C.D(Pn(x), P ′n(x)) =

an(x− rk)mk−1(x− rk−1)mk−1−1 . . . (x− r1)

m1−1

por lo que el polinomio

Q(x) =P (x)

M.C.D(Pn(x), P ′n(x))

= an(x− rk)(x− rk−1) . . . (x− r1)

tiene las mismas raıces que P (x) pero todas con multiplicidad 1. Bastarıaaplicar los metodos a Q(x) para localizar las raıces de Pn(x). El problema esel calculo de Q(x) y de M.C.D(Pn(x), P ′

n(x)).

32 Capıtulo 2

2.4.1. Separacion de raıces. Sucesion de Sturm.

La sucesion de Sturm nos permite calcular el numero de raıces de un po-linomio en un intervalo [a, b]. Sea Pn(x) un polinomio de grado n y P ′

n(x) suderivada. Para construir la sucesion de Sturm aplicamos la division de poli-nomios repetidas veces a Pn(x), P ′

n(x) y a los restos obtenidos cambiados designo:

Pn(x) = c1(x)P ′n(x) + r1(x)

P ′n(x) = c2(x)(−r1(x)) + r2(x)

−r1(x) = c3(x)(−r2(x)) + r3(x). . .

−rk−2(x) = ck(x)(−rk−1(x)) + rk(x)

Puede ocurrir lo siguiente:1. Hay un resto rk(x) = 0. En este caso rk−1(x) es el maximo comun divisor

de Pn(x) y de P ′n(x). Construimos entonces la sucesion de Sturm correspon-

diente al polinomio

Q(x) =Pn(x)

rk−1(x)

que es un polinomio que tiene las mismas raıces que Pn(x) pero todas simples.2. Llegamos a un resto rk(x) constante distinto de 0. En este caso la suce-

sion de Sturm es

Pn(x), Pn−1(x),−r1(x), . . . ,−rk(x).

La diferencia entre el numero de cambios de signos de las sucesiones siguientes

Pn(a), Pn−1(a),−r1(a), . . . ,−rk(a) Pn(b), Pn−1(b),−r1(b), . . . ,−rk(b)

es el numero de raıces de Pn(x) en el intervalo [a, b].


Ejemplo: Localiza en intervalos disjuntos las raıces del P (x) = 36x4 −12x3 − 11x2 + 2x + 1.

La derivada de P (x) es P ′(x) = 144x3−36x2−22x+2. Podemos multiplicarlos polinomios implicados por numeros positivos para facilitar los calculos sinque ello afecte al objetivo perseguido. Comenzamos los calculos:

1. Multiplicamos Pn(x) por 4 y efectuamos la primera division:

144x4−48x3−44x+8x+4 = (x−1)(144x3−36x2−22x+2)+(−300x2+50x+50).

2. El tercer miembro de la sucesion de Sturm es 300x2 − 50x − 50. Efec-tuamos la siguiente division:

144x3 − 36x2 − 22x + 2 = 2/50(12x− 1)(300x2 − 50x− 50) + 0

3. Como hay un resto que es 0 el maximo comun divisor de P (x) y P ′(x) es300x2−50x−50 o bien dividiendo por 50 6x2−x−1. calculamos el polinomio

Q(x) =P (x)

6x2 − x− 1= 6x2 − x− 1.

4. Construimos la sucesion de Sturm para Q(x).

6x2 − x− 1 = (− 1

24+

x

2)(12x− 1) +

−25

24

Con la sucesion es

6x2 − x− 1, 12x− 1,25

24

Estudiamos algunos cambios de signo:

−∞ -2 0 2 ∞6x2 − x− 1 + + - + +

12x− 1 - - - + +2524 + + + + +

cambios de signo 2 2 1 0 0

En el intervalo (−∞,−2) y en (2,∞) no hay ninguna raız porque en los ex-tremos de los intervalos la sucesion de Sturm tiene los mismos cambios designo.

En el intervalo (−2, 0) hay una raız y hay otra en el intervalo (0, 2) porquela diferencia de cambios de signo en los extremos en ambos casos es 1.

34 Capıtulo 2

Ejemplo: Localiza en intervalos disjuntos las raıces del P (x) = x4 +2x3−3x2 − 4x− 1.

La derivada es P ′(x) = 4x3 + 6x2 − 6x− 4 = 2(2x3 + 3x2 − 3x− 2). Parafacilitar el calculo multiplicamos P (x) por 2, P ′(x) por 1

2 y calculamos es tercermiembro de la sucesion de Sturm:

2x4 + 4x3 − 6x2 − 8x− 2 = (x +1

2)(2x3 + 3x2 − 3x− 2) + (−9

2x2 − 9

2x− 1).

Podemos multiplicar por 2 el resto con lo que −r1(x) = 9x2+9x+2. Para faci-litar los calculos multiplicamos P ′(x) por 9

2 . Calculamos el siguiente miembrode la sucesion:

18x3 + 27x2 − 27x− 18 = (2x + 1)(9x2 + 9x + 2) + (−40x− 20).

Dividiendo por 20 podemos suponer que −r2(x) = 2x + 1. Multiplicamos−r1(x) por 2 y calculamos −r3(x):

18x2 + 18x + 4 = (9x +9

2)(2x + 1) +

−1

2.

Con lo que la sucesion queda ası:

x4 + 2x3 − 3x2 − 4x− 1, 2x3 + 3x2 − 3x− 2, 9x2 + 9x + 2, 2x + 1,1

2.

Separamos las raıces en intervalos disjuntos:

−∞ -3 -2 -1 -1/2 0 1 2 ∞x4 + 2x3 − 3x2 − 4x− 1 + + - - - - - + +

2x3 + 3x2 − 3x− 2 - - 0 + - - 0 + +9x2 + 9x + 2 + + + + + + + + +

2x + 1 - - - - 0 + + + +12 + + + + + + + + +

cambios de signo 4 4 3 3 2 1 1 0 0

En el intervalo (−∞,−3), (−2,−1), (0, 1) y en (2,∞) no hay ninguna raızporque en los extremos de los intervalos la sucesion de Sturm tiene los mismoscambios de signo.

En cada uno de los intervalos (−3,−2), (−1,−1/2), (−1/2, 0) y (1, 2) hayuna raız porque la diferencia de cambios de signo en los extremos en amboscasos es 1.


2.4.2. Acotacion de raıces

Proposicion: Sea P (x) = anxn + an−1x

n−1 + · · ·+ a0. Si r es una raız deP (x), entonces

|r| < 1 +A

|an|siendo A = max0≤i≤n−1 |ai|.

Demostracion. Si A = 0, entonces el polinomio solo tendrıa la raız 0 y secumplirıa el resultado. Supongamos que A = 0 y que |r| > 1.

0 = |P (r)| = |anrn + an−1rn−1 + · · ·+ a0| ≥

≥ |anrn| − |an−1rn−1 + · · ·+ a0|

≥ |anrn| −[|an−1||rn−1|+ · · ·+ |a0|

]≥

≥ |anrn| −A[|r|n−1 + · · ·+ |r|+ 1

]=

|anrn| −A|r|n − 1

|r| − 1> |an||r|n −A

|r|n

|r| − 1

> |r|n(|an| − A

|r| − 1

)es decir,

0 > |an| − A

|r| − 1⇒ |an| < A

|r| − 1⇒ |r| − 1 <

A

|an|con lo que concluimos que

|r| < 1 +A

|an|.

Si |r| ≤ 1, es trivial que |r| < 1 + A|an| .

Ejemplo: Acota las raıces del polinomio P (x) = x4 + 2x3 − 3x2 − 4x− 1.

Com A = max1≤i≤n−1 |ai| = 4, usando la proposicion anterior obtenemosque cualquier raız r del polinomio verifica lo siguiente

|r| < 1 +A

|an|= 1 +

4

1= 5⇒ −5 < r < 5.

36 Capıtulo 2

2.4.3. Raıces de polinomios con el algoritmo de Horner

Si usamos el metodo de Newton-Raphson para encontrar las raıces de unpolinomio P (x) hemos de calcular P (xn) y P ′(xn). El algoritmo de Hornerhace estos ultimos calculos de forma sencilla para polinomios.

Sea x0 ∈ R y P (x) = anxn + an−1x

n−1 + · · ·+ a1x + a0. Sea

P (x) = Q(x)(x− x0) + R

(R sera constante) el resultado de dividir P (x) entre (x − x0). Supongamosque

Q(x) = bn−1xn−1 + bn−2x

n−2 + · · ·+ b1x + b0.

Entonces P (x0) = R y P (x) = Q(x)(x− x0) + P (x0), con lo que

P (x)− P (x0) = Q(x)(x− x0)

es decir

anxn + an−1x

n−1 + · · ·+ a1x + a0 − P (x0) =

= (bn−1xn−1 + bn−2x

n−2 + · · ·+ b1x + b0)(x− x0)

Desarrollando e igualando los coeficientes de ambos terminos se deduceque

bn−1 = an, bn−2 = an−1+bn−1x0, . . . , b0 = a1+b1x0 P (x0) = a0+b0x0.

La obtencion de los coeficientes bi se pueden representar graficamente ası

an an−1 an−2 . . . a1 a0x0 bn−1x0 bn−2x0 . . . b1x0 b0x0

bn−1 bn−2 bn−3 . . . b0 P (x0)

que tambien es una forma sencilla de calcular P (x0). Ademas,

P (x) = Q(x)(x−x0) +R⇒ P ′(x) = Q′(x)(x−x0) +Q(x)⇒ P ′(x0) = Q(x0).

De modo que para calcular P ′(x0) podemos aplicar el algoritmo grafico anteriora Q(x) :

an an−1 an−2 . . . a1 a0x0 bn−1x0 bn−2x0 . . . b1x0 b0x0

bn−1 bn−2 bn−3 . . . b0 P (x0)x0 cn−2x0 cn−3x0 . . . c0x0

cn−2 cn−3 cn−4 . . . P ′(x0)


Ejemplo: Siendo P (x) = x4 + 2x3− 3x2− 4x− 1, calcula con el algoritmode Horner P (2) y P ′(2).

2 1 −3 4 −52 4 10 14 36

2 5 7 18 312 4 18 50

2 9 25 68

Se concluye que P (2) = 31 y P ′(2) = 68.

2.4.4. Raıces multiples

Si una ecuacion f(x) = 0 tiene soluciones multiples, el calculo puede serproblematico. Si f(x) es un polinomio, podemos intentar calcular las raıces de

Q(x) =P (x)

MCD(P (x), P ′(x))

que sera un polinomio con las mismas raıces que f(x) pero todas con multi-plicidad 1.

Sea o no f(x) un polinomio, se puede intentar aplicar el metodo de Newton-Raphson. Si la sucesion generada converge muy lentamente, la causa puede en-contrarse en esa multiplicidad de las raıces. Para detectar de que multiplicidades la raız en cuestion, para cada termino de la sucesion xn generada calculamosf ′(xn), f ′′(xn), f ′′′(xn), . . . , fk)(xn). Si las derivadas hasta orden k−1 son casi0 en los valores xn y fk)(xn) no es cercano a 0, entonces puede ocurrir que lasolucion r verifique que

fn)(r) = 0 ∀ 0 ≤ n < k y fk)(r) = 0,

es decir que r sea una solucion de f(x) = 0 de multiplicidad k. En este caso,se puede acelerar la convergencia de la sucesion calculandola ası

xn+1 = xn − kf(xn)

f ′(xn).

38 Capıtulo 2

El error se obtendrıa ası

|en| <|f(xn)|

mınθ∈[a,b] |f ′(θ)|

Ejercicio: Sea f(x) = (x − 1)3. Aplica Newton-Raphson para encontraruna solucion de f(x) = 0 empezando en x0 = 0.

Capıtulo 3

Sistemas Lineales

El objetivo de este tema es el desarrollo de tecnicas de busqueda de solu-ciones de sistemas de ecuaciones lineales como el siguiente:

a11x1 + a12x2 + · · ·+ a1nxn = b1a21x1 + a22x2 + · · ·+ a2nxn = b2

· · ·an1x1 + an2x2 + · · ·+ annxn = bn

Denotando A = (aij)1≤i,j≤n, b = (bi)1≤i≤n y x = (xi)1≤i≤n, el sistema anteriorse puede expresar

A · x = b

3.1. Algebra de matrices

Sea A una matriz.

A es de orden m×n si tiene tiene m filas y n columnas (A = (aij)1≤i≤m,1≤j≤n ∈Mm×n). Por ejemplo (

3 5 22 1 3

)tiene orden 2× 3.

La matriz traspuesta de A es la matriz At que resulta de intercambiar lasfilas de A por las columnas de A.

39

40 Capıtulo 3

A es simetrica si At = A.Si α ∈ R y B es una matriz, podemos realizar las siguientes operaciones:- Multiplicar una matriz por un numero:

α ·A = (αaij)1≤i≤m,1≤j≤n.

- Sumar dos matrices A y B si A,B ∈Mm×n:

A + B = (ai,j + bij)1≤i,j≤m.

- Multiplicar dos matrices A ∈Mm×p y B ∈Mp×n:

A ·B = (

p∑k=1

aikbkj)1≤i≤m,1≤j≤n ∈Mm×n.

La matriz identidad de orden n× n es

Id =

1 0 0 . . . 00 1 0 . . . 0

. . .0 0 . . . 0 1

Se verifica que A · Id = Id ·A = A, ∀A ∈Mn×n.

B es una matriz inversa de A por la derecha si B ·A = IdB es una matriz inversa de A por la izquierda si A ·B = Id.Si A ∈ Mn×n, entonces la inversa por la izquierda es tambien inversa por

la derecha. En este caso la inversa es unica y se denota A−1

A−1 ·A = A ·A−1 = Id.

Si A ∈Mn×n tiene inversa, entonces A es regular, invertible o no singular.La matriz inversa es

A−1 =1

|A|(Aij)

t,

siendo Aij el adjunto del elemento aij .En este caso el sistema A · x = b tiene como solucion x = A−1 · b.Se llaman operaciones elementales de una matriz al intercambio de dos

filas o columnas, al producto de una fila o columna por un numero distintode 0, a la suma de una fila o columna otra fila o columna multiplicada porun numero, o cualquier combinacion finita de las operaciones anteriores. Cada

Sistemas lineales 41

operacion fundamental en una matriz A se puede expresar mediante productode A por otra matriz, que se llama matriz fundamental. Por ejemplo, dada lamatriz

A =

a11 a12 a13a21 a22 a23a31 a32 a33

,

el intercambio de la segunda y tercera filas se puede expresar ası:1 0 00 0 10 1 0

·a11 a12 a13a21 a22 a23a31 a32 a33

=

a11 a12 a13a31 a32 a33a21 a22 a23

.

El producto de la segunda fila por un numero α se puede expresar ası:1 0 00 α 00 0 1

·a11 a12 a13a21 a22 a23a31 a32 a33

=

a11 a12 a13αa21 αa22 αa23a31 a32 a33

Y la suma a la tercera fila de la segunda fila multiplicada por un numero α es:1 0 0

0 1 00 α 1

·a11 a12 a13a21 a22 a23a31 a32 a33

=

a11 a12 a13a21 a22 a23

αa21 + a31 αa22 + a32 αa23 + a33

El producto de matrices1 0 0

0 0 α0 1 αλ

·a11 a12 a13a21 a22 a23a31 a32 a33

representa aplicar a la segunda matriz los siguientes cambios:

- Sustituir la fila segunda por la tercera multiplicada por α.- Multiplicar la tercera fila por αλ y sumarle la segunda.Si A ∈ Mn×n tiene inversa, entonces existe un conjunto de matrices fun-

damentales Ek con k = 1, . . . , p, tales que

EpEp−1 . . . E1A = Id,

con lo que multiplicando por A−1 por la derecha en ambos terminos se obtieneque

EpEp−1 . . . E1Id = A−1,

es decir, aplicando las mismas operaciones fundamentales a la matriz Id ob-tendremos la matriz inversa A−1.

42 Capıtulo 3

Ejemplo: Calcula la inversa de la matriz

A =

1 2 31 3 32 4 7

.

Partimos de A y la matriz Id. El primer cambio consiste en restar a la segundafila la primera (correspondiente al producto por una cierta matriz fundamentalE1):

E1 ·A =

1 2 30 1 02 4 7

E1 · Id =

1 0 0−1 1 00 0 1

El segundo cambio (multiplicar por E2) consiste en restar a la tercera dosveces la primera:

E2 · E1 ·A =

1 2 30 1 00 0 1

E2 · E1 · Id =

1 0 0−1 1 0−2 0 1

A la primera le resto dos veces la segunda (E3):

E3 · E2 · E1 ·A =

1 0 30 1 00 0 1

E3 · E2 · E1 · Id =

3 −2 0−1 1 0−2 0 1

A la primera le resto tres veces la tercera (E4):

E4·E3·E2·E1·A =

1 0 00 1 00 0 1

E4·E3·E2·E1·Id =

9 −2 −3−1 1 0−2 0 1

= A−1.

Esta ultima matriz es la inversa de A.


Para toda matriz A ∈Mn×n son equivalentes:

1. A−1 existe.

2. |A| = 0.

3. Los vectores fila de la matriz A son linealmente independientes.

4. Los vectores columna de la matriz A son linealmente independientes.

5. Para cada b ∈ Rn, el sistema A · x = b tiene una unica solucion.

6. A es producto de matrices elementales.

3.1.1. Valores propios y vectores propios

Sea A ∈Mn×n. Si existe un vector x ∈ Rn no nulo y un numero λ ∈ R talque A · x = λx, entonces, λ es un valor propio o autovalor de la matriz A, yx es un vector propio o autovector para el autovalor λ. Los valores propios deuna matriz A son las raıces del llamado polinomio caracterıstico |A − λ · Id|.Los vectores propios de A para el autovalor λ, son los elementos de nucleo dela aplicacion lineal A− λ · Id.

3.1.2. Matriz definida

Sea A ∈ Mn×n. Se dice que A es definida positiva (equivalentemente,negativa) si para todo x ∈ Rn no nulo, se verifica que

xt ·A · x > 0.

Ejemplo. La matriz A =

(2 11 2

)es definida positiva.

Sea (x1, x2) un vector no nulo.

(x1, x2) ·(

2 11 2

)·(x1x2

)= (x1, x2) ·

(2x1 + x2x1 + 2x2

)=

2x21 + x1x2 + x1x2 + 2x22 = (x1 + x2)2 + x21 + x22 > 0.

Si A es una matriz definida positiva y simetrica, entonces los autovalores

de A son todos numeros reales y positivos.

44 Capıtulo 3

3.2. Resolucion de sistemas de ecuaciones lineales:metodo de Gauss

Dado el sistema A · x = b, si |A| = 0, puede no haber solucion. Si |A| = 0,existe una unica solucion y se puede calcular ası:

x1 =

∣∣∣∣∣∣∣∣b1 a12 . . . a1nb2 a22 . . . a2n

. . .bn an2 . . . ann

∣∣∣∣∣∣∣∣|A|

, · · · , xn =

∣∣∣∣∣∣∣∣a11 . . . a1n−1 b1a21 a2n−1 b2

. . .a2n . . . ann−1 bn

∣∣∣∣∣∣∣∣|A|

Los calculos anteriores son poco practicos cuando el numero de ecuacioneses grande. Por ejemplo, para n = 50 se necesitan 1064 operaciones. Para evi-tar este problema se utilizan metodos numericos de resolucion de ecuacioneslineales que pueden ser:

• Directos: proporcionan un resultado que sera exacto salvo errores deredondeo tras un numero determinado de operaciones.

• Iterativos: construyen una sucesion de soluciones aproximadas, de modoque en cada paso se mejora la aproximacion anterior.

Se dice que dos sistemas de ecuaciones son equivalentes si tienen las mismassoluciones. En un sistema de ecuaciones se pueden efectuar operaciones ele-mentales en la matriz de los coeficientes y en la de los terminos independientesy resulta un sistema equivalente. Es decir, resulta un sistema equivalente siefectuamos alguna de las siguientes operaciones:

1. Multiplicamos una ecuacion por un numero distinto de cero.

2. Sumamos a una ecuacion otra ecuacion multiplicada por un numero.

3. Intercambiamos dos ecuaciones.


3.2.1. Metodo de Gauss

El metodo de Gauss clasico es un metodo directo que transforma medianteoperaciones elementales el sistema de ecuaciones que se pretende resolver enun sistema equivalente triangular de facil resolucion.

El metodo de Gauss puede resumirse en lo siguiente:

Paso 1. Si a11 = 0, se elige a11 como elemento pivote. Si a11 = 0, se inter-cambian las filas para que esto no suceda. Sumando a cada una las filasrestantes la fila del elemento pivote multiplicada por ciertos numeros,llamados multiplicadores, se convierten en 0 los elementos que se en-cuentran en la primera columna distintos del pivote. En las operacionesde sumas de filas hay que incluir los terminos independientes.

Paso i (i ≥ 2). Si aii = 0, se elige aii como elemento pivote. Si aii = 0,se intercambian las filas de orden mayor o igual que i para que esto nosuceda. Sumando a cada una de las filas que ocupan el lugar j (j > i) lafila del elemento pivote multiplicada por numeros, llamados multiplica-dores, se convierten en 0 los elementos que se encuentran en la columnai de la fila j (j > i). En las operaciones de sumas de filas hay que incluirlos terminos independientes.

Se efectuan los pasos necesarios hasta que la matriz de los coeficientes seatriangular. Conseguido esto, el calculo de la solucion es sencillo siguiendo elsiguiente orden: xn, xn−1, . . . , x1.

Ejemplo: Resuelve el siguiente sistema de ecuaciones:

6x1 − 2x2 + 2x3 + 4x4 = 1212x1 − 8x2 + 6x3 + 10x4 = 343x1 − 13x2 + 9x3 + 3x4 = 27−6x1 + 4x2 + x3 − 18x4 = −38.

En forma matricial es sistema serıa6 −2 2 412 −8 6 103 −13 9 3−6 4 1 −18

x1x2x3x4

=

123427−38

.

En el sistema anterior el metodo de Gauss podrıa aplicarse ası:

46 Capıtulo 3

1. Paso 1. Elegimos como elemento pivote el termino 6 de la primera filay efectuamos las siguientes operaciones:

(fila 2a)− 2(fila 1a) (fila 3a)− 1/2(fila 1a) (fila 4a)− (−1)(fila 1a)

obtenemos el sistema6 −2 2 40 −4 2 20 −12 8 10 2 3 −14

x1x2x3x4

=

121021−26

.

Paso 2. Elegimos como elemento pivote el termino −4 de la segunda fila yefectuamos las siguientes operaciones:

(fila 3a)− 3(fila 2a) (fila 4a)− (−1/2)(fila 2a)

obtenemos el sistema6 −2 2 40 −4 2 20 0 2 −50 0 4 −13

x1x2x3x4

=

1210−9−21

.

Paso 3. Elegimos como elemento pivote el termino 2 de la tercera fila yefectuamos la siguiente operacion:

(fila 4a)− 2(fila 3a)

obtenemos el sistema6 −2 2 40 −4 2 20 0 2 −50 0 0 −3

x1x2x3x4

=

1210−9−3

.

que es un sistema de ecuaciones triangular superior de facil resolucion:

6x1 − 2x2 + 2x3 + 4x4 = 12−4x2 + 2x3 + 2x4 = 10

2x3 − 5x4 = −9−3x4 = −3

⇒x1 = 1/6(12− 4x4 − 2x3 + 2x2) = 1x2 = −1/4(10− 2x4 − 2x3) = −3x3 = 1/2(−9 + 5x4) = −2x4 = 1


Se verifica que A = L · U , siendo A la matriz de los coeficientes inicial, Ula matriz de los coeficientes final y L la matriz formada por los numeros quehemos utilizado para multiplicar las filas (multiplicadores):

L =

1 0 0 02 1 0 0

1/2 3 1 0−1 −1/2 2 1

3.2.2. Metodo de Gauss con pivoteo

En un sistema como este(0,001 1

1 1

)(x1x2

)=

(12

)el metodo de Gauss no funcionarıa correctamente por errores de redondeo. Engeneral, si ϵ = 0 es pequeno, el sistema(

ϵ 11 1

)(x1x2

)=

(12

)aplicando el metodo de Gauss se transformarıa en(

ϵ 10 1− 1

ϵ

)(x1x2

)=

(1

2− 1ϵ

)con lo que las soluciones, si consideramos el redondeo, serıan:

x2 =2− 1

ϵ

1− 1ϵ

≃ 1 ⇒ x1 = (1− x2)1

ϵ≃ 0.

Si ϵ = 0,001, obtendrıamos que x2 = 0,998998998999, que, redondeando a dosdecimales, serıa x2 = 0, por lo que x1 = 0. Sin embargo las soluciones exactasson x2 = 998

999 ≈ 1 y x1 = 1000/999 ≈ 1. De hecho, al hacer los calculos con leprograma Maxima o Mathematica, con ϵ = 10−17 se obtiene x2 = 1 y x1 = 0.

Esto sucede cuando el elemento pivote en algun paso es muy pequenocomparado con el resto de los elementos a los que divide. En este caso el

48 Capıtulo 3

multiplicador correspondiente sera muy grande y la ecuacion que origina escasi un multiplo de la ecuacion del elemento pivote:

(fila 2a)− 1

ϵ(fila 1a) ≈ −1

ϵ(fila 1a).

Este problema se puede resolver de dos formas:

1. Pivoteo parcial: si estamos en el paso i, se elige como elemento pivote elelemento aki que tenga mayor valor absoluto entre los de la columna i,con k ≥ i, independientemente de si aii es cero o no, y se reordenan lasfilas.

2. Pivoteo total: si estamos en el paso i, se elige como elemento pivote elelemento akj que tenga mayor valor absoluto entre los de la columna iy los de la fila i, con k, j ≥ i, independientemente de si aii es cero ono, y se reordenan las filas y columnas. Esta opcion tiene un importanteinconveniente pues hay que reordenar las incognitas.

Ejemplo: Resuelve el siguiente sistema de ecuaciones:

2x2 + x4 = 02x1 + 2x2 + 3x3 + 2x4 = −2

4x1 − 3x2 + x4 = −76x1 + x2 − 6x3 − 5x4 = 6.

aplicando el pivoteo parcial.

La matriz de los coeficientes y de los terminos independientes serıa0 2 0 1 | 02 2 3 2 | −24 −3 0 1 | −76 1 −6 −5 | 6

El pivoteo parcial serıa este proceso:

Paso 1. Elegimos como elemento pivote el termino 6 de la ultima fila porser el de mayor valor absoluto de entre los elementos de la primera columna.Obtenemos la matriz

6 1 −6 −5 | 62 2 3 2 | −24 −3 0 1 | −70 2 0 1 | 0


Efectuamos las operaciones indicadas obtenemos la matriz siguiente:

(fila 2a)− 26 (fila 1a)

(fila 3a)− 46 (fila 1a)

6 1 −6 −5 | 60 1,6667 5 3,6667 | −40 −3,6667 4 4,3333 | −110 2 0 1 | 0

Paso 2. Elegimos como elemento pivote el termino −3,667 de la tercera fila porser el de mayor valor absoluto de entre los elementos de la segunda columna.Obtenemos la matriz

6 1 −6 −5 | 60 −3,6667 4 4,3333 | −110 1,6667 5 3,6667 | −40 2 0 1 | 0

Efectuamos las operaciones indicadas obtenemos la matriz siguiente:

(fila 3a)− 1,6667−3,6667 (fila 2a)

(fila 4a)− 2−3,6667 (fila 2a)

6 1 −6 −5 | 60 −3,6667 4 4,3333 | −110 0 6,8182 5,6364 | −9,00010 0 2,1818 3,3636 | −5,9999

Paso 3. El elemento pivote sera 6,8182 de la tercera fila y no es necesariointercambiar filas. Efectuamos la operacion indicada y obtenemos la matrizsiguiente:

(fila 4a)− 2,18186,8182 (fila 3a)

U =

6 1 −6 −5 | 60 −3,6667 4 4,3333 | −110 0 6,8182 5,6364 | −9,00010 0 0 1,5600 | −3,1199

que se corresponde con el sistema de ecuaciones

6x1 + x2 +−6x3 − 5x4 = 6−3,6667x2 + 4x3 + 4,3333x4 = −11

6,8182x3 + 5,6364 = −9,00011,5600x4 = −3,1199.

.

Este sistema se resuelve facilmente por ser triangular despejando las incognitasen orden inverso. Se obtiene la solucion:

50 Capıtulo 3

x4 = −1,9999, x3 = 0,33325, x2 = 1,0000, x1 = −0,50000,

que es una respuesta aceptable teniendo en cuenta que la solucion exactaes

x4 = −2, x3 = 1/3, x2 = 1, x1 = −1/2.

La matriz L de los multiplicadores en el ejemplo anterior es la siguiente:

L =

1 0 0 0

0,6667 1 0 00,3333 −0,45454 1 0

0 −0,54545 0,32 1

y se verifica que

L · U =

6 1 −6 54 −3 0 12 2 3 20 2 0 1

:= A′.

El resultado no es exactamente A porque al aplicar el pivoteo hemos permu-tado filas. En todo caso se verifica que existe una matriz P , que representapermutaciones de filas si empleamos pivoteo parcial, tal que

P ·A = A′ = L · U.

Esto ultimo nos permite calcular el determinante de A de forma sencilla alser |L| = 1, |P | = (−1)p (p es el numero de permutaciones de filas) y U unamatriz triangular:

|A| = |P−1 ·A′| = |P−1 · L · U | = |P−1||L||U | = (−1)p|U |.

El metodo de Gauss permite:

1. Encontrar las soluciones de un sistema de ecuaciones.

2. Calcular el determinante de la matriz de los coeficientes.

3. Descomponer la matriz de los coeficientes A en un producto de matricesP · A′ donde P representa una permutacion de filas de A y A′ = L · U ,siendo L y U matrices triangulares inferior y superior, respectivamente.A = A′ si no se utiliza el pivoteo.


3.3. Factorizacion LU. Factorizacion de Cholesky

Dado el sistema de ecuaciones Ax = b, si aplicamos el metodo de Gauss sinpivoteo solo a la matriz de los coeficientes y no a los terminos independientesb, obtendremos una descomposicion de la matriz A en dos matrices A = L ·U ,siendo L triangular inferior y U triangular superior. Esto permite, por ejemplo,calcular de forma sencilla el determinante de A:

|A| = |L||U | = (

n∏i=1

lii)(

n∏i=1

uii).

Ademas, podemos resolver el sistema de ecuaciones siguiendo estos dos pasos:Paso 1. Resolvemos el sistema Ly = b y obtenemos una solucion y.Paso 2. Resolvemos el sistema Ux = y. La solucion x obtenida sera la

solucion de Ax = b porque

Ax = LUx = Ly = b.

Los metodos de factorizacion LU consisten en descomponer la matriz enproducto de dos matrices triangulares, para despues aplicar los pasos 1 y 2anteriores y ası resolver el sistema.

Si A = L · U tendremos quea11 a12 . . . a1na21 a22 . . . a2n

. . .an1 an2 . . . ann

=

l11 0 0 . . . 0l21 l22 0 . . . 0

. . .ln1 ln2 . . . lnn

u11 u12 . . . u1n0 u22 . . . u2n

. . . u3n0 0 . . . unn

.

Haciendo el producto de matrices e igualando termino a termino se obtienenn2 ecuaciones lineales con n2 + n incognitas. Si imponemos la condicion

l11 = l22 = · · · = lnn = 1

obtendremos un sistema de ecuaciones facilmente resoluble. Las soluciones danlugar a una descomposicion L·U que es la misma que la obtenida por el metodode Gauss sin intercambio de filas.

Si A es no singular, una condicion necesaria y suficiente para que A admitauna descomposicion L·U de la forma anterior es que los menores fundamentales(todos los determinantes Ak de las matrices formadas con las primeras k filasy columnas) sean no nulos. En este caso la factorizacion es unica.

52 Capıtulo 3

Ejemplo: Resuelve el siguiente sistema mediante factorizacion L · U deGauss 3 1 2

6 3 2−3 0 −8

·x1x2x3

=

012

.

Como A1 = 3, A2 =

∣∣∣∣ 3 16 3

∣∣∣∣ = 0 y A3 =

∣∣∣∣∣∣3 1 26 3 2−3 0 −8

∣∣∣∣∣∣ = 0, el sistema

admite la factorizacion L·U buscada. Igualando termino a termino las matricessiguientes 3 1 2

6 3 2−3 0 −8

=

1 0 0l21 1 0l31 l32 1

u11 u12 u130 u22 u230 0 u33

.

obtenemos

u11 = 3u12 = 1u13 = 2

l21u12 = 6l21u12 + u22 = 3l21u13 + u23 = 2

l31u11 = −3l31u12 + l32u22 = 0

l31u13 + l32u23 + u33 = −8

⇒

u11 = 3u12 = 1u13 = 2l21 = 2u22 = 1u23 = −2l31 = −1l32 = 1u33 = −4

Paso 1.

Ly = b⇒

1 0 02 1 0−1 1 1

·y1y2y3

=

012

⇒ y1 = 0y2 = 1y3 = 1

.

Paso 2.

Ux = y ⇒

3 1 20 1 −20 0 −4

·x1x2x3

=

011

⇒ x1 = 0x2 = 1

2x3 = −1

4

,

que es la solucion del sistema.


3.3.1. Metodo de Crout

Si en la factorizacion LU fijamos

u11 = u22 = · · · = unn = 1

se obtiene la descomposicion por el metodo de Crout, que tiene el inconvenienteque puede producir grandes errores en la resolucion del sistema de ecuaciones.

3.3.2. Metodo de Cholesky

Si en la factorizacion LU exigimos que U = Lt de modo que

l11 = u11, l21 = u12, l31 = u13, . . . , lij = uji, . . . , lnn = unn,

se obtiene la descomposicion por el metodo de Cholesky que funciona si y solosi la matriz inicial A es simetrica y definida positiva. Por tanto, si consigo ladescomposicion de una matriz por el metodo de Cholesky, puedo asegurar quedicha matriz es simetrica y definida positiva.

La resolucion de un sistema de n = 50 ecuaciones con el metodo de Cho-lesky requiere 503/3 operaciones.

3.3.3. Sistemas triangulares

Son sistemas de ecuaciones lineales en los que la matriz de los coeficienteses de la forma:

A =

a11 a12 0 0 0 . . . 0a21 a22 a23 0 0 . . . 00 a32 a33 a34 0 . . . 0

. . .0 0 . . . 0 ann−1 ann

.

Aparecen frecuentemente en la resolucion numerica de ecuaciones diferencialesy en la aproximacion por splines cubicos. En la mayorıa de los casos, la matrizadmite una descomposicion de la forma

1 0 0 . . . 0l21 1 0 . . . 00 l32 1 0

. . .0 0 . . . lnn−1 1

u11 u12 0 . . . 00 u22 u23 0 00 0 u32 u33 . . . 0

. . .0 0 . . . unn−1 unn

.

54 Capıtulo 3

Ejemplo: Dada la matriz 60 30 2030 20 1520 15 12

a) Obten la descomposicion de Cholesky.b) Obten la factorizacion de Gauss sin intercambio de filas y, a partir de ella,las descomposiciones de Crout y de Cholesky.

a) 60 30 2030 20 1520 15 12

= LLt =

l11 0 0l21 l22 0l31 l32 l33

l11 l21 l310 l22 l320 0 l33

por tanto

l211 = 60l11l21 = 30l11l31 = 20l221l

222 = 20

l21l31 + l22l32 = 15l231 + l232 + l233 = 12

⇒

l11 =√

60

l21 =√

60/2

l31 =√

60/3

l22 =√

5

l32 =√

5

l33 =√

3/3

b) La factorizacion de Gauss se obtiene resolviendo el sistema de ecuacionesque se deduce de60 30 20

30 20 1520 15 12

=

1 0 0l21 1 0l31 l32 1

u11 u12 u130 u22 u230 0 u33

.

La descomposicion que resulta es

A =

60 30 2030 20 1520 15 12

=

1 0 01/2 1 01/3 1 1

60 30 200 5 50 0 1/3

.

Descomponiendo la segunda matriz de la parte derecha de la igualdad en unproducto de una matriz diagonal por otra triangular superior con 1 en cada


elemento de la diagonal, se obtiene que

A =

1 0 01/2 1 01/3 1 1

60 0 00 5 00 0 1/3

1 1/2 1/30 1 10 0 1

= (∗)

Multiplicando las dos primeras matrices anteriores

(∗) =

60 0 030 5 020 4 1/3

1 1/2 1/30 1 10 0 1

que es la descomposicion de Crout.

Si en (∗) descomponemos la matriz diagonal

(∗) =

1 0 01/2 1 01/3 1 1

√

60 0 0

0√

5 0

0 0√

1/3

√

60 0 0

0√

5 0

0 0√

1/3

1 1/2 1/30 1 10 0 1

=

=

√

60 0 0√60/2

√5 0√

60/3√

5√

3/3

√

60√

60/2√

60/3

0√

5√

5

0 0√

3/3

que es la descomposicion de Cholesky.

3.4. Normas y analisis del error

Definicion: Sea V un espacio vectorial. Una norma sobre V es una apli-cacion ∥ · ∥ : V → R tal que ∀v, w ∈ V y ∀α ∈ R se verifica que:

1. ∥v∥ ≥ 0 (v = 0⇔ ∥v∥ = 0)

2. ∥v + w∥ ≤ ∥v∥+ ∥w∥

3. ∥αv∥ = |α|∥v∥.

56 Capıtulo 3

Ejemplo: las siguientes son normas definidas sobre el espacio vectorialR3:

1. Norma euclıdea o norma 2: ∥(x, y, z)∥2 =√

x2 + y2 + z2

2. Norma del maximo o norma infinito: ∥(x, y, z)∥∞ = max(|x|, |y|, |z|)

3. Norma 1: ∥(x, y, z)∥1 = |x|+ |y|+ |z|

La norma de un vector mide la ”distancia”de ese vector al origen. La

distancia entre dos vectores se mide con la norma del vector diferencia. Lasnormas del ejemplo son diferentes formas de ”medir”las distancias entre losvectores de R3.

Puesto que el conjunto Mn×n de las matrices cuadradas de orden n, conlas operaciones suma y producto por un numero real es un espacio vectorial,para ”medir”matrices tambien se utilizaran normas que, por tanto, tendranque verificar las condiciones de la definicion de norma. Si ademas de esascondiciones se verifica que

∥A ·B∥ ≤ ∥A∥∥B∥ A,B ∈Mn×n

entones dicha norma es una norma matricial.Una norma matricial ∥ · ∥M sobre Mn×n y una norma vectorial ∥ · ∥V sobre

Rn se dice que son compatibles si se verifica que

∥Av∥V ≤ ∥A∥M∥v∥V ∀A ∈Mn×n∀v ∈ Rn.

Dada cualquier norma vectorial ∥ · ∥V sobre Rn, es posible definir una normamatricial sobre Mn×n de la siguiente forma:

∥A∥M = max∥Av∥V : ∥v∥V = 1

Esta nueva norma se llama norma matricial inducida o subordinada ala norma vectorial correspondiente. La norma vectorial y su norma matricialsubordinada son siempre compatibles.


Ejemplos:1. La norma matricial inducida por la norma vectorial 1 en R3 es

∥A∥1 =

∥∥∥∥∥∥a11 a12 a13a21 a22 a23a31 a32 a33

∥∥∥∥∥∥1

= max1≤j≤n

n∑i=1

|aij |.

De modo que

∥∥∥∥∥∥ 3 1 2

6 3 2−3 0 −8

∥∥∥∥∥∥1

= max(12, 4, 12) = 12.

2. La norma matricial inducida por la norma vectorial 2 en R3 es

∥A∥2 =√

ρ(At ·A).

ρ(At ·A) es el radio espectral de At ·A que es el maximo de los valores absolutosde los autovalores de la matriz At ·A.

Para calcular ∥A∥2 =

∥∥∥∥∥∥2 3 0

0 −1 00 0 1

∥∥∥∥∥∥2

debemos calcular primero los au-

tovalores de ∥At ·A∥ que son las raıces de

|AtA− αId| =

∣∣∣∣∣∣13− α −3 0−3 1− α 00 0 1− α

∣∣∣∣∣∣ = (1− α)(α2 − 14α + 4).

Los autovalores de AtA son 7±√

5 y 1 y como consecuencia

∥A∥2 =

∥∥∥∥∥∥2 3 0

0 −1 00 0 1

∥∥∥∥∥∥2

=

√7 +√

5

3. La norma matricial inducida por la norma infinito en R3 es

∥A∥∞ =

∥∥∥∥∥∥a11 a12 a13a21 a22 a23a31 a32 a33

∥∥∥∥∥∥∞

= max1≤i≤n

n∑j=1

|aij |.

58 Capıtulo 3

De modo que

∥∥∥∥∥∥ 3 1 2

6 3 2−3 0 −8

∥∥∥∥∥∥∞

= max(6, 11, 11) = 11.

4. Se define la norma matricial de Frobenius como ∥A∥ =√∑

1≤i,j≤n |ai,j |2.

3.4.1. Numero condicion de una matriz

Se define el numero de condicion de una matriz A como

k(A) = ∥A∥∥A−1∥.

El numero de condicion siempre es mayor o igual que 1 y se utiliza para estimarsi un sistema de ecuaciones Ax = b esta bien o mal condicionado, es decir, sipequenos cambios en los datos (A o b) pueden producir grandes cambios en lasolucion del sistema.

Supongamos que en vez de b utilizamos b de modo que la solucion delsistema es x y Ax = b. Esto puede ser debido a que hay un pequeno errorcometido al hacer alguna medida en b o bien porque en vez de la solucionexacta x hemos obtenido una solucion aproximada x. Si denotamos e al erroren la solucion y r al error en b o error residual, se verifica que Ae = Ax−Ax =b− b = r.

Teorema:1

k(a)

∥r∥∥b∥≤ ∥e∥∥x∥≤ k(a)

∥r∥∥b∥

.

Demostracion. Veamos que 1k(a)

∥r∥∥b∥ ≤

∥e∥∥x∥ . Se verifica que

∥r∥∥x∥ = ∥Ae∥∥A−1b∥ ≤ ∥A∥∥e∥∥A−1∥∥b∥

Por tanto1

∥A∥∥A−1∥∥r∥∥b∥≤ ∥e∥∥x∥

Vemos que ∥e∥∥x∥ ≤ k(a)∥r∥∥b∥ . Se verifica que

∥e∥∥b∥ ≤ ∥A−1r∥∥Ax∥ ≤ ∥A−1∥∥r∥∥A∥∥x∥.


Por tanto∥e∥∥x∥≤ ∥A∥∥A−1∥∥r∥

∥b∥.

Si k(A) estan cerca de 1 el sistema estara bien condicionado porque el error

relativo en la solucion ∥e∥∥x∥ sera similar al error relativo ∥r∥

∥b∥ en los datos. Si por

el contrario k(A) es muy grande, el sistema puede estar mal condicionado ypequenos cambios en los datos podran producir grandes cambios en la solucion.

Ejemplo: Dado el sistema

3,02 −1,05 2,534,33 0,56 −1,78−0,83 −0,54 1,47

x1x2x3

=

235

1. Encuentra el numero de condicion de la matriz A de los coeficientes con∥ · ∥∞.

2. Acota el error relativo de las posibles soluciones en funcion del error enlos datos.

a) La matriz inversa redondeada a dos decimales es

A−1 =

5,66 −7,27 −18,55200,51 −268,26 −669,91−76,85 −102,65 −255,88

con lo que

k(a) = ∥A∥∥A−1∥ = 6,67 · 1138,68 = 7595.

b) Se verifica que

1

7595

∥r∥5≤ ∥e∥∥x∥≤ 7595

∥r∥5

es decir

2,63 · 10−5∥r∥ ≤ ∥e∥∥x∥≤ 1591∥r∥

60 Capıtulo 3

3.5. Mejora de soluciones

3.5.1. Refinamiento iterativo

Es un metodo para mejorar una solucion aproximada ya obtenida. Sea xuna solucion aproximada de Ax = b tal que Ax = b. Se verifica que

x = A−1b = x+x− x = x+A−1b+A−1b = x+A−1(b− b) = x+A−1r = x+e.

El metodo consiste en aplicar los siguientes pasos:1. Calculo r = b−Ax con doble precision.2. Calculo e resolviendo el sistema Ae = r.3. La nueva solucion aproximada sera ˜x = x + e.

Los pasos anteriores se repiten las veces necesarias para mejorar la solucion.

Ejemplo: Si como solucion aproximada del sistema420 210 140 105210 140 105 84140 105 84 70105 84 70 60

x1x2x3x4

=

875539399319

se obtiene x = (0,999988, 1,000137, 0,999670, 1,000215), aplica el metodo derefinamiento iterativo para mejorarla.1. Calculo r = b−Ax con doble precision.

r =

875539399319

−

420 210 140 105210 140 105 84140 105 84 70105 84 70 60

0,9999881,0001370,9996701,000215

=

−1,05000000 · 10−4

−7,0000000 · 10−5

−3,5000000 · 10−5

−4,8000000 · 10−5

.

2. Resuelvo 420 210 140 105210 140 105 84140 105 84 70105 84 70 60

e1e2e3e4

=

−1,05 · 10−4

−7,0 · 10−5

−3,5 · 10−5

−4,8 · 10−5

y obtengo e = (1,1 · 10−5,−1,4 · 10−5, 3,30 · 10−4,−2,15 · 10−4)

3. La nueva solucion es ˜x = (0,999999, 0,999997, 1, 1).


3.5.2. Escalamiento

Si en un sistema de ecuaciones los coeficientes de las incognitas son de muydiferente magnitud, la solucion numerica del sistema puede contener errores.

Ejemplo: Resolvamos el siguiente sistema por eliminacion gaussiana:

2x1 + 100000x2 = 100000x1 + x2 = 2(

2 100000 1000001 1 2

)(fila 2)− 1/2 (fila 1)

(

2 100000 1000000 −999991 −99998

)De modo que las soluciones son

x1 = 100000−100000x22 ≈ 0

x2 = −99998−99999 ≈ 1( redondeado a 4 cifras)

Sin embargo el resultado correcto es x1 = 1,00002 y x2 = 0,99998.

Para solucionar este problema se puede recurrir al escalamiento que con-siste en multiplicar cada ecuacion por un numero para que el coeficiente masgrande de las incognitas en valor absoluto sea 1. Al aplicar el escalamientopodemos cometer a su vez errores de redondeo, pero puede ayudar en casosextremos cuando hay mucha diferencia entre los coeficientes.

Ejemplo: En el caso anterior

0,00002x1 + x2 = 1x1 + x2 = 2

redondeando a 4 decimales

x2 = 1x1 + x2 = 2

aplicando el pivoteox1 + x2 = 2

x2 = 1

y las soluciones son x1 = y x2 = 1.

62 Capıtulo 3

3.6. Metodos iterativos

Los metodos directos estudiados anteriormente requieren aproximadamen-te n3/3 operaciones. Son sensibles a errores de redondeo que se acrecientanal aumentar n. De modo que, aunque teoricamente conducen a una solucionexacta, en la practica la solucion obtenida puede ser peor que la obtenidaaproximadamente por un metodo iterativo.

Los metodos iterativos estan especialmente indicados para la resolucionde sistemas con una matriz de gran dimension pero dispersa, es decir, conmuchos ceros, que suelen aparecer por ejemplo en la resolucion de ecuacionesdiferenciales en derivadas parciales.

Un metodo iterativo de resolucion de un sistema de n ecuaciones Ax = b esaquel a partir de un x0 ∈ Rn genera una sucesion de posibles soluciones apro-ximadas x1, x2, . . . . El metodo es convergente si la sucesion generada convergea la solucion del sistema a partir de cualquier vector inicial x0. Es consistentesi, en caso de converger, el lımite es la solucion.

Todos las metodos iterativos de resolucion de un sistema Ax = b se basanen la descomposicion de la matriz A en diferencia de dos matrices M − N ,siendo M una matriz regular. De modo que

Ax = b⇔ (M −N)x = b⇔Mx = Nx + b⇔ x = M−1(Nx + b)

Si definimos G(x) = M−1(Nx+b), veremos que para encontrar la solucion delsistema basta encontrar un vector fijo de la funcion G(x). Para ello se generauna sucesion xk = G(xk−1) partiendo de una determinada solucion aproximadax0. La sucesion convergera a la solucion en determinadas condiciones.

En la practica, dadas las matrices M y N , puesto que

xk = G(xk−1) = M−1(Nxk−1) + b)⇒Mxk = Nxk−1 + b,

para calcular xk a partir de xk−1 se resuelve el sistema

Mxk = Nxk−1 + b.

Es esencial que M sea una matriz regular y es conveniente que sea una matrizsencilla para facilitar el calculo.

La condicion que asegura la convergencia es la siguiente.


Teorema: En las condiciones anteriores, si existe una norma matricialtal que ∥M−1N∥ < 1, entonces la sucesion xk+1 = G(xk) converge a un puntofijo de la funcion G(x).

Demostracion. Teniendo en cuenta que

∥G(y)−G(y′)∥ = ∥M−1(Ny + b)−M−1(Ny′ + b)∥= ∥M−1Ny −M−1Ny′∥ = ∥M−1N(y − y′)∥≤ ∥M−1N∥∥y − y′∥.

Si ∥M−1N∥ < λ < 1 para un cierto λ ∈ R, entonces ∥G(y)−G(y′)∥ < λ∥y−y′∥.A partir de un vector x0 vamos generando la sucesion xk+1 = G(xk) y severifica que

∥xk+1 − xk∥ = ∥G(xk)−G(x)k−1∥ =∥G(xk − xk−1)∥ <λ∥xk − xk−1∥ < . . . < λk∥x1 − x0∥.

Por tanto

∥xk+p − xk∥ < ∥xk+p − xk+p−1 + xk+p−1 − . . .− xk∥

≤p∑

i=1

∥xk+i − xk+i−1∥ <p−1∑i=0

λk+i−1∥x1 − x0∥

= ∥x1 − x0∥p−1∑i=0

λk+i−1 = ∥x1 − x0∥λk−1 − λk+p−1

1− λ

k→∞→ 0.

De modo que la sucesion xk+1 = G(xk) es una sucesion convergente a un vectorx ∈ Rn si ∥M−1N∥ < λ < 1 para un cierto λ ∈ R. Dicho vector x verifica que

x−G(x) = lımk→∞

xk−G( lımk→∞

xk) = lımk→∞

G(xk−1)−G( lımk→∞

xk)

= G( lımk→∞

(xk−1 − xk)) = G(0) = 0⇒ x = G(x).

64 Capıtulo 3

En resumen, un metodo iterativo de resolucion de un sistema Ax = bconsiste en generar una sucesion xk = G(xk−1) siendo G(x) = M−1(Nx + b)y A = M −N . Hay que tener en cuenta lo siguiente:

1. Si |M | = 0, puedo calcular xk a partir de xk−1.2. Si el sistema Mx = c se resuelve facilmente, el metodo iterativo es mas

rapido.3. Si ∥M−1N∥ < 1 para alguna norma matricial, la sucesion xk converge

a la solucion del sistema Ax = b.

3.6.1. Metodo de Jacobi

Dada la matriz A = (aij)1≤i,j≤n, se definen:

D =

a11 0 · · · 00 a22 0 · · · 0

· · ·0 · · · 0 ann

, L =

0 · · · 0a21 0 · · · 0a31 a32 0 · · · 0

· · ·an1 an2 · · · ann−1 0

U =

0 u12 · · · u1n0 0 u23 · · · u2n· · ·

0 · · · 0 un−1n

0 · · · 0 0

En el metodo de Jacobi M = D y N = −(L + U). Para calcula xk a partir dexk−1 se resuelve el sistema

Dxk = −(L + U)xk−1 + b.

Para aplicar el metodo se requiere que aii = 0 para todo i. Si |A| = 0, esto sepuede conseguir intercambiando filas si fuese necesario.

El metodo converge si ∥M−1N∥ < 1 para alguna norma matricial. Puestoque

M−1N =

0 −a12

a11−a13a11

· · · −a1na11−a21

a220 −a23

a22· · · −a2n

a22· · ·

−an1ann

· · · −ann−1

ann0

y


∥M−1N∥∞ = max1≤i≤n

(∑nj=1,j =i |aij |

)1aii

,

una condicion suficiente para que el metodo de Jacobi converja es que la matrizA sea diagonal dominante, es decir que |aii| >

∑nj=1,j =i |aij | ∀i = 1, 2, . . . , n.

Como

∥M−1N∥1 = max1≤j≤n∑n

i=1,i=j

(|aij | 1aii

),

otra condicion suficiente es que lo anterior sea menor que 1.Por ultimo, una tercera condicion suficiente es que la matriz D − L − U

sea simetrica y definida positiva.

Ejemplo: Aplica una iteracion del metodo de Jacobi al siguiente sistema6 5 00 2 12 0 3

x1x2x3

=

010

.

La matriz es diagonal dominante de modo que la sucesion xk de vectoresgenerados mediante el metodo de Jacobi convergera a la solucion. Se tiene que

D =

6 0 00 2 00 0 3

− (L + U) =

0 −5 00 0 −1−2 0 0

Si elegimos como vector inicial el x0 = (0, 0, 0), habra que resolver el sistema6 0 0

0 2 00 0 3

x11x12x13

=

0 −5 00 0 −1−2 0 0

000

+

010

Por tanto

6x11 = 02x12 = 13x13 = 0

⇒x11 = 0x12 = 1

2x13 = 0

El siguiente vector serıa x1 = (0, 1/2, 0).

66 Capıtulo 3

3.6.2. Metodo de Gauss-Seidel

En este metodo M = D+L y N = −U . El sistema a resolver para calcularxk a partir de xk−1 sera (D + L)xk = −Uxk−1 + b. Distintas condicionessuficientes para que el metodo funcione son:

1. que A sea diagonal dominante2. que At sea diagonal dominante3. que A sea simetrica y definida positiva.El metodo de Gauss-Seidel es dos veces mas rapido que el de Jacobi.

Ejemplo: Aplica el metodo de Gauss-Seidel al siguiente sistema6 2 00 5 12 0 3

x1x2x3

=

010

.

Tendremos que

D + L =

6 0 00 5 02 0 3

U =

0 2 00 0 10 0 0

.

Si partimos de x0 = (0, 0, 0), habra que resolver el sistema6 0 00 5 02 0 3

x11x12x13

=

0 −2 00 0 −10 0 0

000

+

010

y se obtiene que x1 = (0, 1/5, 0)


3.6.3. Metodos de relajacion

En el metodo de Jacobi, para calcular cada componente de xk se utili-zan todos los valores de xk−1. Sin embargo, en el metodo de Gauss-Seidel,para calcular la componente i-esima de xk (xki ) se utilizan las componen-tes xk1, x

k2, . . . , x

ki−1 del vector xk y las componentes xk−1

i+1 , . . . , xk−1n del vector

xk−1. De este modo, el calculo de las componentes de xk se divide en etapas,utilizando en cada una de ellas un vector diferente:

Etapa 1: calculo xk1 usando (xk−12 , xk−1

3 , . . . , xk−1n ).

Etapa 2: calculo xk2 usando (xk1, xk−13 , . . . , xk−1

n ).

Etapa 3: calculo xk3 usando (xk1, xk2, x

k−14 , . . . , xk−1

n ).

Etapa n: calculo xkn usando (xk1, xk2, . . . , x

kn−1).

En cada etapa se manejan tan solo vectores de n − 1 variables en lugar den con el consiguiente ahorro en tiempo. Ademas, la informacion obtenida encada etapa se incorpora a la siguiente. Los metodos que utilizan esto ultimo sellaman ”metodos de relajacion sucesiva 2se desarrollaron para la resolucion desistemas con matrices de dimension grande pero con casi todos los elementosnulos.

En un metodo de relajacion se considera

M =1

wD + L N =

1− w

wD − U

con lo que M −N = A. w se llama factor de relajacion. En funcion de su valorse tiene

Subrelajacion si 0 < w < 1.

Super-relajacion si w > 1

Gauss-Seidel si w = 1.

Si w ∈ (0, 2) el metodo no converge. Una condicion suficiente para la conver-gencia es que A sea simetrica y definida positiva.

Capıtulo 4

Aproximacion de funciones

4.1. Aproximacion alrededor de un punto

4.1.1. Polinomio de Taylor

Dada una funcion f(x), ¿como construir una funcion sencilla que ”se pa-rezca” a f(x) alrededor de un cierto punto a ∈ R?

Sean f(x) y g(x) dos funciones reales definidas alrededor de un puntoa ∈ R. Diremos que f(x) y g(x) tienen un contacto de orden r ≥ 0 en a,si se verifica que:

1. f(a) = g(a)

2. lımx→a

f(x)− g(x)

(x− a)r= 0.

Intuitivamente, que dos funciones tenga un punto de contacto en a significaque las dos funciones son ”parecidas” alrededor de a. El parecido sera mayorcuanto mayor sea el orden r del contacto. Si tienen un contacto de orden men a, entonces tiene contacto para cualquier orden r ≤ m.

69

70 Capıtulo 4

-1 -0.5 0.5 1 1.5 2

1

2

3

4

5

6

7

Figura 4.1: ex y sus polinomios de Taylor en 0 de orden 0, 1, 2 y 3.

Teorema 4.1. (local de Taylor) Sea B una bola abierta centrada en unpunto a ∈ R, y f : B ⊂ R → R es una funcion diferenciable m veces en a.Existe una unica funcion polinomica de grado menor o igual que m que tieneun contacto en a con f(x) de orden m. Dicha funcion es:

Pma f(x) = f(a) + f ′(a)(x− a) +

f ′′(a)

2!(x− a)2 + · · ·+ fm)(a)

m!(x− a)m

Pma f(x) es el polinomio de Taylor de f(x) de grado m en el punto a.

Demostracion. Sea P (x) = c0 + c1(x− a) + c2(x− a)2 + · · ·+ cm(x− a)m.Si f(x) y P (x) tienen un contacto de orden m en a, entonces se cumplira losiguiente:

1. f(a) = P (a)⇒ c1 = f ′(a).

2. lımx→a

f(x)− P (x)

(x− a)r= 0, ∀ 0 ≤ r ≤ m.

Aproximacion de funciones 71

En particular, si r = 1,

0 = lımx→a

f(x)− P (x)

(x− a)=

= lımx→a

f(x)− [f(a) + c1(x− a) + · · ·+ cm(x− a)m]

(x− a)=

lımx→a

f(x)− f(a)

x− a− c1 ⇒ c1 = f ′(a).

Si r = 2,

0 = lımx→a

f(x)− P (x)

(x− a)2=

= lımx→a

f(x)− [f(a) + c1(x− a) + · · ·+ cm(x− a)m]

(x− a)2=

= lımx→a

f(x)− f(a)− f ′(a)(x− a)

(x− a)2= (L’Hopital) =

= lımx→a

f ′(x)− f ′(a)

2(x− a)− c2 ⇒ c2 = f ′′(a)

2 .

Razonando de forma similar hasta r = m se obtiene que los unicos coeficientesposibles para el polinomio son los siguientes

c0 = f(a) c1 = f ′(a) c2 = f ′′(a)2! · · · cm = fm)(a)

m! .

Ejercicio 4.1. Calcula los polinomios de Taylor para la funciones ex ysen(x) en el punto 0 de orden 1, 2, 3, 4. Deduce la expresion para orden n.

Si dos funciones son derivables hasta un cierto orden m en un punto a ytiene un contacto de orden m en a, entonces las derivadas de ambas funcionesde cualquier orden r ≤ m valen lo mismo en a. Sin embargo, hay funcionesque tiene un contacto de orden m en un punto y no son derivables hasta ordenm en dicho punto.

Ejemplo 4.1. La funcion

f(x) =

x3 sen( 1x) si x = 0

0 si x = 0

y el polinomio P (x) = 0 tiene un punto de contacto de orden 3 en el punto 0.Sin embargo, la funcion f(x) no es 3 veces derivable en 0.

72 Capıtulo 4

Ejercicio 4.2. Si f(x) es un polinomio de grado n, ¿cuales son sus poli-nomios de Taylor de grado 0, 1, 2, . . . ,m?

4.1.2. Teorema de Taylor. Formula de Taylor

Ejercicio 4.3. Calcula aproximadamente el valor de e0,2 haciendo uso delos polinomios de Taylor de grado 1,2,3 y 4.

¿Como estimar el error cometido al aproximar un cierto valor haciendo usode los polinomios de Taylor? Si f(x) es m + 1 veces derivable en [a, b) y mveces derivable en (a, b), entonces existe c ∈ (a, b) tal que

f(b) = Pma f(b) + fm+1)(c)

(m+1)! (b− a)m+1.

La justificacion de lo anterior se basa en el siguiente teorema.

Teorema 4.2. Sean f(x) y g(x) dos funciones definidas en [a, b], continuasen [a, b],m veces derivables en [a, b) ym+1 veces derivables en (a, b). Entonces:

1. Existe c ∈ (a, b) tal que

[f(b)− Pma f(b)] · gm+1)(c) = [g(b)− Pm

a g(b)] · fm+1)(c)

2. Existe c ∈ (a, b) tal que

f(b)− Pma f(b) = fm+1)(c)

(m+1)! (b− a)m+1.

Demostracion. Caso 1. Consideramos las funciones

F (x) = f(x)− Pma f(x) G(x) = g(x)− Pm

a g(x).

A partir de ellas construimos la funcion H(x) = F (b)G(x)−G(b)F (x).

Paso 1. H(x) verifica las hipotesis del teorema de Rolle en el intervalo [a, b]:es continua en [a, b], es derivable en (a, b) y 0 = H(a) = H(b). Por tanto existeun punto c1 ∈ (a, b) tal que H ′(c1) = 0.


Paso 2. H ′(x) verifica las hipotesis del teorema de Rolle en el intervalo[a, c1]: es continua en [a, c1], es derivable en (a, c1) y 0 = H ′(a) = H ′(c1). Portanto existe un punto c2 ∈ (a, c1) tal que H ′′(c2) = 0.

Y ası sucesivamente hasta el paso m.

Paso m+1. Hm(x) verifica las hipotesis del teorema de Rolle en el intervalo[a, cm]: es continua en [a, cm], es derivable en (a, cm) y 0 = Hm)(a) = Hm)(cm).Por tanto existe un punto c ∈ (a, cm) tal que Hm+1(c) = 0. Pero

0 = Hm+1)(c) = F (b) ·Gm+1)(c)−G(b) · Fm+1)(c) =

[f(b)− Pma f(b)] · gm+1)(c)− [g(b)− Pm

a g(b)] · fm+1)(c).

Por lo tanto, hemos encontrado el punto c buscado.

Caso 2. Basta aplicar el caso 1 a f(x) y a g(x) = (x− a)m+1.

La expresion

f(b)− Pma f(b) =

fm+1)(c)

(m + 1)!(b− a)m+1

sirve para estimar la diferencia entre el valor de f(x) en un punto b cuandolo aproximamos por el valor del polinomio de Taylor en b, es decir, se utilizapara acotar el error cometido en la aproximacion.

Ejercicio 4.4. Calcula el error cometido al aproximar e0,2 haciendo usode los polinomio de Taylor de grado 1,2,3 y 4.

Ejercicio 4.5. Encuentra un polinomio que sirva para aproximar el valorde e1 con seis cifras decimales exactas.

Ejercicio 4.6. Encuentra un polinomio para aproximacion el valor sen(3)con diez cifras decimales exactas.

Para obtener el polinomio de Taylor de orden 4 alrededor del punto 0 deuna funcion f(x), con SAGE podemos usar el comando f(x).taylor(x,0,4).

74 Capıtulo 4

4.2. Aproximacion simultanea en varios puntos

Dada una tabla de datos como esta

x −1 0 1 2

y −2 −2 0 4

que representa unos valores y en funcion de otros x, pretendemos encontraruna funcion y = f(x) que pase por esos puntos (xi, yi). Los puntos xi se llamannodos.

4.2.1. Polinomio interpolador

Teorema 4.3. Dada una tabla de n + 1 puntos (x0, y0), (x1, y1), . . . ,(xn, yn) con xi = xj , existe un unico polinomio Pn(x) de grado menor o igualque n de modo que

Pn(xi) = yi ∀i = 0, 1, . . . , n.

Dicho polinomio se llama polinomio interpolador de los puntos.

Demostracion. Consideremos un polinomio cualquiera de grado menor oigual que n :

Pn(x) = a0 + a1x + a2x2 + · · ·+ anx

n.

Puesto que queremos que Pn(xi) = yi ∀i = 0, 1, . . . , n tendremos que

a0 + a1x0 + a2x20 + · · ·+ anx

n0 = y0

a0 + a1x1 + a2x21 + · · ·+ anx

n1 = y1

· · ·a0 + a1xn + a2x

2n + · · ·+ anx

nn = yn

que es un sistema de n + 1 ecuaciones lineales con n + 1 incognitas en el queel determinante de la matriz de los coeficientes es distinto de 0 :∣∣∣∣∣∣∣∣

1 x0 x20 · · · xn01 x1 x21 · · · xn1· · ·

1 xn x2n · · · xnn

∣∣∣∣∣∣∣∣ = 0.


Por lo tanto, el sistema es compatible determinado y tiene una unica solucion.

Ejemplo 4.2. Demuestra que existe un polinomio que pasa por los puntosde la tabla

x −1 0 1 2

y −2 −2 0 4

Como son 4 puntos buscaremos un polinomio de grado menor o igual que3 de la forma Pn(x) = a0 + a1x + a2x

2 + a3x3. El sistema que resulta es

a0 + a1(−1) + a2(−1)2 + a3(−1)3 = −2a0 + a1(0) + a2(0)2 + a3(0)3 = −2a0 + a1(1) + a2(1)2 + a3(1)3 = −2a0 + a1(2) + a2(2)2 + a3(2)3 = 4

Que es un sistema compatible determinado y que, por tanto, tiene una unicasolucion, es decir, existe un unico polinomio de grado menor o igual que 3 quepasa por los puntos de la tabla.

Ejercicio 4.7. Resuelve con software informatico el sistema de ecuacionesanterior y comprueba que el polinomio buscado es P (x) = −2 + x + x2.

Metodo de Lagrange

Dados los puntos (x0, y0), (x1, y1), . . . , (xn, yn), el polinomio interpoladorPn(x) se puede calcular ası:

Pn(x) =n∑

i=0

Li(x)yi

donde

Li(x) =Πj =i(x− xj)

Πj =i(xi − xj)=

=(x− x0)(x− x1) . . . (x− xi−1)(x− xi+1) . . . (x− xn)

(xi − x0)(xi − x1) . . . (xi − xi−1)(xi − xi+1) . . . (xi − xn).

Se verifica que

Li(xi) = 1 Li(xj) = 0 ∀i = j

76 Capıtulo 4

de modo que

Pn(xj) =

n∑i=0

Li(xj)yi = yj .

Ejercicio 4.8. Con los datos del ejercicio anterior, comprueba con un soft-ware matematico que los polinomios de Lagrange son los siguientes y calculael polinomio interpolador en funcion de dichos polinomios.

L0(x) = x(x−1)(x−2)(−1)(−2)(−3) L1(x) = (x+1)(x−1)(x−2)

1(−1)(−2)

L2(x) = (x+1)x(x−2)2·1(−1) L3(x) = (x+1)x(x−1)

3·2·1

El polinomio interpolador de un conjunto de puntos se puede calcular direc-tamente con SAGE usando el metodo de lagrange. Para aplicarlo es necesarioprimero escribir R=PolynomialRing(RDF,’x’). Despues el polinomio interpo-lador se calcula como R.lagrange_polynomial(puntos), donde puntos es lalista de puntos a interpolar.

Metodo de Newton

Dados los puntos (x0, y0), (x1, y1), . . . , (xn, yn), se pueden calcular lasdiferencias divididas• de primer orden: [xi, xi+1] = yi+1−yi

xi+1−xi

• de segundo orden: [xi, xi+1, xi+2] = [xi+1,xi+2]−[xi,xi+1]xi+2−xi

· · ·• de orden k [xi, xi+1, . . . , xi+k] =

[xi+1,··· ,xi+k]−[xi,··· ,xi+k−1]xi+k−xi

Las diferencias divididas son invariantes frente a permutaciones. Por ejemplo:

[x0, x1, x2] = [x1, x0, x2].

La tabla siguiente permite construir de forma sencilla las diferencias divididas:

x0 y0 [x0, x1] [x0, x1, x2] [x0, x1, x2, x3]x1 y1 [x1, x2] [x1, x2, x3]x2 y2 [x2, x3]x3 y3. . .

El polinomio interpolador se construye a partir de las diferencias divididas ası:

Pn(x) = y0 + [x0, x1](x− x0) + [x0, x1, x2](x− x0)(x− x1) + . . .

· · ·+ [x0, x1, x2, . . . , xn](x− x0)(x− x1) . . . (x− xn−1).


Ejercicio 4.9. Comprueba con software matematico que las diferenciasdivididas de los datos anteriores son las siguientes y expresa el polinomiointerpolador en funcion de dichas diferencias divididas:

−1 −2 0 1 00 −2 2 11 0 42 4

Para obtener las diferencias divididas son SAGE, primero tecleamos lasentencia R = PolynomialRing(RDF, ’x’). Despues, calculamos la tabla dediferencias divididas para el conjunto de puntos mediante mediante la ordenR.divided_difference(puntos, full_table=True), que devuelve una listade vectores. Si tomamos las primeras posiciones de cada vector obtenemos losvalores de y (primera columna), si tomamos las segundas posiciones obtenemoslas diferencias divididas de primer orden (segunda columna en nuestra tablade diferencias divididas) y ası sucesivamente.

4.3. Error del polinomio interpolador

Dados los puntos (x0, y0), (x1, y1), . . . , (xn, yn) que provienen de una ciertafuncion y = f(x), si conocemos cierta informacion sobre las derivadas de lafuncion f(x), podemos estimar el error que se comete al aproximar f(x) porsu polinomio interpolador Pn(x) de grado menor o igual que n.

Teorema 4.4. Sea f(x) una funcion n + 1 veces diferenciable en un in-tervalo (a, b) y sea Pn(x) su polinomio interpolador en los nodos x0, x1, . . . , xncontenidos en [a, b]. Para cada x ∈ [a, b] existe un punto θx ∈ (a, b) tal que

f(x)− Pn(x) =fn+1)(θx)

(n + 1)!(x− x0)(x− x1) . . . (x− xn).

Demostracion. Sea w(t) = (t − x0)(t − x1) . . . (t − xn). Para un valor fijox ∈ [a, b] existe un λ ∈ R tal que f(x)−Pn(x) = λw(x). De modo que tenemosla funcion

Φ(t) = f(t)− Pn(t)− λw(t)

que vale 0 en x, x0, x1, . . . , xn. Por tanto, Φ1)(t) se anulara en n + 1 puntos,Φ2)(t) se anulara en n puntos y Φn+1)(t) se anulara en 1 punto θx ∈ (a, b):

Φn+1)(θx) = fn+1)(θx)− Pn+1)(θx)− λwn+1)(θx) = 0

78 Capıtulo 4

es decir

fn+1)(θx)− λ(n + 1)! = 0⇒ λ =fn+1)(θx)

(n + 1)!.

Como consecuencia

f(x)− Pn(x) =fn+1)(θx)

(n + 1)!w(x).

Ejemplo 4.3. El error al aproximar la funcion f(x) = sen(x) en un puntox ∈ [0, 1] mediante un polinomio interpolador en 9 nodos contenidos en elintervalo [0, 1] se puede acotar ası:

|f(x)− P9(x)| ≤ |f10)(θx)|10!

|Π9i=0(x− xi)| ≤

1

10!< 2,8 · 10−7

Ejercicio 4.10. Elige 9 puntos del intervalo [0, 1]; obten el polinomio in-terpolador para f(x) = sen(x) en dichos puntos; calcula el valor de dichopolinomio interpolador en otro punto de [0, 1]; calcula el error, es decir, ladiferencia entre el valor exacto de la funcion y el valor en el polinomio inter-polador; comprueba que se verifica la cota de error prevista por el teorema enel ejemplo anterior.

4.3.1. Eleccion de nodos. Polinomios de Chebyshev

Si pretendemos encontrar un polinomio que interpole a una determinadafuncion y los nodos xi no estan previamente fijados, ¿que nodos podemoselegir? Elegir los nodos proporcionados por los polinomios de Chebyshev tienenalguna ventaja.

Definicion 4.1. Los polinomios de Chebyshev se definen ası:

T0(x) = 1 T1(x) = x Tn+1(x) = 2xTn(x)− Tn−1(x) si n ≥ 1

Proposicion 4.5. Las raıces del polinomio de Chebyshev Tn(x) si n ≥ 1son los valores

cos

(2j − 1

2nπ

)∀1 ≤ j ≤ n.


Ejercicio 4.11. Calcula el polinomio de Chebyshev de grado 10 y suscorrespondientes raıces. Comprueba la proposicion anterior.

Como consecuencia, todas las raıces de cualquier polinomio de Chebyshevestan contenidas en el intervalo [−1, 1].

Teorema 4.6. Si los nodos (xi), con i = 0, 1, . . . , n, de interpolacion deuna funcion f(x) en x ∈ [−1, 1] son las raıces del polinomio de ChebyshevTn+1(x), entonces

|f(x)− Pn(x)| ≤ 1

2n(n + 1)!max|t|≤1|fn+1)(t)|.

La justificacion del teorema se basa en que la expresion

max|x|≤1

Πni=0|(x− xi)|

que proviene de la formula general del error en la interpolacion, se hace mıni-ma cuando xi son las raıces del polinomio de Chebyshev correspondiente. Elresultado del teorema permite encontrar una cota del error del polinomio in-terpolador mas pequena que la establecida con caracter general.

Si en lugar del intervalo [−1, 1] tenemos un intervalo [a, b], si aplicamos alos nodos de Chebyshev la transformacion afın que lleva el intervalo [−1, 1] enel intervalo [a, b], obtenemos unos nuevos nodos tales que:

|f(x)− Pn(x)| ≤(b−a2

)n+1

2n(n + 1)!maxt∈[a,b]

|fn+1)(t)|.

De modo que el uso de los nodos que proporcionan los raıces de los polino-mios de Chebyshev tiene como consecuencia la obtencion de mejores resultadosen la interpolacion, debido a que los errores de aproximacion, en general, sonmenores.

Ejemplo 4.4. Calcula los cuatro nodos para la interpolacion de una fun-cion en el intervalo [−1, 1] y en el intervalo [2, 8] haciendo uso de los polino-mio de Chebyshev. Establece una cota del error en cada caso para la funcionf(x) = sin(x).

Calculamos las raıces de T4(x) = 8x4 − 8x2 + 1 :

80 Capıtulo 4

x0 = −0,9239 x1 = −0,3827 x2 = 0,3827 x3 = 0,9239.

Estos serıa los nodos en el intervalo [−1, 1]. Para calcular los nodos en elintervalo [2, 8] necesitamos una transformacion afın αx+β que lleve el intervalo[−1, 1] en el intervalo [2, 8] :

[−1, 1] → [2, 8]−1 7→ α(−1) + β = 2

1 7→ α(1) + β = 8

Resolviendo el sistema formado por las dos ultimas ecuaciones, obtenemosque la transformacion es 3x + 5. Aplicamos esta transformacion a los nodosde Chebyshev y obtenemos los nuevos nodos en el intervalo [2, 8] :

x′0 = 2,2284 x′1 = 3,8519 x′2 = 6,1481 x′3 = 7,7716.

Veamos las cotas de los errores para la funcion f(x) = sin(x) en cada caso.Para el intervalo [−1, 1] :

|f(x)− P (x)| ≤ 1

2n(n + 1)!max|t|≤1|fn+1)(t)| ≤ 1

23(4)!

Para el intervalo [2, 8] :

|f(x)− Pn(x)| ≤(b−a2

)n+1

2n(n + 1)!maxt∈[a,b]

fn+1)(t) ≤ 34

234!.

¿Es cierto que si n aumenta, entonces el polinomio Pn(x) se ”parece”mas a

la funcion f(x)? En general la respuesta es que no. Por ejemplo, dada funcionf(x) = 1

1+x2 en [−5, 5], si elegimos los nodos xi igualmente espaciados, se puededemostrar que los polinomios no convergen a la funcion fuera del intervalo[−3,63, 3,63]. (Ver figura 4.2)

Otro ejemplo es la funcion f(x) = |x| en [−1, 1], en la cual la convergenciase produce solo en los puntos −1, 0, 1.

Sin embargo, si elegimos los nodos de Chebyshev y la funcion f(x) a in-terpolar es continua y con derivada continua, entonces Pn(x) converge uni-formemente a f(x) en todo el intervalo. No obstante, para cualquier eleccion


-4 -2 2 4

0.5

1

1.5

2

-4 -2 2 4

-0.5

0.5

1

1.5

2

Figura 4.2: 11+x2 y sus polinomios interpoladores de grado 10 y 20.

de nodos es posible encontrar una funcion continua tal que no se produce esaconvergencia uniforme.

Con SAGE podemos obtener los polinomios de Chebyshev con el comandochebyshev_T(n,x), donde n es el grado del polinomio.

4.4. Interpolacion a trozos y con condiciones sobrela derivada

4.4.1. Interpolacion a trozos

Si disponemos de ciertos valores de una funcion f(x) en un intervalo [a, b],se puede obtener una aproximacion a la funcion f(x) dividiendo el interva-lo en varios intervalos mas pequenos (subintervalos) y utilizando polinomiosde pequeno grado (1,2,3, etc.) distintos para interpolar los valores de cadasubintervalo. La siguiente grafica muestra la funcion seno y una funcion quela interpola con trozos de recta en varios puntos:

Ejercicio 4.12. Dada la tabla

x −1 0 1 2

y −2 −2 0 4

encuentra una funcion que pase por los puntos anteriores formada pora) por trozos de rectab) por trozos de parabola

82 Capıtulo 4

1 2 3 4 5 6

-1

-0.5

0.5

1

4.4.2. Interpolacion con condiciones sobre la derivada

Los siguientes son ejemplos de problemas de interpolacion con condicionessobre la derivada.

Ejemplo 4.5. Encuentra un polinomio P (x) tal que P (0) = 0, P (1) = 1y P ′(0,5) = 2.

Puesto que son tres las condiciones que tenemos, cabe pensar que basta buscarentre los polinomios con tres parametros, es decir, los polinomios de grado 2de la forma P (x) = a + bx + cx2 :

P (0) = 0P (1) = 1

P ′(0,5) = 2⇒

a = 0b + c = 1

b + 2 · c · 0,5 = 2⇒

a = 0b + c = 1b + c = 2

El sistema anterior no tiene solucion, con lo que habrıa que buscar entre lospolinomios de grado 3. Al plantear las ecuaciones correspondientes resulta unsistema con solucion no unica.

Ejercicio 4.13. Comprueba lo afirmado en el ejemplo anterior mediantesoftware matematico.

Ejemplo 4.6. (Interpolacion de Hermite): Encuentra un polinomio P (x)que interpole a una funcion y = f(x) sabiendo que


x −1 0 1

y 1 1 −1y′ 2 0 0

Puesto que hay seis condiciones para la funcion y sus derivadas buscamosentre los polinomios de la forma P (x) = a0 + a1x + · · · + a5x

5 que tienenseis coeficientes. Planteamos el sistema correspondiente segun las condicionesy obtenemos lo siguiente:

P (−1) = 1 P ′(−1) = 2P (0) = 1 P ′(0) = 0P (1) = 0 P ′(1) = 0

⇒a0 = 1 a3 = −3a1 = 0 a4 = 1

2a2 = −3

2 a5 = 2

Con lo que el polinomio buscado es P (x) = 1− 32x

2 − 3x3 + 12x

4 + 2x5.


Ejemplo 4.7. (Interpolacion de Hermite cubica a trozos): Encuentra unafuncion formada por trozos de polinomios de grado 3 que verifique las siguien-tes condiciones

x −1 0 1

y 1 1 −1y′ 2 0 0

Buscamos una funcion de la forma

g(x) =

P (x) si x ∈ [−1, 0]Q(x) si x ∈ [0, 1]

donde P (x) = a + bx + cx2 + dx3 y Q(x) = a′ + b′x + c′x2 + d′x3. A partir delas condiciones planteamos el sistema

P (−1) = 1 P ′(−1) = 2 Q(0) = 1 Q′(0) = 0P (0) = 1 P ′(0) = 0 Q(1) = −1 Q′(1) = 0

con lo que

g(x) =

1 + 2x2 + 2x3 si x ∈ [−1, 0]1− 6x2 + 4x3 si x ∈ [0, 1]

(Ver Figura 4.3)

84 Capıtulo 4

-1 -0.5 0.5 1

-1

-0.5

0.5

1

Figura 4.3: La funcion g(x)


Ejemplo 4.8. (Splines cubicos): Encuentra una funcion polinomica cubi-ca a trozos que sea de clase C2 (continua, derivable y con derivada segundacontinua) que verifique las siguientes condiciones

x −1 0 1

y 1 1 −1

Buscamos una funcion de la forma

h(x) =

P (x) si x ∈ [−1, 0]Q(x) si x ∈ [0, 1]

donde P (x) = a+ bx+ cx2 + dx3 y Q(x) = a′ + b′x+ c′x2 + d′x3 de modo que

P (−1) = 1 P (0) = 1 Q(0) = 1 Q(1) = −1 P ′(0) = Q′(0) P ′′(0) = Q′′(0)

Puesto que tenemos 6 ecuaciones y 8 incognitas necesitamos 2 condicionesmas. Por ejemplo, anadimos las siguientes

P ′′(−1) = 0 Q′′(1) = 0


-1 -0.5 0.5 1

-1

-0.5

0.5

1

Figura 4.4: La funcion h(x)

y obtenemos un sistema compatible determinado de modo que

h(x) =

1− x− 3

2x2 − 1

2x3 si x ∈ [−1, 0]

1− x− 32x

2 + 12x

3 si x ∈ [0, 1]

(Ver Figura 4.4)


Con SAGE podemos manejar los spline de un conjunto de puntos te-cleando previamente from sage.gsl.all import spline. Despues el splinese calcula mediante la orden S=spline(puntos), que posteriormente nos per-mitira calcular los valores del spline. Por ejemplo, S(2.3) nos darıa el valordel spline en 2.3.

Interpolacion de Hermite

Un problema de interpolacion de Hermite consiste en encontrar un polino-mio P (x) (Polinomio de Hermite) tal que

P (x0) = y0 P ′(x0) = y′0 P ′′(x0) = y′′0 · · · P k0−1)(x0) = yk0−10

P (x1) = y1 P ′(x1) = y′1 P ′′(x1) = y′′1 · · · P k1−1)(x1) = yk1−11

· · ·P (xn) = yn P ′(xn) = y′n P ′′(xn) = y′′n · · · P kn−1)(xn) = ykn−1

n

86 Capıtulo 4

siendo yi, y′i, · · · , y

k)i los valores de una cierta funcion y sus derivadas en n+ 1

puntos. Si el numero de condiciones impuestas es k1 + k2 + · · ·+ kn = m + 1entonces

Teorema 4.7. Existe un unico polinomio de grado menor o igual que mque satisface las condiciones de interpolacion de Hermite anteriores.

Demostracion. Sea P (x) un polinomio de grado menor o igual que m. Lascondiciones P j(xi) = yji dan lugar a un sistema de m + 1 ecuaciones de laforma Au = y. Este sistema tiene solucion unica si el sistema Au = 0 tieneunicamente la solucion 0, es decir, si el sistema formado por las condicionesP j)(xi) = 0 tiene como solucion unica el polinomio 0. Veamos que esto ultimoes cierto. Un polinomio P (x) que sea solucion del sistema P j)(xi) = 0 verificaque tiene la raız xi repetida ki veces, con lo que si P (x) no es cero, serıa degrado mayor o igual que m + 1, en contradiccion con la hipotesis inicial.

Ejercicio 4.17. ¿Que otro nombre recibe el polinomio de Hermite cuandosolo hay un nodo?

Teorema 4.8. Sean x0, x1, . . . , xn nodos distintos de [a, b] y f ∈ C2[a, b].Si P (x) es el polinomio de grado menor o igual que 2n+1 tal que P (xi) = f(xi)y P ′(xi) = f ′(xi), con 0 ≤ i ≤ n, entonces ∀x ∈ [a, b] ∃θ ∈ [a, b] tal que

f(x)− P (x) =f2n+2)(θ)

(2n + 2)!Πn

i=0(x− xi)2

El polinomio de Hermite se puede calcular usando las diferencias divididas.Para ello definimos

[x0, x0] = limx→x0 [x0, x] = lımx→x0

y − y0x− x0

= lımx→x0

f(x)− f(x0)

x− x0= f ′(x0)

[x0, x0, x1] =[x0, x1]− [x0, x0]

x1 − x0

[x0, x0, x0] =1

2!f2)(x0) [x0, x0, x0, x0] =

1

3!f3)(x0)

Y de forma similar el resto de diferencias. En estas condiciones el polinomiode interpolacion de Hermite se calcula de forma parecida a como se calculabael polinomio de interpolacion sin condiciones sobre la derivada.


El calculo es mas sencillo haciendo uso de una tabla. Por ejemplo, si elproblema de Hermite incluye datos sobre la funcion y su derivada en trespuntos la tabla quedarıa ası:

x y 10 20 30 40 50

x0 y0 [x0, x0] [x0, x0, x1] [x0, x0, x1, x1] [x0, x0, x1, x1, x2] [x0, . . . , x2]x0 y0 [x0, x1] [x0, x1, x1] [x0, x1, x1, x2] [x0, x1, x1, x2, x2]x1 y1 [x1, x1] [x1, x1, x2] [x1, x1, x2, x2]x1 y1 [x1, x2] [x1, x2, x2]x2 y2 [x2, x2]x2 y2

siendo[x0, x0, . . . , x0] = 1

k!fk)(x0)

k veces.

El polinomio serıa

P (x) = y0 + [x0, x0](x− x0) + [x0, x0, x1](x− x0)2+

+ [x0, x0, x1, x1](x− x0)2(x− x1)+

+ [x0, x0, x1, x1, x2](x− x0)2(x− x1)

2

+ [x0, x0, x1, x1, x2, x2](x− x0)2(x− x1)

2(x− x2).

Ejemplo 4.9. Calcula el polinomio que verifica estas condiciones:

x −1 0 1

y 1 1 −1

y′ 2 0 0

En este caso la tabla quedarıa ası

x y 10 20 30 40 50

-1 1 2 −2 2 −3/2 2-1 1 0 0 −1 5/20 1 0 −2 40 1 −2 21 -1 01 -1

El polinomio serıa

88 Capıtulo 4

P (x) = 1+2(x+1)−2(x+1)2+2(x+1)2x−3/2(x+1)2x2+2(x+1)2x2(x−1)

Ejemplo 4.10. Calcula el polinomio que verifica las siguientes condicio-nes:

x 1 2

y 2 6

y′ 3 7

y′′ 8

En este caso la tabla quedarıa ası

x y 10 20 30 40

1 2 3 1 2 −11 2 4 3 12 6 7 8/2!=42 6 72 6

El polinomio serıa

P (x) = 2 + 3(x− 1) + 1(x− 1)2 + 2(x− 1)2(x− 2)− 1(x− 1)2(x− 2)2.

Interpolacion cubica de Hermite a trozos

Dada la tabla de valores

x x0 x1 · · · xny y1 y2 · · · yny′ y′1 y′2 · · · y′n

correspondientes a una cierta funcion y sus derivadas en los nodos xi, el metodode interpolacion cubica de Hermite a trozos consiste en encontrar una unicafuncion derivable definida a trozos, tal que en cada intervalo [xi, xi+1] (i =0, . . . , n− 1), la funcion es un polinomio Pi(x) de grado menor o igual a 3 queverifica lo siguiente

Pi(xi) = yi Pi(xi+1) = yi+1 P ′(xi) = y′i P ′i (xi+1) = y′i+1.

Planteando los correspondientes sistemas de ecuaciones, se puede encontrar launica funcion que interpola en esas condiciones.


Splines

Una funcion spline de grado k con nodos x0 ≤ x1 ≤ · · · ≤ xn es una funcionS(x) definida a trozos tal que:

a) S(x) es un polinomio Si(x) de grado menor o igual que k en cada inter-valo [xi, xi+1], con i = 0, 1, . . . , n− 1.

b) S(x) admite derivada segunda continua de orden k − 1 en [x0, xn].

El spline cubico (k = 3) es el mas usado en interpolacion. Un spline cubico deinterpolacion en x0, x1, . . . , xn es una funcion que cumple las propiedadesa), b) y S(xi) = f(xi), con i = 0, 1, . . . , n. El numero de condiciones que sehan de cumplir es el siguiente:

condicion ecuaciones num. ec.

continuidad Si−1(xi) = Si(xi) con i = 1, . . . , n− 1 n− 1derivada 1a continua S′

i−1(xi) = S′i(xi) con i = 1, . . . , n− 1 n− 1

derivada 2a continua S′′i−1(xi) = S′′

i (xi) con i = 1, . . . , n− 1 n− 1interpolacion Si(xi) = yi con i = 0, . . . , n n + 1

total 4n− 2

Cada polinomio cubico Si(x), con i = 0, . . . , n − 1, tendra 4 parametros,con lo que el numero de incognitas sera 4n−4. Por lo tanto, para poder resolverel correspondiente sistema de ecuaciones hemos de anadir dos condiciones masy podemos hacerlo de distintas formas:

Si exigimos S′′(x0) = 0 y S′′(xn) = 0 obtendremos el spline cubiconatural.

Si exigimos S′(x0) = S′(xn) y S′′(x0) = S′′(xn) obtendremos el splineperiodico.

Se pueden elegir otras condiciones distintas del spline natural o del periodico.Sin embargo, los splines natural y periodico existen y son unicos para cualquierconjunto de nodos y de datos. Ademas, si la funcion a interpolar f(x) tienederivada segunda continua en [x0, xn], entonces∫ b

a[S′′(x)]2dx ≤

∫ b

a[f ′′(x)]2dx,

90 Capıtulo 4

siendo S(x) el spline natural. Es decir, de entre todas las funciones con derivadasegunda continua que pasan por ciertos puntos (x0, y0), (x1, y1), . . . , (xn, yn) elspline natural es el que tiene menor energıa elastica (menor tension1).

Por otra parte, si aumentamos el numero de puntos a interpolar en elintervalo [a, b] (n→∞) disminuyendo simultaneamente el tamano maximo delos subintervalos [xi, xi+1] (max0≤i≤n[xi, xi+1]→ 0), entonces el spline naturaly su derivada convergen uniformemente a f(x) y f ′(x) respectivamente.

4.5. Aproximacion por mınimos cuadrados

(Material didactico obtenido de

http://www.ugr.es/ lorente/APUNTESMNQ/capitulo4/ajuste5.pdf )

La busqueda de una funcion que represente lo mejor posible a unos datosexperimentales (xi, yi) (i = 1, . . . , n) puede abordarse, entre otras, de las dosformas siguientes:

1. Mediante interpolacion (expuesta en las secciones anteriores), es decir,obtencion de una funcion y = f(x) que pase exactamente por todos losdatos (xi, yi);

2. mediante la obtencion de una curva, y = f(x), que se aproxime a losdatos sin que, necesariamente, pase por ellos.

El caso 2 es el que a continuacion abordamos. Para precisar el metodo autilizar hemos de tener en consideracion dos aspectos de importancia:

¿Que tipo de funcion usaremos para aproximar?

¿Cual es el criterio de aproximacion que utilizaremos?

En la practica para aproximar se usan funciones polinomicas de grado bajo, obien de tipo exponencial, potencial, e incluso funciones tipo spline o polinomiosa trozos. En cuanto al criterio a considerar para obtener el modelo concretoconsiste, geometricamente, en hacer mınima la suma de los cuadrados de lasdiferencias (distancia euclıdea) entre los valores yi experimentales y los valoresf(xi) (que denotaremos tambien yi) que tiene la funcion en el mismo puntoxi.

1La tension o energıa elastica en [a, b] es proporcional a∫ b

a[f ′′(x)]2dx.


El problema quedarıa planteado en los terminos siguientes: dados los datos(xi, yi) : i = 1, 2, . . . , n, hallar f(x) que verifique que la siguiente suma

n∑i=1

(yi − f(xi))2 (4.1)

es mınima cuando buscamos entre las funciones f(x) de una cierta clase.La expresion 4.1 dependera de unos parametros a, b, . . . segun la clase de

funciones f(x) elegida. Una solucion del problema anterior recibe el nombrede ajuste mınimos cuadrados de los datos.

Dependiendo del modelo de funcion f(x) utilizado aparecen distintos tiposde ajuste.

1. Modelos lineales: aquellos que utilizan funciones f(x) de la forma:

f(x) = a1q1(x) + a2q2(x) + · · ·+ amqm(x)

donde ai (i = 1, . . . ,m, y normalmente m < n) son los parametros adeterminar en el problema de ajuste, y qi(x) (i = 1, . . . ,m) son funcioneslinealmente independientes de cierto espacio de funciones (polinomios,spline, etc.).

Como casos particulares de ajustes de este tipo tenemos:

a) Ajuste lineal (recta de mınimos cuadrados) si utilizamos: f(x) =a+ bx. Encontrar la funcion f(x) que hace mınima la expresion 4.1para este caso supone buscar el mınimo de la expresion F (a, b) =∑n

i=1(yi − a − bxi)2 que depende de a y de b. Dicho mınimo se

alcanza en un punto donde las derivadas de F (a, b) respecto de ay respecto de b se anulan, lo que conlleva a plantear un sistema dedos ecuaciones lineales en las incognitas a y b. Su solucion son losvalores de a y b que minimizan 4.1.

b) Ajuste polinomial si utilizamos el modelo: f(x) = a0 + a1x + · · ·+akxk (llamado cuadratico para k = 2 y cubico para k = 3).

2. Modelos no lineales: aquellos que utilizan funciones f(x) que no sonlineales respecto los parametros del ajuste. Algunos modelos clasicosson:

a) Ajuste exponencial: f(x) = aebx.

92 Capıtulo 4

b) Ajuste potencial: f(x) = axb.

c) Ajuste racional: f(x) = 1a+bx .

Ejercicio 4.18. Encuentra el ajuste lineal, cuadratico y cubico para lossiguientes datos:

x −1 0 1 2

y −2 −2 0 4

La funcion de SAGE find_fit encuentra el mejor ajuste por mınimoscuadrados de un conjunto de datos (xi, yi) para un cierto modelo f(x) especi-ficado.

4.6. Aproximacion de funciones periodicas. Seriesde Fourier

4.6.1. Funciones periodicas. Funciones pares e impares

Sea f : A ⊂ R→ R. Se dice que

f(x) es periodica si existe T > 0 tal que f(x + T ) = f(x) ∀x ∈ A. Ental caso, el periodo de la funcion es el menor valor de T que verifica loanterior.

f(x) es par si f(−x) = f(x) ∀x ∈ A,

f(x) es impar si f(−x) = −f(x) ∀x ∈ A.

Las funciones sen(x), cos(x) y tan(x) son periodicas de perido 2π, 2π y π,respectivamente. La funcion cos(x) es par. Las funciones sen(x) y tan(x) sonimpares.

Propiedades

1. Si f(x) es una funcion impar, entonces∫ a−a f(x)dx = 0 ∀a ∈ R.

2. Si f(x) es una funcion par, entonces∫ a−a f(x)dx = 2

∫ a0 f(x)dx ∀a ∈ R.

3. Si f(x) es una funcion periodica de periodo 2l, entonces∫ l−l f(x)dx =

∫ a+2la f(x)dx ∀a ∈ R.


-2 -1 1 2

-1

-0.5

0.5

1

-2 -1 1 2

-1

-0.5

0.5

1

Integral de una funcion impar

-2 -1 1 2

-0.75

-0.5

-0.25

0.25

0.5

0.75

1

-2 -1 1 2

-0.75

-0.5

-0.25

0.25

0.5

0.75

1

Integral de una funcion par

-2 2 4 6 8

-1

-0.5

0.5

1

-2 2 4 6 8

-1

-0.5

0.5

1

Integral de una funcion periodica

4. Si f(x) es impar, entonces f(x) cos(nx) es impar y f(x) sen(nx) es par∀n ∈ Z.

5. Si f(x) es par, entonces f(x) cos(nx) es par y f(x) sen(nx) es impar∀n ∈ Z.

6. Si f(x) es periodica de periodo 2π, entonces f(x) cos(nx) y f(x) sen(nx)tambien son periodicas de periodo 2π ∀n ∈ Z.

7. Si f(x) es periodica de periodo 2l, entonces f(x) cos(nπl x) y f(x) sen(nπ

l x)tambien son periodicas de periodo 2l ∀n ∈ Z.

94 Capıtulo 4

4.6.2. Desarrollo de una funcion en serie de Fourier

Dada una funcion periodica f(x) de periodo T = 2l, pretendemos encontraruna serie de la forma:

S(x) =a02

+

∞∑n=1

an cos(nπ

lx)

+ bn sen(nπ

lx)

(4.2)

que se parezca a f(x) en el sentido en el que

1

2l

∫ a+2l

a[f(x)− S(x)]2 dx (4.3)

sea lo mınimo posible.

Cuando en la expresion 4.2 manejamos un numero finito de sumandos, ha-blamos de polinomios trigonometricos en lugar de serie. Es decir, una expresiondel tipo

S(x) =a02

+k∑

n=1

an cos(nπ

lx)

+ bn sen(nπ

lx)

se denomina polinomio trigonometrico, y k es el grado del polinomio. Lospolinomios trigonometricos que mejor aproximan (en el sentido de minimizarla integral 4.3) a una funcion periodica son los polinomios de Fourier.

-7.5 -5 -2.5 2.5 5 7.5

-3

-2

-1

1

2

3

-7.5 -5 -2.5 2.5 5 7.5

-3

-2

-1

1

2

3

Figura 4.5: f(x) = x con x ∈ [π, π] y su polinomio de Fourier de grado 3.


Para obtener los polinomios trigonometricos o el desarrollo en serie deFourier de una funcion f(x) periodica de perıodo 2l, los coeficientes an y bnse han de calcular ası:

a0 =1

l

∫ l

−lf(x)dx

an =1

l

∫ l

−lf(x) cos

(nπlx)dx

bn =1

l

∫ l

−lf(x) sen

(nπlx)dx.

Proposicion 4.9. Sea f(x) es una funcion periodica de perıodo 2l. Si

1. f(x) es par, entonces en su desarrollo en serie de Fourier, los coeficientesbn seran nulos.

2. f(x) es impar, entonces en su desarrollo en serie de Fourier, los coefi-cientes an seran nulos.

Desarrollo de una funcion no periodica

Sea f(x) una funcion cualquiera (no necesariamente periodica) en [a, b],monotona a trozos y acotada.

1. Si la completamos definiendola de forma periodica, podremos obtener suserie de Fourier.

2. Si la completamos definiendola de forma periodica y par, podremos obte-ner su serie de Fourier que solo tendra terminos con cosenos (desarrolloen cosenos de f(x)).

3. Si la completamos definiendola de forma periodica e impar, podremosobtener su serie de Fourier que solo tendra terminos con senos (desa-rrollo en senos de f(x)).

Teorema 4.10. (de Dirichlet): Sea f(x) una funcion periodica de pe-riodo 2l, monotona a trozos y acotada en (−l, l). En estas condiciones la seriede Fourier

S(x) =a02

+∞∑n=1

an cos(nπ

lx)

+ bn sen(nπ

lx)

96 Capıtulo 4

es convergente en todo punto x ∈ R. La suma de la serie obtenida es igual a

S(x) =

f(x) si f(x) es continua en x

lımt→x− f(x)+lımt→x+ f(x)2 si f(x)no es continua en x.

Teorema 4.11. (de mınimo desvıo cuadratico): Dada una funcionf(x) periodica de perıodo 2l, de entre todos los polinomios trigonometricos degrado k de la forma

Sk(x) =a02

+

k∑n=1

an cos(nπ

lx)

+ bn sen(nπ

lx),

el obtenido mediante los coeficientes de Fourier para f(x), hace mınimo eldesvıo medio cuadratico:

1

2l

∫ a+2l

a[f(x)− Sk(x)]2 dx

Con SAGE es posible definir un funcion f en un intervalo de amplitud T ycalcular con el comando f.fourier_series_partial_sum(k,T) su desarrollode Fourier de orden k considerando la funcion periodica de periodo T .

Ejercicio 4.19. Obten el polinomio de Fourier de grado 10 para las si-guientes funciones y representalo graficamente:

1. f(x) = x en el intervalo [0, π] y de periodo π.

2. f(x) = x2 en el intervalo [0, 2] y de period 2.

Ejercicio 4.20. Obten el desarrollo en senos y cosenos de las funcionesdel ejercicio anterior.

4.7. Practicas con SAGE

Practicas sobre el tema disponibles para los alumnos en

http://sage.unex.es

1. MUI_Interpolacion2. MUI_Interpolation and Aproximation

3. MUI_Fourier

-7.5 -5 -2.5 2.5 5 7.5

-3

-2

-1

1

2

3

f(x) = x si x ∈ [−π, π]

-7.5 -5 -2.5 2.5 5 7.5

-2

-1

1

2

Polinomio de Fourier de orden 1 def(x)

-7.5 -5 -2.5 2.5 5 7.5

-2

-1

1

2


-7.5 -5 -2.5 2.5 5 7.5

-3

-2

-1

1

2

3


-7.5 -5 -2.5 2.5 5 7.5

-3

-2

-1

1

2

3

Polinomio de Fourier de orden 10de f(x)

-7.5 -5 -2.5 2.5 5 7.5

-3

-2

-1

1

2

3

Polinomio de Fourier de orden 30de f(x)

Capıtulo 5

Diferenciacion e integracionnumerica

5.1. Diferenciacion numerica y extrapolacion de Ri-charson

¿Puede calcularse la derivada de una funcion f(x) en un punto si se conocensolo los valores de la funcion f(x0), f(x1),. . . ,f(xn) en n + 1 puntos? Si lafuncion es un polinomio de grado n, entonces seguro que se puede porquepodremos obtener la expresion explıcita del polinomio. Si de la funcion solosabemos que es derivable, entonces no sera posible, pues hay muchas funcionesderivables que pasan por los mismos puntos. El objetivo de esta seccion esestimar el valor de f ′(x) a partir de los valores de la funcion en ciertos puntosf(x0), f(x1),. . . ,f(xn) dando una cota del error cometido. Dado que

f ′(x) = lımh→0

f(x + h)− f(x)

h

una primera estimacion es

f ′(x) ≈ f(x + h)− f(x)

h

que es aplicable cuando conocemos los valores de la funcion en dos puntos:x0 = x y x1 = x+h. La formula anterior es exacta (es decir, el valor aproximadocoincide con el valor exacto) para funciones lineales (por ejemplo, para f(x) =

99

100 Capıtulo 5

3x + 2) y en otros casos de manera fortuita. Calculemos una cota del error.Por el teorema de Taylor:

f(x + h) = f(x) + f ′(x)h +f ′′(θ)

2!h2

es decir

f ′(x) =f(x + h)− f(x)

h− h

2f ′′(θ).

El error al aproximar f ′(x) por f ′(x) = f(x+h)−f(x)h , llamado error de trunca-

miento, es

error = −h

2f ′′(θ) (formula de orden h)

siendo θ un numero entre x y x + h. Cuanto menor sea h, mejor sera laaproximacion de f ′(x), hasta que el error de redondeo impida la mejora.

Ejemplo 5.1. Estima el valor de la derivada de cos(x) en el punto π4 y

calcula una cota del error cometido.

Tomando como h = 0,01

f ′(π

4) ≈

f(π4 + h)− f(π4 )

h=

cos(π4 + 0,01)− cos(π4 )

0,01= −0,71063050.

La cota del error serıa

|error| = | − h

2f ′′(θ)| = |0,01

2cos(θ)| ≤ 0,01

2= 0,005.

Ejercicio 5.1. Calcula la estimacion del ejemplo 5.1 anterior mediantesoftware matematico tomando h = 0,01. Calcula el error mediante la diferenciaentre la aproximacion y el valor exacto obtenido al derivar la funcion.

Otra estimacion de f ′(x) es

f ′(x) ≈ f(x + h)− f(x− h)

2h.

Veamos la estimacion del error. Del teorema de Taylor se deduce que

f(x + h) = f(x) + f ′(x)h + h2

2 f ′′(x) + h3

6 f ′′′(θ1)

f(x− h) = f(x)− f ′(x)h + h2

2 f ′′(x)− h3

6 f ′′′(θ2).

Diferenciacion e integracion 101

Si restamos las ecuaciones anteriores

f(x + h)− f(x− h) = 2f ′(x)h +h3

6[f ′′′(θ1) + f ′′′(θ2)].

Despejando

f ′(x) =f(x + h)− f(x− h)

2h− h2

12[f ′′′(θ1) + f ′′′(θ2)]

con lo que el error de truncamiento serıa

−h2

12[f ′′′(θ1) + f ′′′(θ2)]

siendo θ1 un numero entre x y x+h y θ2 un numero entre x y x−h. Si f ′′′(x)existe y es continua en [x−h, x+h], entonces existe un numero θ ∈ [x−h, x+h]

tal que f ′′′(θ) = f ′′′(θ1)+f ′′′(θ2)2 , con lo que error de truncamiento quedarıa

error = −h2

6f ′′′(θ) (formula de orden h2).

Ejercicio 5.2. Calcula de nuevo la estimacion del ejemplo 5.1 anteriormediante software matematico tomando h = 0,01 y usando la formula deorden h2. Calcula el error mediante la diferencia entre la aproximacion y elvalor exacto obtenido al derivar la funcion.

Para estimar el valor de una derivada segunda podemos usar la formula

f ′′(x) =f(x + h)− 2f(x) + f(x− h)

h2.

Veamos la estimacion del error. Por el teorema de Taylor

f(x + h) = f(x) + f ′(x)h + h2

2 f ′′(x) + h3

6 f ′′′(x) + h4

24f4)(θ1)

f(x− h) = f(x)− f ′(x)h + h2

2 f ′′(x)− h3

6 f ′′′(x) + h4

24f4)(θ2)

.

Si sumamos las ecuaciones anteriores y despejamos

f ′′(x) =f(x + h)− 2f(x) + f(x− h)

h2− h2

24[f4)(θ1) + f4)(θ2)].

102 Capıtulo 5

Siendo θ1 y θ2 numeros que estan entre x y x + h y entre x y x− h. Si f4)(x)existe y es continua, entonces existe un numero θ ∈ [x− h, x + h] tal que

error =h2

12f4)(θ) (formula de orden h2 )

Ejercicio 5.3. Calcula la estimacion de la derivada segunda de cos(x) enx = π/4 mediante software matematico tomando h = 0,01 y usando la formulade orden h2. Calcula el error mediante la diferencia entre la aproximacion y elvalor exacto obtenido al derivar la funcion.

5.1.1. Diferenciacion mediante interpolacion

Dada una funcion f(x) y recordando como se calcula el polinomio de in-terpolacion en los puntos x0, x1, . . . ,xn, se tiene que :

f(x) =n∑

i=0

f(xi)Lxi(x) +1

(n + 1)!fn+1)(θx)w(x)

donde w(x) = Πni=0(x− xi). Si derivamos

f ′(x) =n∑

i=0

[f(xi)L

′xi

(x) +1

(n + 1)!fn+1) (θx)w′(x)

+1

(n + 1)!w(x)

d

dx(fn+1)(θx))

]Si x = xk

f ′(xk) =

n∑i=0

[f(xi)L

′xi

(xk) +1

(n + 1)!fn+1)(θxk

)w′(xk)

].

Y como

w′(x) =

n∑i=0

Πnj=0,j =i(x− xj)⇒ w′(xk) = Πn

j=0,j =k(xk − xj)

nos queda

f ′(xk) =n∑

i=0

[f(xi)L

′xi

(xk)]

+1

(n + 1)!fn+1)(θxk

)Πnj=0,j =k(xk − xj)


que es una forma de calcular un valor aproximado de f ′(xk) a partir de losvalores x0, x1, . . . ,xn.

Ejemplo 5.2. a) Aplica la formula anterior para calcular f ′(x1) con n = 2.b) Calcula lo anterior con x1 − x0 = x2 − x1 = h

Apartado a). Puesto que

Lx0 =(x− x1)(x− x2)

(x0 − x1)(x0 − x2)Lx1 =

(x− x0)(x− x2)

(x1 − x0)(x1 − x2)

Lx2 =(x− x0)(x− x1)

(x2 − x0)(x2 − x1)

se tiene que

L′x0

=2x− x1 − x2

(x0 − x1)(x0 − x2)L′x1

=2x− x0 − x2

(x1 − x0)(x1 − x2)

L′x2

=2x− x0 − x1

(x2 − x0)(x2 − x1).

Por tanto

f ′(x1) ≈ f(x0)x1 − x2

(x0 − x1)(x0 − x2)+ f(x1)

2x1 − x0 − x2(x1 − x0)(x1 − x2)

+ f(x2)x1 − x0

(x2 − x0)(x2 − x1)

siendo el error de esta aproximacion

1

3!f3)(θx1)(x1 − x0)(x1 − x2).

Apartado b). Siendo x0 = x1 − h y x2 = x1 + h se obtiene la formula yaconocida

f ′(x1) ≈ f(x1 − h)−h

(−h)(−2h)+ f(x1)

2x1 − (x1 − h)− (x1 + h)

(h)(−h)

+ f(x1 + h)h

(2h)(h)

f(x1 + h)− f(x1 − h)

2h,

siendo el error de esta aproximacion

1

3!f3)(θx1)(h)(−h) = − 1

3!f3)(θx1)h2.

104 Capıtulo 5

5.1.2. Extrapolacion de Richarson

La extrapolacion de Richarson sirve para generar resultados de gran exac-titud cuando se usan formulas de bajo orden. La extrapolacion puede aplicarsesiempre a cualquier metodo de aproximacion en el que sepamos el termino deerror de una forma previsible y se basa en un parametro que generalmente esel tamano de paso h.

Por el teorema de Taylor sabemos que

f(x + h) =∞∑k=0

1

k!hkfk)(x) y f(x− h) =

∞∑k=0

1

k!(−1)khkfk)(x).

Restando las formulas anteriores

f(x + h)− f(x− h) = 2hf ′(x) +2

3!h3f3)(x) +

2

5!f5)(x) + . . .

y despejando f ′(x)

f ′(x) =f(x + h)− f(x− h)

2h−

[f3)(x)

3!h2 +

f5)(x)

5!h4 +

f7)(x)

7!h6 + . . .

].

Si definimos Φ(h) = f(x+h)−f(x−h)2h y denotamos ai al coeficiente de hi en la

formula anterior, tenemos

f ′(x) = Φ(h) + a2h2 + a4h

4 + a6h6 + . . . (5.1)

que es la formula ya conocida de aproximacion de f ′(x)

f ′(x) ≈ Φ(h) con error = h2[a2 + a4h

2 + a6h4 + . . .

]Si sustituimos h por h/2

f ′(x) = Φ(h

2) +

a222

h2 +a424

h4 +a626

h6 + . . .

y multiplicando por 4

4f ′(x) = 4Φ(h

2) + a2h

2 +a44h4 +

a624

h6 + . . . (5.2)


con lo que (5.2)-(5.1) es

3f ′(x) = 4Φ(h

2)− Φ(h)− 3a4

4h4 − 15a6

16h6

que, despejando f ′(x), da origen a la formula de aproximacion

f ′(x) ≈ 4

3Φ(

h

2)− 1

3Φ(h) con error = h4

[−a4

4− 5a6

16h2 − . . .

](5.3)

que tambien puede expresarse

f ′(x) ≈ Φ(h

2) +

Φ(h2 )− Φ(h)

3con error = h4

[−a4

4− 5a6

16h2 − . . .

]Ejemplo 5.3. Siendo f(x) = arctan(x), calcula f ′(

√2) haciendo uso de

la extrapolacion de Richarson.

Tomamos h = 0,1.

f ′(x) ≈ Φ(h) =f(x + h)− f(x− h)

2h

=f(√

2 + 0,1)− f(√

2− 0,1)

2 · 0,1= 0,3339506968

y

f ′(x) ≈ Φ(h/2) = 0,3334876594

Usando la extrapolacion

f ′(x) ≈ Φ(h/2) +Φ(h/2)− Φ(h)

3= 0,3333333136

(El valor exacto es 1/3).

Si en la formula (5.3) denotamos Ψ(h) = 43Φ(h2 ) − 1

3Φ(h), b4 = −a44 , b6 =

−5a616 , . . . obtenemos la expresion

f ′(x) = Ψ(h) + b4h4 + b6h

6 + . . . (5.4)

Sustituyendo h por h/2

106 Capıtulo 5

f ′(x) = Ψ(h/2) + b4

(h

2

)4

+ b6

(h

2

)6

+ . . .

y multiplicando por 24

16f ′(x) = 16Ψ(h/2) + b4h4 + b6

h6

4+ . . . (5.5)

Restando (5.5)-(5.4) nos queda

15f ′(x) = 16Ψ(h/2)−Ψ(h)− 3b6h6

4+ . . .

y al despejar f ′(x)

f ′(x) =16

15Ψ(h/2)− 1

15Ψ(h)− b6

h6

20+ . . .

da origen a la formula

f ′(x) ≈ 16

15Ψ(h/2)− 1

15Ψ(h) con error =h6

[−b6

1

20+ . . .

]que tambien puede expresarse

f ′(x) ≈ Ψ(h/2) +Ψ(h/2)−Ψ(h)

15con error =h6

[−b6

1

20+ . . .

]El proceso puede continuar de forma similar tomando ahora otra funcion

Ω(h) = 1615Ψ(h/2)− 1

15Ψ(h). Se puede demostrar el siguiente resultado:

Teorema 5.1. El algoritmo de Richarson genera aproximaciones a f ′(x)cada vez de un orden superior.

Ejemplo 5.4. Aplica dos pasos del algoritmo de Richarson para calcularf ′(√

2) siendo f(x) = arctan(x) y tomando h = 0,1.

Aplicamos la formula (5.1.2) con h y h/2 para obtener las dos primerasaproximaciones de orden 2

Φ(0,1) =arctan(

√2 + 0,1)− arctan(

√2− 0,1)

2 · 0,1= 0,33395069677432


Φ(0,05) =arctan(

√2 + 0,05)− arctan(

√2− 0,05)

2 · 0,05= 0,33348765942111

A partir de esos dos valores y aplicando la formula (5.3) calculamos la apro-ximacion Ψ(h) que es de orden 4

Ψ(0,1) = Φ(0,05) +Φ(0,05)− Φ(0,1))

3= 0,33333331363671

Calculando Φ(h/4) (Φ(0,025) = 0,33337191390106) podremos tambien calcu-lar Ψ(h/2)

Ψ(0,05) = Φ(0,025) +Φ(0,025)− Φ(0,05)

3= 0,33333333206104

y a partir de estos dos valores obtenemos la aproximacion Ω(h)

Ω(0,1) = Ψ(0,05) +Ψ(0,05)−Ψ(0,1)

15= 0,33333333328933.

En general, el proceso de extrapolacion de Richarson puede resumirse ası:

Etapa 0. Calculo distintas aproximaciones R(0, h) en funcion de un parame-

tro h (en nuestro caso R(0, h) = Φ(h) = f(x+h)−f(x−h)2h ).

Etapa 1. A partir de R(0, h) y de R(0, h/2) calculo las aproximaciones

R(1, h) = R(0, h/2) +R(0, h/2)−R(0, h)

41 − 1.

Etapa 2. A partir de R(1, h) y de R(1, h/2) calculo las aproximaciones

R(2, h) = R(1, h/2) +R(1, h/2)−R(1, h)

42 − 1.

Etapa n. A partir de R(n− 1, h) y de R(n− 1, h/2) calculo las aproxima-ciones

R(n, h) = R(n− 1, h/2) +R(n− 1, h/2)−R(n− 1, h)

4n − 1.

Etapa 0 Etapa 1 Etapa 2 Etapa 3 . . .

R(0, h)R(0, h/2) R(1, h)R(0, h/4) R(1, h/2) R(2, h)R(0, h/8) R(1, h/4) R(2, h/2) R(3, h)

. . . . . . . . . . . . . . .

108 Capıtulo 5

5.2. Integracion numerica mediante interpolacion

El objetivo es obtener un valor aproximado de∫ ba f(x)dx sin calcular la

primitiva de f(x). La idea es sustituir∫ ba f(x)dx por

∫ ba p(x)dx, siendo p(x)

una funcion que se parezca a f(x), por ejemplo un polinomio interpoladoren unos nodos x0, x1, . . . , xn del intervalo [a, b]. Si utilizamos el metodo deLagrange para calcular el polinomio interpolador

p(x) =

n∑i=0

f(xi)lxi(x) con lxi(x) =

n∏j=0

x− xjxi − xj

,

tendrıamos que∫ b

af(x)dx ≈

∫ b

ap(x)dx =

n∑i=0

f(xi)

∫ b

alxi(x)

es decir, obtenemos una formula como esta

∫ b

af(x)dx ≈

n∑i=0

Aif(xi). (5.6)

En el ejemplo anterior Ai =∫ ba lxi(x)dx. Una expresion como (5.6) se

denomina formula de cuadratura y los puntos xi son los nodos de cuadratura.Las formulas ası construidas son exactas para los polinomios de grado menor oigual que n (¿por que?). Si en el ejemplo utilizado los nodos estan igualmenteespaciados, estaremos hablando del metodo de Newton-Cotes.

5.2.1. Regla del trapecio

El caso particular del metodo de Newton-Cotes con n = 1 se conoce comoel metodo del trapecio. Puesto que

lx0(x) =b− x

a− x, lx1(x) =

x− a

b− a,

se tiene que

A0 =

∫ b

alx0(x)dx =

1

2(b− a)


-1 1 2

0.2

0.4

0.6

0.8

1

-1 1 2

0.2

0.4

0.6

0.8

1

-3 -2 -1 1 2

-1

-0.5

0.5

1

Figura 5.1: formulas de cuadratura graficamente

y

A0 =

∫ b

alx1(x)dx =

1

2(b− a).

Por tanto, la formula de cuadratura queda ası∫ b

af(x)dx ≈

n∑i=0

Aif(xi) =b− a

2[f(a) + f(b)]

La expresion del error del polinomio de interpolacion p(x) para n = 1 es

f(x)− p(x) = f ′′(ξx)2 (x− a)(x− b), de modo que∫ b

af(x)dx−

∫ b

ap(x)dx =

∫ b

af(x)− p(x)dx =

∫ b

a

f ′′(ξx)

2(x− a)(x− b)dx

Puesto que (x − a)(x − b) no cambia de signo en todo [a, b] y suponiendoque f ′′(x)/2 es una funcion continua en [a, b], aplicando la generalizacion delteorema de valor medio del calculo integral 1 obtenemos que el error es

f ′′(ξ)

2

∫ b

a(x− a)(x− b)dx =

f ′′(ξ)

2

(a− b)3

6

1Teorema: Sean h(x) y g(x) dos funciones reales e integrables en un intervalo [a, b]. Sig(x) no cambia de signo en todo x ∈ [a, b] y h(x) es continua en [a, b], entonces existe unpunto ξ ∈ [a, b] tal que: ∫ b

a

h(x)g(x)dx = h(ξ)

∫ b

a

g(x)dx,

.

110 Capıtulo 5

para un cierto valor ξ ∈ (a, b). En resumen, si h = b − a, la Regla delTrapecio quedarıa ası:∫ b

af(x)dx ≈ h

2[f(a) + f(b)] con error=− f ′′(ξ)

h3

12(5.7)

para un cierto ξ ∈ [a, b].

Ejemplo 5.5. Calcula∫ 0,20 ex

2dx con la formula de Newton-Cotes con n =

1 (regla del trapecio).

∫ 0,2

0ex

2dx ≈ 0,2− 0

2[e0 + e0,2

2] = 0,20408107741924

Sabiendo que si x ∈ [0, 0,2] entonces

|f ′′(x)| = |2ex2+4x2ex

2 | ≤ |f ′′(0,2)| = |2e0,22+4·0,2e0,22 | = 2,248151272255559,

una cota del error sera

| − 2,2481512722555590,23

12| = 0,00149876751484.

Ejercicio 5.4. Calcula de nuevo la estimacion del ejemplo 5.5 anteriormediante software matematico usando la regla del trapecio.

Ejercicio 5.5. Compara la estimacion anterior con las obtenidas con loscomandos de SAGE f.integrate( , , ).n( ), que obtiene la expresiondecimal de la integral exacta, y numerical_integral(f, , ), que obtiene unaproximacion numerica de la integral.

5.2.2. Regla de Simpson

La formula de Newton-Cotes para n = 2 o n = 3 se llama Regla de Simpson.El caso n = 2, h = b−a

2 y xi = x0 + ih, se conoce como la regla de Simpson1/3: ∫ b

af(x)dx ≈ h

3[f(x0) + 4f(x1) + f(x2)] con error− f4)(ξ)

h5

90


para un cierto ξ ∈ [a, b]. Por ser una formula de Newton-Cotes con tres puntoses exacta para polinomios de grado menor o igual que 2; pero ademas es exactapara los polinomios de grado menor o igual que 3, ya que el termino del errorcontiene el factor f4)(ξ).

El caso n = 3, h = b−a3 y xi = x0 + ih, se conoce como la regla de

Simpson 3/8:

∫ b

af(x)dx ≈ 3h

8[f(x0) + 3f(x1) + 3f(x2) + f(x3)] con error− f4)(ξ)

3h5

80

para un cierto ξ ∈ [a, b]. Es mejor usar la regla de Simpson 1/3 que la de3/8, puesto que ambas son de orden 5 y | − 1/90| < | − 3/80|. En generales preferible usar una regla con n par (numero impar de puntos) que con nimpar.

Ejemplo 5.6. Calcula∫ 0,20 ex

2dx con la formula de Newton-Cotes con n =

2 (regla de Simpson 1/3).

∫ 0,2

0ex

2dx ≈ 0,2/2

3[e0 + e0,1

2+ e0,2

2] = 0,2027

Sabiendo que si x ∈ [0, 0,2] entonces

|f4)(x)| = |12ex2

+ 48ex2x2 + 16ex

2x4| ≤ |f4)(0,2)

= |12e0,22

+ 48e0,220,22 + 16e0,2

20,24| = 14,5147,

una cota del error sera

| − 14,5147(0,2/2)5

90| = 1,6 ∗ 10−6

Ejercicio 5.6. Calcula de nuevo la estimacion del ejemplo 5.6 mediantesoftware matematico usando las reglas de Simpson. Compara la estimacion an-terior con las obtenidas con los comandos de SAGE f.integrate( , , ).n()

y numerical_integral(f, , ).

112 Capıtulo 5

5.2.3. Reglas compuestas

Este metodo consiste en dividir el intervalo [a, b] en subintervalos y apli-car en cada uno de ellos una regla simple. Por ejemplo, la regla del trapeciocompuesta consiste en fijar unos puntos a = x0 < x1 < . . . < xn−1 < xn = b yaplicar la regla del trapecio a cada intervalo [xi−1, xi]:∫ b

af(x)dx =

n∑i=1

∫ xi

xi−1

f(x)dx ≈ 1

2

n∑i=1

(xi − xi−1) [f(xi−1) + f(xi)] .

Es equivalente a sustituir f(x) por una funcion que interpole a f(x) y queeste formada por trozos de lınea. Si los puntos estan igualmente espaciados(h = xi − xi−1 = b−a

n , xi = a + ih con i = 1, 2, . . . , n), la regla del trapeciocompuesta nos queda

∫ b

af(x)dx ≈ h

2[f(x0) + 2f(x1) + · · ·+2f(xn−1) + f(xn)]

con error− h2

12(b− a)f ′′(ξ)

para un cierto ξ ∈ [a, b]. La formula es exacta para polinomios de gradomenor o igual que 1.

-1 1 2

0.2

0.4

0.6

0.8

1

-1 1 2

0.2

0.4

0.6

0.8

1

-1 1 2

0.2

0.4

0.6

0.8

1

Figura 5.2: Regla del trapecio compuesta


Ejemplo 5.7. Calcula un valor aproximado de∫ 0,20 ex

2dx mediante la re-

gla del trapecio compuesta.

Si tomamos tres puntos x0 = 0, x1 = 0,1 y x2 = 0,2∫ 0,2

0ex

2dx ≈ 0,1

2

[e0 + 2e0,1

2+ e0,2

2]

= 0,201085

Una cota del error serıa∣∣∣∣−0,12

12(0,2− 0)(ex

2(2 + 4ex

2))

∣∣∣∣ ≤ ∣∣∣∣−0,12

120,2 · f ′′(0,2)

∣∣∣∣ = 0,00187346

Ejercicio 5.7. Calcula de nuevo la estimacion del ejemplo 5.7 anteriormediante software matematico usando la regla del trapecio compuesta. Com-para la estimacion anterior con las obtenidas con los comandos de SAGEf.integrate( , , ).n() y numerical_integral(f, , ).

Si el numero n de intervalos en los que dividimos [a, b] es par podemosaplicar de forma sencilla la regla de Simpson compuesta con h = b−a

n , xi =a+ ih (0 ≤ i ≤ n), aplicando a cada dos subintervalos la regla de Simpson 1/3∫ b

af(x)dx =

∫ x2

x0

f(x)dx +

∫ x4

x2

f(x)dx + · · ·+∫ xn

xn−2

f(x)dx

con lo que obtenemos la formula

∫ b

af(x) ≈ h

3[f(x0) + 4f(x1) + 2f(x2) + 4f(x3) + 2f(x4) + . . .

+ 4f(xn−1) + f(xn)]

o lo que es lo mismo

∫ b

af(x)dx ≈ f(x0) +

n2∑

i=1

[2f(x2i)+4f(x2i−1)] + f(xn)

con error− b− a

180h4f4)(ξ)

114 Capıtulo 5

con ξ ∈ [a, b].

5.2.4. Metodo de los coeficientes indeterminados

Las reglas de Newton-Cotes con n + 1 puntos son exactas para los po-linomios de grado n, es decir, si en la formula de cuadratura

∫ ba f(x)dx ≈∑n

i=0Aif(xi) tomamos los coeficientes Ai =∫ ba lxi(x)dx, el resultado que ob-

tenemos al calcular con esa regla la integral de un polinomio de grado n es elvalor exacto. El objetivo de esta seccion es, dados un conjunto de n + 1 pun-tos x0, x1, . . . , xn, encontrar el valor de los coeficientes A0, A1, . . . , An queaseguren que la formula de cuadratura

∫ ba f(x)dx ≈

∑ni=0Aif(xi) es exacta

para los polinomios de grado n. En cada caso, el problema se puede resolverde forma similar al siguiente ejemplo.

Ejemplo 5.8. Encuentra los valores A0, A1 y A2 para que la formula∫ 1

0f(x)dx ≈ A0f(0) + A1f(

1

2) + A2f(1)

sea exacta para los polinomios de grado menor o igual que 2.

En particular la formula tendra que ser exacta para los polinomios 1, x y x2

que forma una base del conjunto de los polinomios de grado menor o igual que2:

1 =∫ 10 1dx = A0 + 1A1 + A2

12 =

∫ 10 xdx = 0A0 + 1

2A1 + A213 =

∫ 10 x2dx = 0A0 + 1

4A1 + A2

y obtenemos que A0 = 16 , A1 = 2

3 y A2 = 16 . Ademas, si P (x) = a + bx + cx2

es cualquier otro polinomio de grado menor o igual que 2, la formula serıatambien exacta porque∫ 1

0P (x)dx = a

∫ 1

01dx + b

∫ 1

0xdx + c

∫ 1

0x2dx =

= a(A0 + 1A1 + A2) + b(1

2A1 + A2) + c(

1

4A1 + A2) =

= A0(a) + A1(a +1

2b + (

1

2)2c) + A2(a + 1b + 12c) =

= A0P (0) + A1P (1

2) + A2P (1)


que serıa el resultado obtenido al aplicar la formula de cuadratura con loscoeficientes calculados. Ademas, es la misma formula que se obtendrıa con elmetodo de Simpson.

5.3. Cuadratura gaussiana

Supongamos que queremos obtener una formula de cuadratura∫ 1

−1f(x)dx ≈ A1f(t1) + A2f(t2)

pero no fijamos previamente los valores de t1, t2, A1 y A2. Para calcular estosvalores y puesto que tenemos 4 incognitas, le imponemos 4 condiciones: que laformula sea exacta para los polinomios de grado menor o igual que 3, es decir,para 1, x, x2 y x3 : ∫ 1

−1 x3dx = 0 = A1t

31 + A2t

32∫ 1

−1 x2dx = 2

3 = A1t21 + A2t

22∫ 1

−1 xdx = 0 = A1t1 + A2t2∫ 1−1 1dx = 2 = A1 + A2

Restando a la ecuacion primera la tercera multiplicada por t21

0 = A2(t32 − t21t2) = A2t2(t2 − t1)(t2 + t1)

se obtiene alguno de los siguientes valores: A2 = 0, t2 = 0, t2 = t1, t2 = −t1.Sustituyendo en el resto de las ecuaciones podemos comprobar que solo esaceptable t2 = −t1 y que en ese caso los parametros son A1 = A2 = 1 y

t2 = −t1 =√

13 . Por tanto, la formula

∫ 1

−1f(x)dx ≈ f(−

√1

3) + f(

√1

3)

es exacta para integrar cualquier polinomio de grado menor o igual que 3 enel intervalo [−1, 1]. Para integrar en [a, b] se utiliza un cambio de variablemediante una transformacion afın que lleve un intervalo en el otro.

116 Capıtulo 5

Los nodos ti (i = 1, 2, . . . , n) en la cuadratura gaussiana coinciden con lasraıces del polinomio de Legendre de grado n. Los polinomios de Legendre seconstruyen ası:

L0(x) = 1 L1(x) = x Ln+1(x) =(2n + 1)xLn(x)− nLn−1(x)

n + 1.

Los polinomios de Legendre se pueden obtener con SAGE con el comandogen_legendre_P(n, 0, x), siendo n el grado del polinomio.

Ejemplo 5.9. Calcula un valor aproximado de∫ π

20 sen(x)dxmediante cua-

dratura gaussiana de dos terminos.

Como el intervalo de integracion es [0, π2 ], hemos de hacer un cambio de variableafın que lleve [−1, 1] en [0, π2 ]. Dicho cambio es x = π

4 t+π4 con lo que dx = π

4dt.Ya hemos deducido cuales son los nodos ti y los coeficientes Ai (i = 1, 2). Si nosupieramos cuales son, podrıamos averiguar primero los nodos calculando lasraıces del polinomio de Legendre de grado 2 (L2(x) = 3x2−x

2 ). Posteriormentecalcularıamos los coeficientes Ai planteando un sistema de ecuaciones. En todocaso se obtendrıa∫ π

2

0sen(x)dx =

∫ 1

−1sen(π

4t +

π

4

) π

4dt ≈

≈ π

4

[sen

(π

4(−√

1

3) +

π

4

)+ sen

(π

4(

√1

3) +

π

4

)]≈ 0,99847.

Ejercicio 5.8. Calcula de nuevo la estimacion del ejemplo 5.9 mediantesoftware matematico usando la regla cuadratura gaussiana de dos terminos.Compara la estimacion anterior con las obtenidas con los comandos de SAGEf.integrate( , , ).n() y numerical_integral(f, , ).

5.4. Integracion de Romberg

Para aplicar el metodo de integracion de Romberg construimos una tablacomo esta


R(0, 0)

R(1, 0) R(1, 1)

R(2, 0) R(2, 1) R(2,2)

R(3, 0) R(3, 1) R(3,2) R(3,3)

. . .

R(M, 0) R(M, 1) R(M, 2) R(M, 3) . . . R(M,M)

donde, si queremos aproximar∫ ba f(x)dx

R(0, 0) =1

2(b− a) [f(a) + f(b)]

(Regla del trapecio para 20 intervalos)

R(n, 0) =1

2R(n− 1, 0) +

b− a

2n

2n−1∑i=1

f(a + (2i− 1)b− a

2n)

(Regla del trapecio para 2n intervalos)

R(n,m) = R(n,m− 1) +1

4m − 1[R(n,m− 1)−R(n− 1,m− 1)]

(Extrapolacion de Richarson)

(5.8)

La ultima formula se obtiene aplicando la extrapolacion de Richarson quees un proceso que ya hemos utilizado en la seccion 5.1.2. El factor 1

4m−1 de-pende del orden de convergencia del metodo empleado, en este caso la regladel trapecio.

Cada valor R(i, j) se considera una aproximacion a∫ ba f(x)dx. A partir

de la aproximacion R(n,m − 1) que se considera que mejora a la anteriorR(n− 1,m− 1), se construye otra mejor aproximacion R(n,m):

R(n,m) = (+exacta) +1

4m − 1[(+exacta)− (-exacta)] .

Teorema 5.2. Si f(x) es una funcion continua en [a, b] entonces la su-cesion que se obtiene en cada columna del metodo de Romberg converge a∫ ba f(x)dx.

Ejemplo 5.10. Aplica el metodo de Romberg para obtener una aproxi-macion de

∫ 20 ex

2dx.

118 Capıtulo 5

R(0, 0) es la regla del trapecio con 20 = 1 intervalo:

R(0, 0) =2− 0

2(e0 + e2

2) ≈ 55,5981

R(1, 0) es la regla del trapecio compuesta con 21 intervalos:

R(1, 0) =1

2

[e0 + 2e1

2+ e2

2]

= 30,5174

R(2, 0) es la regla del trapecio compuesta con 22 intervalos:

R(2, 0) =2/4

2

[e0 + 2e0,5

2+ 2e1

2+ 2e1,5

2+ e2

2]

= 20,6446

y de modo similar los terminos R(n, 0) de la primera columna. El resto determinos se calculan segun la formula 5.8 La mejor aproximacion obtenidaserıa

∫ 20 ex

2dx ≈ 16,4529.

m = 0 m = 1 m = 2 m = 3 m = 4

n = 0 55,5981

n = 1 30,5174 22,1572

n = 2 20,6446 17,3537 17,0334

n = 3 17,5651 16,5386 16,4843 16,4756

n = 4 16,7354 16,4588 16,4535 16,4530 16,4529

Ejemplo 5.11. Aplica el metodo de Romberg para obtener una aproxi-macion de

∫ 20 x2dx.

Obtenemos la tabla siguiente en donde la mejor aproximacion obtenida serıa∫ 20 x2dx ≈ 2,6667.

m = 0 m = 1 m = 2 m = 3 m = 4

n = 0 4

n = 1 3 2,6667

n = 2 2,75 2,6667 2,6667

n = 3 2,6875 2,6667 2,6667 2,6667

n = 4 2,6719 2,6667 2,6667 2,6667 2,6667


5.5. Cuadratura adaptativa

A grandes rasgos, el metodo consiste en calcular un valor aproximadode una integral

∫ ba f(x)dx dividiendo el intervalo [a, b] en intervalos disjun-

tos [xi, xi+1] con las siguientes condiciones:

Se fija un parametro llamado tolerancia T.

Se obtienen dos estimaciones de∫ xi+1

xif(x)dx mediante una formula de

cuadratura.

Se fija la tolerancia en el intervalo [xi, xi+1]:

T[ xi, xi+1] =T

(b− a)(xi+1 − xi).

Si la diferencia entre las dos estimaciones de∫ xi+1

xif(x)dx es menor que

T[xi,xi+1] se estima definitivamente∫ xi+1

xif(x)dx por extrapolacion. En

caso contrario se divide el intervalo en dos subintervalos disjuntos.

Se estima∫ ba f(x)dx:∫ b

af(x)dx =

n∑i=0

∫ xi+1

xi

f(x)dx

Ejemplo 5.12. Fijado un valor de tolerancia T = 0,02, calcula∫ 10,2

1x2dx

usando la regla de Simpson.

Paso 1. Aplicamos la regla de Simpson para calcular la integral:

h1 =1− 0,2

2= 0,4⇒ S1[0,2, 1] = 4,94814815.

Paso 2. Aplicamos la regla de Simpson compuesta en dos subintervalos de[0,2, 1], en concreto en [0,2, 0,6] y en [0,6, 1]:

h2 =0,6− 0,2

2= 0,2 ⇒ S2[0,2, 0,6] = 3,51851852.

h2 =1− 0,6

2= 0,2 ⇒ S2[0,6, 1] = 0,66851852.

120 Capıtulo 5

Paso 3. Hasta ahora tenemos dos aproximaciones a∫ 10,2

1x2dx, a saber:

S1[0,2, 1] y S2[0,2, 0,6] + S2[0,6, 1]. Calculo la diferencia entre ambas estima-ciones:

S1[0,2, 1]− (S2[0,2, 0,6] + S2[0,6, 1]) = 0,7611111 > T = 0,02.

Si la diferencia hubiera sido menor que T, el proceso finalizarıa efectuandouna extrapolacion del tipo:

(+exacta) +1

15[(+exacta)− (-exacta)]

entendiendo que la primera estimacion es menos exacta que la segunda. Puestoque la diferencia es mayor que T, el proceso continua fijando la tolerancia deforma proporcional al tamano de cada subintervalo y de modo que sume latolerancia total T = 0,02. Con estas condiciones la tolerancia sera 0.01 en[0,2, 0,6] y 0.01 en [0,6, 1].

Paso 4. Estudiamos el intervalo [0,6, 1]. Tenemos la estimacion S2[0,6, 1]de∫ 10,6

1x2dx y obtenemos otra sumando los resultados obtenidos al aplicar la

regla de Simpson en los intervalos [0,6, 0,8] (S3[0,6, 0,8] = 0,41678477) y en[0,8, 1] (S3[0,6, 1] = 0,25002572). Comparamos ambas estimaciones:

S2[0,6, 1]− (S3[0,6, 0,8] + S3[0,8, 1]) =

= 0,66851852− 0,66681049 = 0,001708 < T[0,6,1] = 0,01

Como la diferencia es menor que la tolerancia, extrapolando obtenemos laestimacion definitiva de

∫ 10,6

1x2dx :∫ 1

0,6

1

x2dx ≈ (+exacta) +

1

15[(+exacta)− (-exacta)]

0,66681049 +1

15(0,66681049− 0,66851852) = 0.66669662

Paso 5. Actuamos ahora sobre el intervalo [0,2, 0,6] donde la toleran-cia tambien es 0,01. Contamos con una estimacion S2[0,2, 0,6] de la integralen el intervalo y aplicamos la regla de Simpson en [0,2, 0,4] (S3[0,2, 0,4] =2,52314815) y en [0,4, 0,6] (S3[0,4, 0,6] = 0,83435926) para obtener una nueva(S3[0,2, 0,4] + S3[0,4, 0,6]). Comparamos ambas estimaciones:

S2[0,2, 0,6]− (S3[0,2, 0,4] + S3[0,4, 0,6]) = 0,161111 > T[0,2,0,6] = 0,01.


Puesto que la diferencia es mayor que la tolerancia seguimos el proceso ydividimos el intervalo [0,2, 0,6] en los intervalos [0,2, 0,4] y [0,4, 0,6]. Fijamosla tolerancia de cada intervalo de forma proporcional al tamano del intervalocorrespondiente y de modo que la suma sea igual a la tolerancia en el intervalo[0.2,0.6] (T[0,2,0,6] = 0,01). Con estas condiciones T[0,2,0,4] = 0,005 y T[0,4,0,6] =0,005.

Paso 6. Estudiamos el intervalo [0,4, 0,6]. Tenemos la estimacion S3[0,4, 0,6]de∫ 0,60,4


regla de Simpson en los intervalos [0,5, 0,6] (S4[0,4, 0,5] = 0,50005144) y en[0,5, 0,6] (S4[0,5, 0,6] = 0,33334864). Comparamos ambas estimaciones:

S3[0,4, 0,6]− (S4[0,4, 0,5] + S4[0,5, 0,6]) = 0,000859 < T[0,4,0,6] = 0,005

Como la diferencia es menor que la tolerancia, extrapolando obtenemos laestimacion definitiva de

∫ 0,60,4

1x2dx :∫ 0,6

0,4

1


1

15[(+exacta)− (-exacta)] = 0.8333428.



regla de Simpson en los intervalos [0,2, 0,3] y [0,3, 0,4]. Comparamos ambasestimaciones:

S3[0,2, 0,4]− (S4[0,2, 0,3] + S4[0,3, 0,4]) > T[0,2,0,4] = 0,005.

Puesto que la diferencia es mayor que la tolerancia seguimos el proceso y di-vidimos el intervalo [0,2, 0,4] en los intervalos [0,2, 0,3] y [0,3, 0,4]. Fijamos latolerancia de cada intervalo de forma proporcional al tamano del intervalocorrespondiente y de modo que la suma sea igual a la tolerancia en el inter-valo [0.2,0.4] (T[0,2,0,4] = 0,005). Con estas condiciones T[0,2,0,3] = 0,0025 yT[0,3,0,4] = 0,0025.




S4[0,3, 0,4]− (S5[0,3, 0,35] + S5[0,35, 0,4]) = 0,0002220 < T[0,3,0,4] = 0,0025.

122 Capıtulo 5

Como la diferencia es menor que la tolerancia extrapolando obtenemos laestimacion definitiva de

∫ 0,40,3

1x2dx :∫ 0,4

0,3

1


1

15[(+exacta)− (-exacta)] = 0.83333492.




S4[0,2, 0,3]− (S5[0,2, 0,25] + S5[0,25, 0,3]) < T[0,2,0,3] = 0,0025.

Como la diferencia es menor que la tolerancia extrapolando obtenemos laestimacion definitiva de

∫ 0,30,2

1x2dx :∫ 0,3

0,2

1


1

15[(+exacta)− (-exacta)] = 1.666686.

Por lo tanto, la estimacion mediante este metodo de∫ 10,2

1x2dx serıa:∫ 1

0,2

1

x2dx =

∫ 0,3

0,2

1

x2dx +

∫ 0,4

0,3

1

x2dx +

∫ 0,6

0,4

1

x2dx +

∫ 1

0,6

1

x2dx

= 1,666686 + 0,83333492 + 0,8333428 + 0,66669662 = 4,00005957.



http://sage.unex.es

1. MUI_Integracion_numerica

Tab

ladiferen

ciacione

integracion

num

erica

Derivacion

f′(x

)=

f(x

+h)−

f(x

)h

−

h2 f′′(ξ)

f′(x

)=

f(x

+h)−

f(x

−h)

2h

−

h2

6f′′′(ξ)

f′′(x

)=

f(x

+h)−

2f(x

)+f(x

−h)

h2

−

h2

12 f

4)(ξ)

f′(x

k )=

∑

ni=0f(x

i )L′i (x

k )1

(n+

1)! f

n+

1)(ξ

xk )Π

nj=0,j6=

k (xk−

xj )

R(0,h

)=

f(x

+h)−

f(x

−h)

2h

R(n

,h)

=R

(n−

1,h/2)

+1

4n−

1 [R(n

−1,h

/2)−

R(n

−1,h

)]

Integracion

∑

ni=0 (

∫

bali (x

)dx)f

(xi )

h2 [f(x

0 )+

f(x

1 )]−

112 h

3f′′(ξ)

h3 [f(x

0 )+

4f(x

1 )+

f(x

2 )]−

190 h

5f4)(ξ)

3h8[f

(x0 )

+3f

(x1 )

+3f

(x2 )

+f(x

3 )]−

380 h

5f4)(ξ)

h2 [f(x

0 )+

2f(x

1 )+

2f(x

2 )+···+

2f(x

n−

1 )+

f(x

n )]−

(b−a)

12

h2f

′′(ξ)h3 [f

(x0 )

+4f

(x1 )

+2f

(x2 )

+4f

(x3 )

+···+

4f(x

n−

1 )+

f(x

n )]−

(b−a)

180

h4f

4)(ξ)

f(−

√

13 )+

f(√

13 )

R(n

,m)

=R

(n,m

−1)

+1

4m−

1 [R(n

,m−

1)−

R(n

−1,m

−1)]

Si+

1 [a,b]+

115

[Si+

1 [a,b]

−S

i [a,b]

]

1

Capıtulo 6

Resolucion de ecuacionesdiferenciales ordinarias deprimer orden

6.1. Existencia y unicidad de soluciones

Una ecuacion diferencial ordinaria (EDO) es una expresion de la forma

F (x, y(x), y′(x), y′′(x), . . . ) = 0,

es decir, una ecuacion que contiene una variable independiente x, otra depen-diente de la primera y(x) y derivadas de distinto orden de la segunda respectode la primera (y′(x), y′′(x), . . . ).

El orden de una EDO es el de la mayor derivada que aparece en la ecuacion.El grado de una EDO es el mayor exponente de la derivada de mayor ordenque aparece en la ecuacion. Por ejemplo,

ecuacion orden grado

y′ = y tan(x + 3) 1 1

(y′)2 = y tan(x + 3) 1 2

y′′ = y tan(x + 3) 2 1

Por supuesto, la notacion de las variables independiente y dependiente nosiempre ha de ser x e y. En general, el contexto de la ecuacion determina quien

125

126 Capıtulo 6

es una y quien es otra. Por ejemplo, en la ecuacion x′(t) = t x(t), la variabledependiente es x y la independiente es t.

La solucion general de una EDO es una familia de funciones que verificanla EDO. En general, la solucion general de una EDO de orden k es una familiade funciones que depende de una parametro k.

Ejemplo 6.1. La solucion general de y′ = sen(x) es y = − cos(x) + k conk ∈ R. La solucion general de y′′ = −y es y = a sen(x) + b cos(x), con a ∈ R yb ∈ R.

Una solucion particular de una EDO es un elemento particular de la solu-cion general. Una solucion singular es una funcion que verifica la EDO peroque no esta recogida en la solucion general.

Ejemplo 6.2. La ecuacion y′ = y1/3 tiene como solucion general y =(23x)

32 +C con C ∈ R. Una solucion particular es y = (23x)

32 +4 y una solucion

singular es y = 0.

Si una EDO de primer orden F (x, y(x), y′(x)) = 0 se puede expresar de laforma y′ = f(x, y(x)), entonces

-10-5

05

10-10

-5

0

5

10

-5-2.5

02.55

-10-5

05

Figura 6.1: La funcion f(x, y) = yx

Ecuaciones diferenciales ordinarias 127

Teorema 6.1. Si f(x, y) y ∂f(x,y)∂y son continuas en un rectangulo D del

plano R2 y (x0, y0) ∈ D, entonces existe una unica solucion de la ecuaciony′ = f(x, y) tal que y(x0) = y0.

Ejemplo 6.3. Consideremos la ecuacion y′ = −yx , es decir, sea f(x, y) =

−yx . La grafica de dicha funcion aparece en la figura 6.1. f(x, y) es continua en

R2 − x = 0. ∂f(x,y)∂y = −1

x es continua tambien en R2 − x = 0.Por el teorema anterior, si existe un rectangulo D ⊂ R2 − x = 0 y

(x0, y0) ∈ D, entonces el problema y′ = f(x, y) tal que y(x0) = y0 tienesolucion unica. Por ejemplo, el problema y′ = −y

x con y(2) = 1, tiene solucionunica porque (2, 1) estan dentro de un rectangulo que a su vez estan contenidoen R2 − x = 0. Dicha solucion es y = 2

x . Sin embargo, si la condicion iniciales y(0) = 2 el problema no tiene solucion.

-2 -1 0 1 2

-2

-1

0

1

2

Figura 6.2: campo de pendientes de y′ = −2x− y

Si la ecuacion se plantea de la forma y′ = f(x, y(x)), podemos dibujar loque se conoce como campo de pendientes, es decir, una grafica en la que en cadapunto del plano (x, y) se traza un pequeno segmento con pendiente f(x, y(x)),que es la pendiente que le corresponde dado que y′ = f(x, y(x)). Por ejemplo,

128 Capıtulo 6

el campo de pendiente de la ecuacion y′ = −2x−y se presenta en la Figura 6.2.Con SAGE, podemos dibujar el campo de pendiente de la ecuacion y′ = f(x, y)con el comando plot_slope_field(f(x,y), (x, , ), (y, , )).

En general, la busqueda de soluciones exactas de una ecuacion diferenciales un problema difıcil. Por ello se usan metodos numericos para encontrarsoluciones aproximadas. Presentamos algunos de ellos en la seccion 6.2. Noobstante, hay algunos tipos de ecuaciones que sı se pueden resolver de formaexacta. Algunas de estas ultimas ecuaciones se presentan en la seccion 6.3

6.2. Resolucion numerica de EDO de primer orden

En esta seccion se presentan algunos metodos para resolver EDO de primerorden con valor inicial, es decir, problemas con la siguiente expresion general:

y′ = f(x, y(x)) y(x0) = y0 (6.1)

Por ejemplo, un problema de valor inicial es el siguiente:

y′ = y tan(x + 3) y(−3) = 1.

La solucion numerica de una EDO sera una tabla de valores aproximadosde la funcion y(x) en un conjunto de puntos xi. Usando esa tabla podemosconstruir una aproximacion de la funcion mediante, por ejemplo, los metodosde interpolacion. En lo que sigue, trabajaremos con puntos xi igualmente es-paciados, es decir, separados cada uno del siguiente por un tamano de paso h,con lo que xi = x0 + ih. y(xi) denotara el valor exacto de la funcion y(x) enel punto xi; yi sera un valor aproximado de y(xi).

6.2.1. Metodo de la serie de Taylor

Para resolver el problema (6.1), el metodo de Taylor de orden k proponesustituir la funcion desconocida y(x) por su polinomio de Taylor de orden ken un entorno de cada punto xi, es decir

y(xi+1) = y(xi + h) ≈ y(xi) + y′(xi)h + y′′(xi)h2

2+ · · ·+ yk)(xi) +

hk

k!

Los valores de la funcion y(x) y de sus derivadas en xi son desconocidos engeneral, por lo que hay que sustituirlos por aproximaciones:


y(xi) ≈ yi

y′(xi) ≈ f(xi, yi)

y′′(xi) ≈ d f(xi, yi)

dx. . .

yk)(xi) ≈ dk-1 f(xi, yi)

dxk−1

De este modo obtenemos la expresion del metodo:

yi+1 = yi + f(xi, yi)h +d f(xi, yi)

dx

h2

2+ · · ·+ dk-1 f(xi, yi)

dxk−1

hk

k!

A la hora de derivar f(x, y) respecto de x hay que tener en cuenta que es unafuncion de dos variables y que y es funcion de x. De modo que, denotando fx yfy a las derivadas parciales de f(x, y) respecto x e y, tendrıamos, por ejemploque

d f(x, y)

dx= fx(x, y) + fy(x, y)y′(x, y) = fx(x, y) + fy(x, y)f(x, y).

Elegido, el tamano del paso h, el metodo de Taylor paso a paso serıa elsiguiente:

Paso 1. Partimos del punto x0, obtenemos una aproximacion de y(x1)

y1 = y0 + f(x0, y0)h +d f(x0, y0)

dx

h2

2+ · · ·+ dk-1 f(x0, y0)

dxk−1

hk

k!

Paso 2. Partimos de x1 = x0 + h y calculamos una aproximacion y(x2):

y2 = y1 + f(x1, y1)h +d f(x1, y1)

dx

h2

2+ · · ·+ dk-1 f(x1, y1)

dxk−1

hk

k!

Paso i+1. Partimos de xi = xi−1 + h y calculamos una aproximacion dey(xi+1):

yi+1 = yi + f(xi, yi)h +d f(xi, yi)

dx

h2

2+ · · ·+ dk-1 f(xi, yi)

dxk−1

hk

k!

130 Capıtulo 6

Ejemplo 6.4. Proporciona tres valores aproximados de la funcion y(x)que sea solucion de y′ = −2x−y con y(0) = −1 mediante el metodo de Taylorcon k = 4 y tamano de paso h = 0,1.

Solucion. Calculamos el desarrollo de Taylor hasta grado 4 de la funciony(x) en un punto cualquiera xi:

y(xi + h) ≈ y(xi) + y′(xi)h +y′′(xi)

2h2 +

y′′′(xi)

3!h3 +

y4)(xi)

4!h4

Denotando yi al valor aproximado de y(xi) y dado que

y′(x) = f(x, y)

y′′(x) = [y′(x)]′ = [−2x− y(x)]′ = −2− y′(x) = −2− f(x, y)

y′′′(x) = [y′′(x)]′ = −y′′(x) = 2 + y′(x) = 2 + f(x, y)

y4)(x) = [y′′′(x)]′ = [−y′′(x)]′ = −y′′′(x) = −2− y′(x) = −2− f(x, y)

definimos el metodo

yi+1 = yi + f(xi, yi)h +−2− f(xi, yi)

2h2 +

2 + f(xi, yi)

3!h3 +

−2− f(xi, yi)

4!h4

Paso 1. Como y0 = −1 y f(x0, y0) = f(0,−1) = 1, obtenemos el valoraproximado

y1 = −1 + 0,1− 1,5 · 0,12 + 0,5 · 0,13 − 0,125 · 0,14 = −0,9645125.

Paso 2. Como y1 = −0,9145125 y f(x1, y1) = f(0,1,−0,9145125) = 0,7145125,obtenemos el valor aproximado

y2 = −0,8561927.

Paso 3. Como y2 = −0,8561927 y f(x2, y2) = f(0,2,−0,8561927) = 0,4561927,obtenemos el valor aproximado

y3 = −0,8224553


Metodo de Euler

Es el metodo anterior usando el polinomio de Taylor de orden k = 1.

yi+1 = yi + f(xi, yi)h

Ejemplo 6.5. Resuelve la ecuacion y′ = −2x−y con y(0) = −1 medianteel metodo de Euler.

Solucion. Con estos datos se tiene que x0 = 0, y0 = −1, f(x0, y0) = 1 con loque

y1 = y0 + f(x0, y0)h = −1 + 1 · 0,1 = −0,9

Puesto que f(x1, y1) = 0,7, se deduce que

y2 = y1 + f(x1, y1)h = −0,9 + 0,7 · 0,1 = −0,83


y3 = y2 + f(x2, y2)h = −0,83 + 0,43 · 0,1 = −0,787


y4 = y3 + f(x3, y3)h = −0,787 + 0,187 · 0,1 = −0,7683

Los comandos de SAGE:

x,y=PolynomialRing(QQ,2,xy”).gens()”eulers_method(f(x,y),x_0,y_0,h,x_1)

aplican el metodo de Euler para la resolucion de la ecuacion y′ = f(x, y) convalor inicial y(x0) = y0, con tamano de paso h y siendo el ultimo valor x1.

Ejercicio 6.1. Aplicando el metodo de Taylor de grado 2 con dos pasos,calcula valores aproximados de y(0,1) y y(0,2), sabiendo que

y′ + 2y = x2 + 2x y(0) = 1

Ejercicio 6.2. Resuelve la ecuacion y′ = y2 − t2 + 1, y(0) = 0,5, siendo0 ≤ t ≤ 1,2, mediante el metodo de Taylor de grado 2 con 5 pasos.

132 Capıtulo 6

6.2.2. Errores

Al resolver numericamente una ecuacion diferencial aparecen varios tiposde errores:

1. error de truncamiento local

2. error de truncamiento global

3. error de redondeo local

4. error de redondeo global

5. error total

El error de truncamiento local es el que aparece en un paso cuando re-emplazamos un proceso infinito por uno finito. Por ejemplo, en el metodo deTaylor, usando el desarrollo hasta grado k, sustituimos la serie de Taylor (deinfinitos terminos) por los primeros k-terminos. En este caso, para un ciertoci ∈ R, se tiene que

y(xi + h) = y(xi) + y′(xi)h +y′′(xi)

2!h2 + · · ·+ yk)(xi)

k!hk +

yk+1)(ci)

(k + 1)!hk+1

y dado que en el algoritmo usamos la aproximacion

y(xi + h) ≈ y(xi) + y′(xi)h +y′′(xi)

2!h2 + · · ·+ yk)(xi)

k!hk

se deduce que el error local en el paso i-esimo es proporcional a hk+1, digamosCih

k+1. Un error de esa magnitud se dice que es O(hk+1).La acumulacion de todos los errores de truncamiento local genera el error

de truncamiento global. Si al aplicar un metodo el error de truncamientolocal es O(hk+1), el tamano de paso es h, el punto inicial es x0 y el punto finales xn, entonces el numero de pasos sera n = xn−x0

h y puesto que

C1hk+1 + C2h

k+1 + · · ·+ Cnhk+1 = hk+1(C1 + C2 + · · ·+ Cn) ≤

≤ hk+1nmaxC1, C2, . . . , Cn = Chk

siendo C = (xn−x0) maxC1 +C2 + · · ·+Cn, el error de truncamiento globalsera O(hk). Si un metodo tiene error de truncamiento global O(hk) decimos


que el procedimiento numerico es de orden k. El metodo de Taylor usando eldesarrollo hasta grado k es un metodo de orden k.

El error de redondeo local es el que en cada paso es provocado por lalimitada precision de nuestras maquinas de computo, como ya se comento eltema inicial del curso. El error de redondeo global es la acumulacion de loserrores de redondeo local en los pasos anteriores. El error total es la sumade los errores de truncamiento global y redondeo global.

6.2.3. Metodos de Runge-Kutta

Los metodos de Taylor tienen el inconveniente de la dificultad de calculode las sucesivas derivadas de la funcion y(x). En los metodos de Runge-Kuttase elimina esa dificultad. Las ideas que dan origen a este conjunto de metodosse resumen en los siguientes pasos:

Paso 1. Consideramos la igualdad

y(xi+1)− y(xi) =

∫ xi+1

xi

y′(x)dx =

∫ xi+1

xi

f(x, y(x))dx.

y despejamos y(xi+1)

y(xi+1) = y(xi) +

∫ xi+1

xi

y′(x)dx = y(xi) +

∫ xi+1

xi

f(x, y(x))dx.

Paso 2. Realizamos el cambio de variable x = xi+θh en la integral anteriory obtenemos

y(xi+1) = y(xi) + h

∫ 1

0f(xi + θh, y(xi + θh))dθ. (6.2)

Paso 3. Consideramos el conjunto de puntos 0 = θ0 ≤ θ1 ≤ . . . θm ≤ 1 delintervalo de integracion [0, 1] y estimamos la integral anterior mediante unaformula de cuadratura del siguiente tipo∫ 1

0F (θ)dθ ≈ A0F (θ0) + A1F (θ1) + · · ·+ AmF (θm)

con lo que, denotando xiθj = xi + θjh, en la ecuacion (6.2) obtenemos

y(xi+1) ≈ y(xi) + h[A0 f

(xiθ0 , y(xiθ0)

)+

A1 f(xiθ1 , y(xiθ1)

)+ · · ·+ Am f

(xiθm , y(xiθm)

)]. (6.3)

134 Capıtulo 6

La formula de cuadratura elegida puede variar, pero se le ha de exigir queal menos sea exacta para las funciones constantes, lo que se traduce en quelos coeficientes han de cumplir la condicion

A0 + A1 + · · ·+ Am = 1

El metodo de Runge-Kutta sera la formula que resulta de sustituir en la ecua-cion (6.3) los valores y(xi), y(xiθ0), y(xiθ1), . . . , y(xiθm) por sus valores apro-ximados. Si bien dispondremos en cada paso del valor yi que aproxima a y(xi),no ocurre lo mismo con las aproximaciones de y(xiθ0), y(xiθ1), . . . , y(xiθm).

Paso 4. En esta etapa conseguiremos obtener un algoritmo para las aproxi-maciones aludidas en el apartado anterior. De forma similar a como se obtuvola ecuacion (6.2), llegamos a la igualdad

y(xiθj ) = y(xi) + h

∫ θj

θ0

f(xi + θh, y(xi + θh))dθ

y aplicando a la integral anterior una formula de cuadratura para cada j deltipo siguiente∫ θj

θ0

F (θ)dθ ≈ Bj0F (θ0) + Bj

1F (θ1) + · · ·+ Bjj−1F (θj−1) (6.4)

se obtienen las aproximaciones buscadas

y(xiθj ) ≈ y(xi) + h

j−1∑k=0

Bjkf(xiθk , y(xiθk)).

Las formulas de cuadraturas usadas en (6.4) han de ser exactas al menos paralas funciones constantes, por lo tanto los coeficientes deben verificar

Bj0 + Bj

1 + · · ·+ Bjj−1 = θj j = 1, . . . ,m.

Paso 5. Si denotamos yiθj a una estimacion de y(xiθj ), un metodo de Runge-Kutta quedarıan ası:

yi+1 = yi + h[A0f(xiθ0 , yiθ0) + A1f(xiθ1 , yiθ1) + · · ·+ Amf(xiθm , yiθm)]

(6.5)


siendo

yiθ0 = yiyiθ1 = yi + h

[B1

0f(xiθ0 , yiθ0)]

yiθ2 = yi + h[B2

0f(xiθ0 , yiθ0) + B21f(xiθ1 , yiθ1)

]. . .yiθm = yi + h

[Bm

0 f(xiθ0 , yiθ0) + Bm1 f(xiθ1 , yiθ1) + . . .· · ·+ Bm

m−1f(xiθm−1 , yiθm−1)]

Previamente deben fijarse los valores

0 = θ0 ≤ θ1 ≤ θ2 ≤ . . . θm ≤ 1

ası como los siguientes parametros que han de cumplir con las condiciones delas sumas correspondientes que aparecen en el cuadro

A0 A1 A2 . . . Am∑m

r=0Ar = 1

B10 B2

0 Bm0

B21 . . . Bm

1

. . .Bm

m−1

B10 = θ1 B2

0 + B21 = θ2 . . .

∑m−1k=0 Bm

k = θm

Un metodo como el anterior se conoce con el nombre de metodo de Runge-Kutta de m+1 etapas, debido a que la formula (6.5) requiere m+1 evaluacionesde la funcion f(x, y).

Para m = 0 (una etapa), el metodo de Runge-Kutta que resulta es elmetodo de Euler.

Para m = 1 (dos etapas), podemos obtener, entre otras posibilidades, elmetodo de Heun:

θ0 = 0, θ1 = 1

A0 = 1/2 A1 = 1/2 A0 + A1 = 1

B10 = 1

B10 = θ1 = 1

yi+1 = yi + h[1

2f(xi0, yi0) +

1

2f(xi1, yi1)

]siendo

yi0 = yiyi1 = yi + hf(xi, yi)

136 Capıtulo 6

con lo que

yi+1 = yi + h[12f(xi, yi) +

1

2f(xi+1, yi + hf(xi, yi))

](6.6)

El llamado metodo de Euler modificado sale con la eleccion

θ0 = 0, θ1 = 1/2

A0 = 0 A1 = 1 A0 + A1 = 1

B10 = 1/2

B10 = θ1 = 1/2

yi+1 = yi + h[f(xiθ1, yiθ1)]

siendoyi0 = yi

yi1/2 = yi + h12f(xi, yi)

con lo que el metodo queda

yi+1 = yi + hf(xi +

h

2, yi +

h

2f(xi, yi)

)(6.7)

Otra eleccion posible es

θ0 = 0, θ1 = 2/3

A0 = 1/4 A1 = 3/4 A0 + A1 = 1

B10 = 2/3

B10 = θ1 = 2/3

yi+1 = yi + h[1

4f(xiθ0 , yiθ0) +

3

4f(xiθ1 , yiθ1)

]siendo

yiθ0 = yiyiθ1 = yi + h

[23f(xiθ0 , yiθ0)

]con lo que el metodo queda ası


yi+1 = yi + h[1

4f(xi, yi) +

3

4f(xi +

2

3h, yi + h

2

3f(xi, yi)

]El metodo clasico de Runge-Kutta es de m = 3 (cuatro etapas) definido

por los parametros

0 = θ0 θ1 = 1/2 θ2 = 1/2 θ3 = 1

A0 = 1/6 A1 = 1/3 A2 = 1/3 A3 = 1/6∑3

r=0Ar = 1

B10 = 1/2 B2

0 = 0 B30 = 0

B21 = 1/2 B3

1 = 0B3

2 = 1

B10 = 1/2 B2

0 + B21 = 1/2

∑2k=0B

3k = 1

que da lugar a

yi+1 = yi +h

6(K1 + 2K2 + 2K3 + K4) (6.8)

donde

K1 = f(xi, yi) K2 = f(xi + h2 , yi + h

2K1)

K3 = f(xi + h2 , yi + h

2K2) K4 = f(xi + h, yi + hK3)

Ejemplo 6.6. Aplica tres pasos del metodo clasico de Runge-Kutta deorden 4 con h = 1/128 para resolver la ecuacion

y′ =yx− y2

x2y(1) = 2.

Solucion. Paso 1. Partimos de x0 = 1, y0 = 2. Para aplicar la ecuacion (6.8)previamente calculamos

K1 = f(x0, y0) = f(1, 2) = −2 K2 = f(x0 + h2 , y0 + h

2K1) = −1,95355

K3 = f(x0 + h2 , y0 + h

2K2) = −1,95409K4 = f(x0 + h, y0 + hK3) = −1,90898.

138 Capıtulo 6

Con estos datos obtenemos

y1 = y0 +h

6(K1 + 2K2 + 2K3 + K4) = 1,98473,

Paso 2. Partimos de x1 = 1 + 1128 , y1 = 1,98473. Calculamos

K1 = f(x1, y1) = −1,90899 K2 = f(x1 + h2 , y1 + h

2K1) = −1,86520

K3 = f(x1 + h2 , y1 + h

2K2) = −1,86570K4 = f(x1 + h, y1 + hK3) = −1,82316.


y2 = y1 +h

6(K1 + 2K2 + 2K3 + K4) = 1,97016.

Paso 3. Partimos de x2 = 1 + 2128 , y2 = 1,97016. Calculamos

K1 = f(x2, y2) = −1,82311 K2 = f(x2 + h2 , y2 + h

2K1) = −1,78186

K3 = f(x2 + h2 , y2 + h

2K2) = −1,78231K4 = f(x2 + h, y2 + hK3) = −1,74215.


y3 = y2 +h

6(K1 + 2K2 + 2K3 + K4) = 1,96088.

Ejercicio 6.3. Aplica el metodo clasico de Runge-Kutta para estimar 10valores aproximados de la solucion de la ecuacion siguiente en el intervalo [0, 1]

y′ + 2y = x2 + 2x y(0) = 1

Ejercicio 6.4. A partir de la identidad

y(xi+1)− y(xi) =

∫ xi+1

xi

y′(x)dx

deduce el metodo de Runge-Kutta de orden 2

yi+1 = yi + h[12f(xi, yi) +

1


]El metodo de Runge-Kutta de orden k = 4 esta implementado en SAGE.

Para aplicarlo usamosdesolve_rk4(f(x,y),y,ics=[x_0,y_0],end_points= ,step=h),que recibe la funcion f(x, y) , la incognita (y), la condicion inicial ([x0, y0]),

el punto final y el paso h, y aplica el metodo.


Errores

Los metodos de Runge-Kutta tienen tambien su error de truncamientoglobal, es decir, su orden. Existen metodos de Runge-Kutta de orden k conk menor o igual que el numero de etapas m + 1. Se da la igualdad parak = 1, 2, 3, 4. Es decir, los metodos de Runge-Kutta de 1,2,3,4 etapas tienenorden 1,2,3,4, respectivamente. Para conseguir orden k = 5 es necesario 6etapas; para k = 6, 7 etapas; en general para k ≥ 7 se necesitan un numerode etapas mayor o igual que k + 2.

6.2.4. Metodos multipaso

Tanto los metodos de Taylor como los de Runge-Kutta son metodos de unpaso, porque para calcular el valor aproximado yi+1 solo se usa, ademas de xiy de h, el valor aproximado anterior yi. Los metodos multipaso se caracterizanpor obtener yi+1 a partir de los valores h, xi, yi, xi−1, yi−1, . . . , xi−k, yi−k.

Partimos siempre de la ecuacion

y′ = f(x, y) y(x0) = y0.

Un primer ejemplo sencillo de un metodo multipaso se puede obtener dela siguiente forma. Sabemos que

y′(x) ≈ y(x + h)− y(x− h)

2h

con lo que

y′(xi) ≈y(xi+1)− y(xi−1)

2h

y se deduce que

y(xi+1) ≈ y(xi−1) + 2 h y′(xi),

lo que sugiere el metodo

yi+1 = yi−1 + 2 h f(xi, yi).

Para poder obtener el punto yi+1 son necesarios los dos valores anterioresyi, yi−1. Es un metodo de dos pasos. Para empezar a aplicarlo tenemos quedisponer de y0 e y1. Este ultimo puede calcularse con un metodo de un paso.

140 Capıtulo 6

El que sigue es un proceso para obtener un conjunto de metodos multipasos.Partimos de la relacion

y(xi+1)− y(xi) =

∫ xi+1

xi

y′(x)dx =

∫ xi+1

xi

f(x, y(x))dx,

y aplicamos el cambio de variable x = xi + θh

y(xi+1)− y(xi) = h

∫ 1

0f(xi + θh, y(xi + θh))dθ.

Aproximamos la integral anterior por una formula de cuadratura∫ 1

0F (θ)dθ = A0F (θ0) + A1F (θ1) + · · ·+ AkF (θk) (6.9)

y nos queda

y(xi+1)− y(xi) ≈ h[A0f(xi + θ0h, y(xi + θh))+

+ A1f(xi + θ1h, y(xi + θh))+

+ · · ·+ Akf(xi + θkh, y(xi + θkh))]

(6.10)

lo cual sugiere un metodo que resultara de sustituir los valores de la funcionpor valores aproximados.

Ejemplo 6.7. Deduzcamos el metodo de orden 2

yi+1 = yi + h[32f(xi, yi)−

1

2f(xi−1, yi−1)

](6.11)

Solucion. Si en la formula de cuadratura (6.9) tomamos los puntos

θ0 = 0 θ1 = −1

y le imponemos la condicion de que sea exacta para los polinomios de gradomenor o igual que 1, obtendremos que los coeficientes A0 y A1 deben verificarlas ecuaciones

A0 + A1 = 1 A1 = −1

2


y la solucion del sistema anterior es A0 = 32 y A1 = −1

2 . Sustituyendo en (6.10)nos queda

y(xi+1)− y(xi) ≈ h[32f(xi, y(xi))−

1

2f(xi−1, y(xi−1)

]lo cual sugiere el metodo

yi+1 = yi + h[32f(xi, yi)−

1

2f(xi−1, yi−1

].

Ejemplo 6.8. Aplica el metodo anterior al calculo de la solucion de

y′ =yx− x2

x2y(1) = 2

Solucion.Para aplicar la formula (6.11) necesitamos dos puntos iniciales y0, y1,y solo disponemos del valor y0 = 2. Es habitual calcular los puntos necesariospara empezar mediante un metodo de un paso del mismo orden. Usamos elmetodo de Euler modificado (ver ecuacion (6.7)), que es un metodo de Runge-Kutta de orden 2:

y1 = y0 + hf(x0 +

h

2, y0 +

h

2f(x0, y0)

)= 1,98474.

Ya disponemos de y0 e y1 y podemos calcular y2 aplicando (6.11)

y2 = y1 + h[32f(x1, y1)−

1

2f(x0, y0

]= 1,97018.

Con y1 e y2 podemos calcular y3

y3 = y2 + h[32f(x2, y2)−

1

2f(x1, y1

]= 1,95627.

Y ası sucesivamente.

142 Capıtulo 6

Formulas de Adams-Bashforth

Un metodo multipaso de la forma

yi+1 = yi + A0f(xi, yi) + A1f(xi−1, yi−1) + · · ·+ Ak−1f(xi−k+1, yi−k+1)

es una formula de Adams-Bashforth de k pasos y orden k. Estos metodosse encuadran dentro de los metodos multipaso explıcitos, debido a que ladefinicion de yi+1 es explıcita en funcion de los valores anteriores, a diferenciade los metodos implıcitos que veremos en el siguiente apartado. Estos sonalgunos ejemplos de formulas de Adams-Bashforth:

1. de orden 2:

yi+1 = yi + h[32f(xi, yi)− 1

2f(xi−1, yi−1)]

2. de orden 4:

yi+1 = yi +h

24

[55f(xi, yi)− 59f(xi−1, yi−1)+

+ 37f(xi−2, yi−2)− 9f(xi−3, yi−3)]

3. de orden 5:

yi+1 = yi +h

720

[1901f(xi, yi)− 2774f(xi−1, yi−1)+

+ 2616f(xi−2, yi−2)− 1274f(xi−3, yi−3) + 251f(xn−4, yn−4)]

Todas ellas se pueden obtener siguiendo el proceso descrito en la pagina140. Dado que en dicho proceso se utiliza el cambio de variable x = xi + θ h,si en la formula final (6.10) deseo que aparezca la funcion f(x, y) valorada enlos puntos

xi, xi−1, ..., xi−k+1,

tendre que elegir los valores

θ0 = 0, θ1 = −1, ..., θk−1 = −k + 1

al construir la formula de cuadratura (6.9).Para aplicar una formula de Adams-Bashforth de k pasos necesitaremos

los valores iniciales y0, y1, ..., yk−1, que podran obtenerse mediante un metodode un paso del mismo orden.


Formulas de Adams-Moulton

Un metodo multipaso de la forma

yi+1 = yi + Bf(xi+1, yi+1) + A0f(xi, yi) + A1f(xi−1, yi−1) + . . .

· · ·+ Ak−1f(xi−k+1, yi−k+1) (6.12)

es una formula de Adams-Moulton de k pasos y orden k + 1. Estos metodosse encuadran dentro de los metodos multipaso implıcitos, debido a que enla formula aparece el valor yi+1 definido de forma implıcita. Al igual quelas formulas de Adams-Bashforth, las formulas de Adams-Moulton se puedenobtener siguiendo el proceso descrito en la pagina 140. Dado que en dichoproceso se utiliza el cambio de variable x = xi + θ h, si en la formula final(6.10) deseo que aparezca la funcion f(x, y) valorada en los puntos

xi+1, xi, xi−1, ..., xi−k+1,

tendre que elegir los valores

θ = 1, θ0 = 0, θ1 = −1, ..., θk = −k + 1

al construir la formula de cuadratura (6.9).

Para aplicar una formula de Adams-Moulton de k pasos necesitaremos losvalores iniciales y0, y1, ..., yk−1, que podran obtenerse mediante un metodo deun paso del mismo orden.

Supongamos que disponemos de los valores yi, yi−1, ..., yi−k+1. ¿Comocalcular yi+1 mediante la formula (6.12)? Se puede hacer de dos formas:

1. Metodo predictor-corrector. Se estima un primer valor predictor y∗i+1

mediante una formula de Adams-Bashforth del mismo orden y despues secalcula el valor corrector yi+1 aplicando

yi+1 = yi + Bf(xi+1, y∗i+1) + A0f(xi, yi) + A1f(xi−1, yi−1) + . . .

· · ·+ Ak−1f(xi−k+1, yi−k+1).

2. Metodo iterativo de punto fijo. En la formula (6.12) todos los terminosde la parte derecha de la igualdad son datos excepto yi+1. En este sentido,

144 Capıtulo 6

yi+1 es un punto fijo de la funcion

Φ(t) = yi + Bf(xi+1, t) + A0f(xi, yi) + A1f(xi−1, yi−1) + . . .

· · ·+ Ak−1f(xi−k+1, yi−k+1).

En consecuencia, yi+1 lo podemos obtener como lımite de la sucesion

t0 = y∗i+1, tn = Φ(tn−1) n ≥ 1

siendo y∗i+1 un valor inicial que se recomienda obtener mediante la formulade Adams-Bashforth del mismo orden. En la practica, es suficiente dos o tresiteraciones para conseguir un valor aceptable de yi+1.

Ejemplo 6.9. Deduzcamos la formula de Adams-Moulton de orden dos

yi+1 = yi + h[12f(xi+1, yi+1) +

1

2f(xi, yi)

](6.13)

Solucion. Aplicaremos el proceso descrito en la pagina 140. Si en la formulade cuadratura (6.9) tomamos los puntos

θ0 = 1 θ1 = 0

y le imponemos la condicion de que sea exacta para los polinomios de gradomenor o igual que 1, obtendremos que los coeficientes A0 y A1 deben verificarlas ecuaciones

A0 + A1 = 1 A0 =1

2

y la solucion del sistema anterior es A0 = 12 y A1 = 1

2 . Sustituyendo en (6.10)nos queda

y(xi+1)− y(xi) ≈ h[12f(xi+1, y(xi+1)) +

1

2f(xi, y(xi)

]lo cual sugiere el metodo

yi+1 = yi + h[12f(xi+1, yi+1) +

1

2f(xi, yi

].


Ejemplo 6.10. Aplica el metodo anterior al calculo de la solucion de

y′ =yx− x2

x2y(1) = 2

Solucion. El metodo a aplicar es de orden 2 y 1 paso, por lo que solo necesita-mos el valor y0 para empezar. Para estimar el resto de valores podemos optarpor el metodo predictor-corrector o por el metodo iterativo de punto fijo.

1. Metodo predictor. Calculamos un valor predictor y∗1. Si es posible, debecalcularse con una formula de Adams-Bashforth del mismo orden, en estecaso la formula (6.11) que es de orden dos. Sin embargo, para aplicar (6.11)necesitamos los dos puntos iniciales y0 e y1, con lo que en este primer paso nopodemos hacer uso de ella. Calculamos y∗1 mediante otro metodo, por ejemplo,por Runge-Kutta de orden 2 (ver 6.6):

y∗1 = y0 + h[12f(x0, y0) +

1

2f(x1, y0 + hf(x0, y0))

]= 1,98474

Con este valor, calculamos el valor corrector con Adams-Moulton

y1 = y0 + h[12f(x1, y

∗1) +

1

2f(x0, y0

]= 1,98473.

Para calcular y2, primero calculo el valor predictor y∗2 con Adams-Bashforthde orden 2 (formula (6.11))

y∗2 = y1 + h[32f(x1, y1)−

1

2f(x0, y0)

]= 1,97017

y despues calculo el corrector con Adams-Moulton

y2 = y1 + h[12f(x2, y

∗2) +

1

2f(x1, y1

]= 1,97015.

En el siguiente paso calculo el predictor y∗3 con Adams-Bashforth

y∗3 = y2 + h[32f(x2, y2)−

1

2f(x1, y1)

]= 1,95624

y despues calculo el corrector con Adams-Moulton

y3 = y2 + h[12f(x3, y

∗3) +

1

2f(x2, y2

]= 1,95622.

El proceso anterior se repite para calcular el resto de valores.

146 Capıtulo 6

2. Metodo iterativo de punto fijo. Dada la funcion

Φ1(t) = y0 + h[12f(x1, t) +

1

2f(x0, y0)

]el valor y1 que buscamos verifica

Φ1(y1) = y1.

Es, por tanto, un punto fijo de la ecuacion

Φ1(t) = t.

Para localizarlo generamos una sucesion iterativa comenzando por un valor y∗1obtenido, si es posible, por Adams-Bashforth del mismo orden. En este casono es posible aplicar Adams-Bashforth puesto que se necesitarıa dos puntosiniciales que no tenemos. Calculamos y∗1 mediante Runge-Kutta de orden 2(y∗1 = 1,98474) y generamos la sucesion iterativa:

t0 = y∗1 = 1,98474t1 = Φ1(1,98474) = 1,98473t2 = Φ1(1,98473) = 1,98473

De forma que tomamos y1 = 1,98473.y2 es un punto fijo de

Φ2(t) = y1 + h[12f(x2, t) +

1

2f(x1, y1)

]que localizamos mediante la sucesion iterativa similar a la anterior que co-mienza en y∗2 = 1,97017, que a su vez es obtenido por Adams-Bashforth deorden 2

t0 = y∗2 = 1,97017t1 = Φ2(1,97017) = 1,97015t2 = Φ2(1,97015) = 1,97015

Tomamos y2 = 1,97015. Y repetimos el proceso para el resto de valores.

Ejercicio 6.5. 1. Deduce la formula de Adams-Bashforth de orden 4:

yi+1 = yi +h

24

[55f(xi, yi)− 59f(xi−1, yi−1)+

+ 37f(xi−2, yi−2)− 9f(xi−3, yi−3)]


y usala para resolver la ecuacion

y′ =yx− x2

x2y(1) = 2

con h = 1/128

2. Repite el apartado anterior para la formula de Adams-Moulton de orden4:

yi+1 = yi +h

24

[9f(xi+1, yi+1 + 19f(xi, yi)−

− 5f(xi−1, yi−1) + f(xi−2, yi−2)]

Ejercicio 6.6. Usa la identidad

y(xi+2)− y(xi) =

∫ xi+2

xi

y′(x)dx

y la formula de cuadratura de Simpson 1/3 para obtener la formula del si-guiente metodo multipaso

yi+2 = yi +h

3

[f(xi, yi) + 4f(xi+1, yi+1) + f(xi+2, yi+2)

].

6.3. Resolucion exacta de EDO de primer orden

En esta seccion plantearemos ecuaciones diferenciales de la forma

y′ = f(x, y(x))

donde y(x) sera una funcion que dependa de x, y buscaremos la soluciongeneral de algunas EDO de primer orden de este tipo.

El comando desolve de SAGE resuelve algunas EDO de primer y segundoorden. desolve(diff(y,x) + y - 1, y) resuelve la ecuacion y′(x) + y− 1 =0. El comando desolve(diff(y,x) + y - 1, y, ics=[10,2]) resuelve lamisma ecuacion con la condicion inicial y(10) = 2.

Ejemplo 6.11. : El modelo de crecimiento de poblaciones de Malthaus seresume en la ecuacion x′(t) = Kx(t), donde t es el tiempo, x es el numero deindividuos y K es una constante que depende del tipo de poblacion en estudio.

148 Capıtulo 6

Resuelve la ecuacion diferencial que aparece en el modelo de crecimiento depoblaciones de Malthaus cuando K = 2 y x(0) = 123:

x′(t) = 2x(t), x(0) = 123.

Solucion. Si x′(t) = 2x(t), entonces

x′(t)

x(t)= 2

e integrando a ambos lados de la igualdad respecto de t∫x′(t)

x(t)dt =

∫2dt

obtenemosL(x(t)) = 2t + C.

Para despejar x(t) usamos la funcion exponencial

eL(x(t)) = e2t+C

de modo que x(t) = e2teC . Si lamamos M = eC , la solucion general sera

x(t) = e2tM

con M > 0. En nuestro caso, x(0) = e0M = 123, por tanto, M = 123 y lasolucion serıa x(t) = e2t123.

Ejemplo 6.12. Resuelve la ecuacion

y′ =−yx

, y(2) = −1.

Solucion. Si y′ = −yx , suponiendo que y = 0 se tiene que

y′

y=−1

x⇒

∫y′

ydx =

∫−1

xdx ⇒ L|y| = − L|x|+ a

con a ∈ R. Denotando b = ea, la solucion general se puede expresar ası:

|y| = 1

|x|b


siendo1 b > 0, es decir

y =c

xc ∈ R− 0.

De todas las soluciones recogidas en la solucion general nos interesa aquellaque verifica y(2) = −1, es decir:

y(2) =c

2= −1 ⇒ c = −2

Nuestra solucion particular serıa y = −2x .

Una solucion singular no recogida en la solucion general serıa la funcion

y = 0,

que habıamos descartado al razonar dividiendo por y en la ecuacion diferencialpero que tambien verifica la ecuacion diferencial.

6.3.1. Ecuaciones en variables separadas

Una EDO en variables separadas se puede expresar y′g(y) = f(x), que enforma normal quedarıa ası

y′ =f(x)

g(y)

y en forma diferencial como

F (x)G(y)dx + H(x)J(y)dy = 0

Se pueden resolver integrando por separado las expresiones que dependen dex y por otro las que dependen de y:∫

y′ g(y) dx =

∫f(x) dx

Ejemplo 6.13. Las ecuaciones x′(t) = 2x(t) y y′ = −yx , y(2) = −1 ya

resueltas son ejemplos de EDO en variables separadas. Otros ejemplos:

1. y′ = cos(3x) + 5 2. 2y dx + 3x dy = 0 3.(x2 + 4)y′ = xy

1Hay que tener en cuenta que b = ea > 0))

150 Capıtulo 6

6.3.2. Ecuaciones diferenciales homogeneas

Una EDO homogenea se puede expresar en forma normal ası

y′ = f(y/x),

es decir, la derivada y′ se puede expresar como una funcion que depende dey/x.

Ejemplo 6.14. Las ecuaciones diferenciales siguientes son homogeneas:

y′ =x2 + y2

xy(x2 − y2) + 3xyy′ = 0.

Las ecuaciones homogeneas se pueden resolver transformandolas en EDOen variables separadas mediante el cambio de variable y = u(x) · x, teniendoen cuenta que u(x) debe ser una funcion de x, por lo que

y′ = u′(x) · x + u

Ejemplo 6.15. Resuelve la ecuacion

(x2 − 3y2) dx + 2xy dy = 0

Solucion. La ecuacion (x2 − 3y2) dx + 2xy dy = 0 se puede expresar

y′ =−x2 + 3y2

2xy.

Dividimos por x2 numerador y denominador:

y′ =−1 + 3( yx)2

2 yx

,

aplicamos el cambio de variable y = u(x) · x, teniendo en cuenta que y′ =u′ · x + u:

u′ · x + u =−1 + 3u2

2u

y la ecuacion obtenida es de variables separadas

(u′ · x + u)2u = −1 + 3u2 ⇒ 2uu′x = −1 + u2 ⇒ 2uu′

−1 + u2=

1

x


que se resuelve integrando ambos terminos

L|u2 − 1| = L|x|+ c ⇒ |u2 − 1| = |x| · ec ⇒ u2 − 1 = ±xec

Si llamamos k = ±ec, la solucion serıa

u = ±√kx + 1, k = 0

Deshaciendo el cambio y = u(x) · x

y

x= ±√kx + 1 ⇒ y = ±x

√kx + 1

6.3.3. Ecuaciones reducibles a homogeneas

Una EDO que en forma normal se puede presentar ası

y′ = f

(ax + by + c

a′x + b′y + c′

)siendo a, b, c, a′, b′ y c′ numeros reales, se puede transformar en una ecuaciondiferencial homogenea mediante alguno de los siguientes cambios de variable:

1. Si las rectas y = ax + by + c, y = a′x + b′y + c′ se cortan en un punto(x0, y0), el cambio serıa

x = x′ + x0, y = y′ + y0.

2. Si las rectas y = ax + by + c, y = a′x + b′y + c′ son paralelas, el cambioserıa

u = ax + by.

Ejercicio 6.7. Resuelve las ecuaciones

1. y′ =6x− 7y − 4

3x + y − 22. y′ = e

4x+3y−25x−y+1 3. (x+y−2) dx+(x−y+4) dy = 0

152 Capıtulo 6

6.3.4. Ecuaciones diferenciales lineales

Una EDO lineal se puede presentar ası

y′ = p(x)y + r(x),

siendo p(x) y r(x) funciones de la variable independiente x.

Si r(x) = 0, la EDO lineal es de variables separadas

y′ = p(x)y

y tiene facil solucion

y′

y= p(x)⇔

∫y′

ydy =

∫p(x) dx⇔ L|y| =

∫p(x) dx + c.

Despejando y y llamando k = ec se obtiene

y = k e∫p(x) dx.

Para encontrar una solucion de la ecuacion lineal completa

y′ = p(x)y + r(x),

primero obtendremos la solucion de la ecuacion lineal correspondiente cuandocuando r(x) = 0, que sera

y = k e∫p(x) dx,

y despues buscaremos la solucion de la EDO lineal completa entre las funcionesde la forma

y = k(x) e∫p(x) dx,

exigiendo que k(x) sea una funcion de la variable x en lugar de constante.

Ejemplo 6.16. Resuelve la ecuacion:

y′ − y tg(x) = cos(x).

Solucion. Como es una ecuacion lineal, primero resolveremos la ecuacion

y′ − y tg(x) = 0


que es de variables separadas

y′

y= tg(x)∫

y′

ydx =

∫tg(x)dx

L|y| = −L| cos(x)|+ c

|y| = e−L| cos(x)|ec =1

eL| cos(x)|ec =

ec

| cos(x)|,

con lo que, llamando k = ±ec, la solucion es

y =k

cos(x)

Para obtener la solucion de la ecuacion lineal completa y′ − y tg(x) = cos(x),

imponemos que la solucion buscada sea de la forma y = k(x)cos(x) , siendo k(x)

una funcion de x. Como tiene que cumplir la ecuacion diferencial completa, setiene que

y′ − y tg(x) = cos(x)

(k(x)

cos(x))′ − k(x)

cos(x)tg(x) = cos(x)

k′(x) cos(x)− k(x)(−sen(x))

cos2(x)− k(x)

cos(x)

sen(x)

cos(x)= cos(x)

k′(x)

cos(x)= cos(x)

k′(x) = cos2(x)

k(x) =x

2+

sin(2x)

4+ C.

Por tanto la solucion general de la ecuacion es

y =k(x)

cos(x)=

x

2 cos(x)+

sin(2x)

4 cos(x)+

C

cos(x).

Ejercicio 6.8. Resuelve la ecuacion:

y′ + 2xy = 4x

154 Capıtulo 6

6.3.5. Ecuaciones de Bernoulli

Una EDO de tipo Bernoulli se puede expresar en forma normal ası:

y′ = p(x)y + q(x)yn

siendo n = 0, 1 porque en estos dos casos la EDO serıa lineal. En formadiferencial quedarıa ası:

dy = (p(x)y + q(x)yn) dx

Si dividimos toda la ecuacion por yn y aplicamos el cambio de variable

u =1

yn−1

la ecuacion de Bernoulli se transformara en lineal2:

dy

yn=

(p(x)

1

yn−1+ q(x)

)dx ⇔ du

1− n= (p(x)u + q(x)) dx

es deciru′

1− n= p(x)u + q(x)

que es la ecuacion lineal u′ = (1− n)p(x)u + (1− n)q(x).

Ejercicio 6.9. : Resuelve la ecuacion:

y′ + xy = x3y3.

6.4. Sistemas de primer orden y ecuaciones diferen-ciales de orden superior

6.4.1. Sistemas de ecuaciones de primer orden

La forma estandar de un sistema de ecuaciones diferenciales de primerorden es:

y′1 = f1(x, y1, . . . , yn)y′2 = f2(x, y1, . . . , yn)...y′n = fn(x, y1, . . . , yn).

2Aplicando el cambio anterior se obtiene que dyyn = du

1−n


En el sistema anterior hay n funciones a determinar y1, y2, . . . , yn, todas ellasdependientes de la variable independiente x. Denotamos y′i a la derivada de yirespecto de x.

Ejemplo 6.17. El siguiente es un sistema de ecuaciones diferenciales or-dinarias, siendo x e y las funciones a determinar y t la variable independiente.

x′ = x + 4y − et

y′ = x + y + 2et.

La solucion general exacta (no numerica) del sistema anterior esx = 2ae3t − 2be−t − 2et

y = ae3t + be−t + 14e

t.

Si al enunciado inicial le imponemos condiciones iniciales para las variablesdependientes, tendremos un problema de valor inicial.

Ejemplo 6.18. El siguiente es un sistema de ecuaciones diferenciales or-dinarias con valor inicial.

x′ = x + 4y − et

y′ = x + y + 2et

x(0) = 4 y(0) = 5/4.

La solucion del problema de valor inicial anterior esx = 4e3t + 2be−t − 2et

y = 2e3t − e−t + 14e

t.

SAGE resuelve sistemas de ecuaciones diferenciales de primer orden usandoel comando desolve_system([de1, de2], [x,y], ics=[0,1,2], ivar=t),siendo[ec1, ec2] las ecuaciones,[x,y] las variables dependientes,ics=[0,1,2]) (opcional) las condiciones iniciales de las variables indepen-diente y dependientes,[ivar=None] (opcional, por si no queda claro en la ecuacion) la variable in-dependiente.

Ejercicio 6.10. Obten la solucion del sistema de los Ejemplo 6.17 usandosoftware informatico.

156 Capıtulo 6

SAGE tambien resuelve sistemas de ecuaciones diferenciales de primer or-den con metodos numericos usando el comando desolve_odeint( ),por ejem-plo, el conjunto de comandosx,y=var(’x,y’)

f=[x*(1-y),-y*(1-x)]

sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y])

proporcional una lista de valores solucion numerica del sistema de ecuacionescon valor inicial siguiente:

x′ = x(1− y)y′ = −y(1− x).x(0) = 0,5 y(0) = 2

estando la variable independiente t valorada entre 0 y 10 con tamano de pasoh = 0,1.

Ejercicio 6.11. Obten una solucion numerica del sistema del Ejemplo6.18 usando software informatico.

6.4.2. Ecuaciones de orden superior

Toda ecuacion diferencial de orden superior del tipo

yn) = f(t, y, y, y′, y′′, . . . , yn−1)

donde la variable independiente es t y las derivadas de la variable dependientey se calculan derivando respecto de t, se puede convertir en un sistema deecuaciones de primer orden. Para ello introducimos unas nuevas variables xidefinidas como

x1 = y, x2 = y′, x3 = y′′, . . . , xn = yn−1)

de modo que las nuevas variables verifican el siguiente sistema de ecuacionesde orden 1

x′1 = x2x′2 = x3x′3 = x4...x′n = f(t, x2, x2, x3, x4, . . . , xn).

Este sistema se puede tratar como un sistema de ecuaciones de primerorden.


Ejercicio 6.12. Transforma el problema de valor inicial(sen(t)y′′′ + cos(ty) + sen(t2 + y′′) + (y′)3 = log ty(2) = 7 y′(2) = 3 y′′(2) = −4.

en un sistema de ecuaciones diferenciales de primer orden con valor inicial yresuelvelo con un metodo numerico.

Ejercicio 6.13. Transforma el problema de valor inicial(x′′)2 + tey + y′ = x′ − xy′y′′ − cos(xy) + sen(tx′y) = x

en un sistema de ecuaciones diferenciales de primer orden con valor inicial.Resuelvelo por metodos numericos con las condiciones iniciales

x(0) = 1, x′(0) = 0, y(0) = 0, y′(0) = 0,1



http://sage.unex.es

1. MUI_Ecuaciones_Diferenciales2. MUI_first order differential equations

3. MUI_Runge Kutta 4th

Problemas

Ecuaciones no lineales

1. Dada la ecuacion √xsen(x)− x3 + 2 = 0

encuentra una aproximacion de una raız

a) Con el metodo de la biseccion con un error menor que 1/30.

b) Con 4 iteraciones del metodo de Newton-Raphson.

2. Determina un intervalo y una funcion para aplicar el metodo del puntofijo a las siguientes ecuaciones:

a) x3 − x− 1 = 0

b) 4− x− tg(x) = 0

Realiza 4 iteraciones del metodo y determina en cada caso una cota delel error cometido.

3. Se considera la ecuacion x2 − 1− sen(x) = 0.

a) Prueba que dicha ecuacion tiene, al menos, una raız positiva.

b) Encuentra un intervalo en el cual la iteracion

xn =√

1 + sen(xn−1) x ∈ N

converja, para cualquier valor inicial x0 de dicho intervalo, a unaraız positiva de la ecuacion anterior. ¿Cuantos pasos deben darse,a partir de x0 = π/2, para obtener una aproximacion de la raız conun error inferior a una milesima?

159

160 Problemas

4. Demuestra que la ecuacion

exL(x) + x3 − 2 = 0

tiene una unica raız positiva. Aplica el metodo de Newton-Raphson paraencontrar una aproximacion a la raız. Realiza iteraciones hasta que serepitan las primeros cuatro decimales.

5. Determina un intervalo y una funcion para poder aplicar el metodo delpunto fijo a las siguientes ecuaciones:

a) x− L(1 + x)− 0,2 = 0.

b) x = −L(x).

Calcula una aproximacion a la solucion con un error inferior a 10−5.

6. Se considera la funcion

F (x) = e−x − sen(x)− 2.

a) Demuestra que F (x) < 0 si x ≥ 0.

b) Prueba que la ecuacion F (x) = 0 tiene una unica raız negativa.

c) Determina un intervalo donde se pueda aplicar el metodo de Newton-Raphson para aproximar dicha raız, ası como los 4 primeros termi-nos de la sucesion definida mediante dicho metodo. Calcula unacota del error cometido.

7. Dada la funcion f(x) = ex − 3x2,

a) Localiza un intervalo que contenga una solucion de f(x) = 0. Enun-cia el resultado teorico que utilices.

b) Encuentra, aplicando el metodo del punto fijo con cuatro iteracio-nes, una aproximacion de un valor donde se anule la derivada def(x) y una cota del error cometido.

8. Dada la ecuacion x = arctg(x+1), aplica el metodo de la biseccion paraencontrar una solucion con un error menor que 1/10.

9. Calcula tres iteraciones del metodo de la secante para encontrar una raızde la ecuacion

cos(x)− x = 0

en el intervalo [0,5, π/4]

Problemas 161

10. Dada la ecuacion 2cos(2x) + 4x− k = 0,

a) Determina el valor de k para que la ecuacion tenga una unica raıztriple en el intervalo [0, 1].

b) Para k = 3, prueba que posee una unica raız simple en el intervalo[0, 1] y calculala con el metodo de Newton con un error menor que10−3.

11. Dado el polinomio P (x) = x4 − 4x3 + 7x2 − 5x− 2

a) Utiliza la sucesion de Sturm para separar las raıces del polinomio.

b) Utiliza el metodo de Newton-Raphson para encontrar una aproxi-macion de la raız real de P (x) de menor valor absoluto.

12. Dado el polinomio P (x) = 8x3−20x2−2x+5, utiliza su sucesion de Sturmpara separar sus raıces en intervalos disjuntos. Encuentra aproximacionesde las raıces con un error menor que 10−3.

13. Dado la ecuacion x4 − 7x3 + 18x2 − 20x + 8 = 0:

a) Demuestra con el teorema de Bolzano que tiene una solucion en elintervalo [1,5, 3]. y utiliza el metodo de Newton para encontrar unaaproximacion con error menor que 10−6

b) Utiliza la sucesion de Sturm para separar las distintas raıces delpolinomio.

14. Cuestion de examen de febrero de 2004:Dada la ecuacion x2 − e2x = −1/2,

a) Demuestra que tiene una unica solucion en el intervalo [−2, 0].

b) Encuentra una aproximacion a la solucion mediante tres iteracionesdel metodo del punto fijo partiendo del punto x0 = −1,5

15. Cuestion de examen de junio de 2004:Dada la funcion f(x) = cos(x)− x,

a) Calcula el numero de pasos necesarios para encontrar una raız def(x) con un error menor que 10−2 con el metodo de la biseccionen el intervalo [0,5, π/4]. Aplica el metodo con el numero de pasoscalculado.

162 Problemas

b) Aplica tres iteraciones del metodo de la secante para encontrar unaraız de f(x) en el intervalo [0,5, π/4].

16. Cuestion de examen de septiembre de 2004:Dado el polinomio P (x) = −x3 + x2 − 2x + 5,

a) Utiliza el metodo de Sturm para separar las raıces en intervalosdisjuntos.

b) Empezando con x = 6 aplica un paso del metodo de Newton ha-ciendo uso del algoritmo de Horner.

Problemas 163

Sistemas de ecuaciones lineales

1. Para el conjunto de ecuaciones.

7x1 − 3x2 + 8x3 = −49x1 − 2x2 − 5x3 = 5

4x1 − 6x2 + 10x3 = −84

a) Calcula la inversa de la matriz de los coeficientes haciendo uso delas operaciones fundamentales. b) Utiliza la inversa para obtener lassoluciones del sistema. c) Calcula el numero de condicion de la matrizcon las distintas normas que conozcas.

2. Dado el sistema de ecuaciones

0,5x1 − x2 = −9,50,28x1 − 0,5x2 = −4,72

a) Resuelvelo por Gauss. b) Sustituye a11 por 0,55 y resuelvelo de nuevo.c) Interpreta los resultados anteriores en terminos del condicionamientodel sistema. d) Calcula el numero de condicion de la matriz del sistemacon la norma infinito.

3. Dado el sistema

−12x1 + x2 − 7x3 = −80x1 − 6x2 + 4x3 = 13−2x1 − x2 + 10x3 = 92

a) Resuelvelo por eliminacion gaussiana simple usando dıgitos con doscifras decimales. b) Encuentra las matrices triangular inferior L, trian-gular superior U y de permutaciones P tal que PA = LU , siendo A lamatriz de los coeficientes. c) Sustituye los resultados en las ecuaciones y,si es necesario, utiliza el metodo del refinamiento iterativo para mejorarla solucion.

4. a) Resuelve los siguientes sistemas de ecuaciones mediante eliminaciongaussiana con pivoteo parcial empleando dıgitos con dos decimales. b)Encuentra las matrices P , L y U tal que PA = LU como en el ejercicioanterior. c) Comprueba los resultados, y si no son exactos, emplea elrefinamiento iterativo para mejorar la solucion.

164 Problemas

1).4x1 + 5x2 − 6x3 = 28

2x1 − 7x3 = 29−5x1 − 8x2 = −64

2).3x2 − 13x3 = −50

2x1 − 6x2 + x3 = 444x1 + 8x3 = 4

5. Dado el sistema 3 2 100−1 3 1001 2 −1

x1x2x3

=

1051022

resuelvelo mediante descomposicion LU de Gauss: a) sin escalamiento yb) con escalamiento previo de filas.

6. Dados

A =

4 −3 02 2 36 1 −6

b =

−7−26

b′ =

149−8

resuelve Ax = b y Ax = b′ mediante una descomposicion LU de Gauss.Utiliza 3 decimales en los calculos.

7. Resuelve el sistema siguiente mediante descomposicion LU :

−4x1 + 2x2 = 0x1 − 4x2 = −4

x2 − 4x3 + x4 = −11x3 − 4x4 + x5 = 5

2x4 − 4x5 = 6

8. Resuelve el sistema4 −1 −1 0−1 4 0 −1−1 0 4 −10 −1 −1 4

x1x2x3x4

=

−404−4

mediante la descomposicion de Cholesky. ¿Es la matriz de los coeficientesdefinida positiva?

Problemas 165

9. Dado el sistema 2 0 −1−2 −10 0−1 −1 4

x1x2x3

=

1−12

2

comprueba si se cumple alguna de las condiciones suficientes para quefuncionen los metodos de Jacobi y Gauss-Seidel. Si se cumplen, realizados iteraciones para cada metodo tomando como valor inicial (0, 0, 0).

10. Repite el ejercicio anterior para el sistema

x1 − 2x2 = −13x1 + x2 = 4

Intenta resolver el sistema con los metodos de Jacobi y Gauss-Seidelhaciendo algun pequeno cambio en las ecuaciones.

11. Cuestion de examen de febrero de 2004:Aplica el metodo de Jacobi a la resolucion del sistema:

7x− y + 4z = 83x− 8y + 2z = −44x + y − 6z = 3

empleando 3 iteraciones y partiendo del punto (0, 0, 0). ¿Es convergenteel metodo?

12. Cuestion de examen de junio de 2004:Resuelve por el metodo de Cholesky el sistema de ecuaciones:

x1 + 2x2 + 3x3 = 72x1 + 5x2 + 4x3 = 9

3x1 + 4x2 + 14x3 = 33

13. Cuestion de examen de septiembre de 2004:Resuelve el siguiente sistema dos veces. Primero, utilizando la elimina-cion gaussiana simple y encontrando la factorizacion A = LU . Despues,utilizando la eliminacion gaussiana con pivoteo y encontrando la facto-rizacion PA = LU.

−x1 + x2 − 4x3 = 02x1 + 2x2 = 13x1 + 3x2 + 2x3 = 1

2

166 Problemas

Aproximacion de funciones. Interpolacion

1. De la funcion f(x, t) que representa una onda que se desplaza sobre elplano se conocen cinco posiciones en un determinado instante:

(−1, 0), (−0,5,−1), (0, 0), (0,5, 1), (1, 0).

a) Encuentra un polinomio que pase por esos puntos.

b) Encuentra una funcion formada por trozos de recta que pase poresos puntos.

c) Sabiendo que las derivadas de f(x) en el mismo instante en lospuntos−1,−0,5, 0, 0,5 y 1 son−π, 0, π, 0 y−π, calcula un polinomioque cumpla todas las condiciones.

d) Calcula el spline natural para ese conjunto de datos.

e) Encuentra un ajuste lineal y cuadratico por mınimos cuadrados.

2. ¿Cual es el polinomio interpolador de f(x) = x3 − 2x− 5 en los puntos−0, 1, 2, 4.? Encuentra el ajuste lineal, cuadratico, cubico y exponencialpor mınimos cuadrados de f(x) para los puntos anteriores.

3. Se ha estudiado el proceso de recepcion de informacion de un ordenadory se ha observado que en los instantes 0, 1, 2 medidos en segundos, lacantidad de informacion recibida era 0, 51, 104 megabites y se recibıaa una velocidad de 0, 52, 54 megabites por segundo, respectivamente.Encuentra una funcion que se aproxime a la que representa el proceso.

4. Dada la funcion f(x) = ex y los puntos 0, 1, 2.,

a) Encuentra su polinomio de interpolacion en los puntos y estima elerror cometido.

b) Repite el apartado anterior con los nodos de Chebyshev en el inter-valo [0, 2].

c) Encuentra una funcion lineal a trozos que se aproxime a f(x) enlos puntos dados.

d) Encuentra un polinomio que interpole a f(x) y a sus derivadas enlos puntos dados.

e) Calcula el spline natural en los puntos dados.

Problemas 167

5. La funcion f(x) = L(x) se ha evaluado en varios puntos, obteniendoselo siguientes datos:

x0 = 1 f(x0) = 0x1 = 4 f(x1) = 1,3863x2 = 6 f(x2) = 1,7918

a) Utilıcense los datos anteriores para calcular el polinomio de inter-polacion y el spline natural. Estımese el valor de L(2).

b) Agregando un cuarto punto x3 = 5, f(x3) = 1,6094, haganse denuevo los calculos y estimaciones.

6. Calculando los valores de la funcion f(x) =√x en 1,4 y 9,

a) Halla el polinomio interpolador de grado menor o igual que 2 porel metodo de Newton que pase por esos puntos.

b) Halla el polinomio interpolador en los nodos de Chebyshev en [1,9].

c) Aplica el metodo de Hermite en los puntos iniciales y en sus deri-vadas.

d) Calcula el spline natural.

7. Dada la funcion f(x) = ex, determina los cuatro nodos de Chebyshev en[0, 1], el polinomio interpolador en dichos nodos y un valor aproximadode e0,5 mediante dicho polinomio. Realiza los calculos con 4 decimales.

Polinomios de Chebysev: T0(x) = 1, T1(x) = x, Tn+1(x) = 2xTn(x) −Tn−1(x).

8. Determina si la funcion

f(x) =

1 + x− x3 si x ∈ [0, 1]

1− 2(x− 1)− 3(x− 1)2 + 4(x− 1)3 si x ∈ [1, 2]4(x− 2) + 9(x− 2)2 − 3(x− 2)3 si x ∈ [2, 3]

es el spline cubico natural que interpola a los puntos de la tabla:

x 0 1 2 3

y 1 1 0 10

168 Problemas

9. La solucion obtenida para una ecuacion diferencial y′ = f(x) medianteun metodo numerico es la siguiente tabla de valores:

x 0 0.1 0.2 0.3 0.4

y 1.2 1.5 1.9 2.1 2.6

Calcula el spline cubico natural y el ajuste exponencial mınimo cuadradopara aproximar la solucion.

Problemas 169

Diferenciacion e integracion numerica

1. De la funcion f(x) = ex − 2x2 + 3x− 1 se obtienen los siguientes datos:

x 0.0 0.2 0.4 0.6

f(x) 0.00000 0.74140 1.3718 1.9021

a) Utilizando la formula mas adecuada de entre las siguientes:

f ′(x) = f(x+h)−f(x)h f ′(x) = f(x+h)−f(x−h)

2h f ′′(x) = f(x+h)−2f(x)+f(x−h)h2 ,

calcula valores aproximados de la primera y segunda derivadas delos puntos.

b) Calcula los errores reales y las cotas de error por medio de lasformulas de error para las estimaciones anteriores.

c) Encuentra una aproximacion de f ′(0,2) haciendo uso de la interpo-lacion de la funcion en los cuatro puntos.

d) Sabiendo que f(−0,2) = −0,861267, aplica la interpolacion de Ri-chardson para calcular una aproximacion de f ′(0,2).

2. De una funcion f(x) se tienen los siguientes datos:

x 0.0 0.5 1 1.5 2 2.5 3

f(x) 2.25 3.75 5 6 5.5 6 7.25

Calcula un valor aproximado de∫ 30 f(x)dx utilizando los siguientes meto-

dos:

a) Trapecios.

b) Simpson 1/3.

c) Simpson 3/8.

d) Trapecios compuestos.

e) Simpson compuesto.

f ) Romberg.

3. Determina el numero de subintervalos que deben tomarse para aproximarel valor de la integral

∫ 10,5 cos(

√x)dx con un error menor que 1

210−2

a) con el metodo de los trapecios compuesto,

170 Problemas

b) con el metodo de Simpson compuesto.

4. Repite el ejercicio anterior para

(a)∫ 21 Ln(x)dx (b)

∫ 32

ex

x dx

5. Calcula∫ 10 sen(x2)dx:

a) con un error menor que 10−2 con el metodo de los trapecios com-puestos

b) con un error menor que 10−2 con el metodo de Simpson compuesto,

c) con el algoritmo de Romberg hasta R(4, 4),

d) con el metodo de cuadratura adaptativa con nivel de tolerancia10−3.

6. Repite el ejercicio anterior con∫ 21 e1/xdx.

7. a) Calcula la integral∫ 10 sen(x2)dx mediante cuadratura gaussiana de

dos terminos.

b) Plantea como resolver el problema de calcular la∫ 10 sen(x2)dx me-

diante cuadratura gaussiana de tres terminos. Calcula la integralsabiendo que los nodos (ti) en este caso son las raıces del polinomio5x3−3x

2 .

8. Aplica algun metodo de integracion numerica para resolver las siguientesintegrales:

(a)∫ 10

sen(x)x dx (b)

∫ 1−1

cos(x)−ex

sen(x) dx (c)∫∞1 (xex)−1dx.

9. Determina los valores de A, B y C que hagan que la formula∫ 2

0xf(x) ≈ Af(0) + Bf(1) + Cf(2)

sea exacta para todos los polinomios de grado tan alto como sea posible.¿Cual es el grado maximo?

10. Evalua la integral de x3sen(x) entre x = 0 y x = 2 usando un valor detolerancia inferior a 0,1 % de la respuesta verdadera (3,791197).

Problemas 171

11. Cuestion de examen de febrero de 2004:Dada la funcion f(x) = L(x) calcula

∫ 32 L(x)dx empleando el metodo de

cuadratura de Gauss con tres terminos y sabiendo que los polinomios deLegendre se calculan ası:

L0(x) = 1, L1(x) = x,

Ln+1(x) = (2n+1)xLn(x)−nLn−1(x)n+1 ∀ n ≥ 1.

12. Cuestion de examen de junio de 2004:Encuentra los coeficientes A0, A1, A2 y A3 para que la formula∫ 3

−1f(x)dx = A0f(−1) + A1f(0) + A2f(1) + A3f(2)

sea exacta para los polinomios de grado menor o igual que 3 y calculaun valor aproximado de

∫ 3−1(x

3 + x)dx haciendo uso de dicha formula.

13. Cuestion de examen de septiembre de 2004:Encuentra los coeficientes A1 y A2 para que la formula∫ 2π

0f(x)dx = A1f(0) + A2f(π)

sea exacta para cualquier funcion que tenga la forma a + b cos(x).

Calcula∫ 2π0 2 + 3 cos(x)dx haciendo uso de dicha formula.

172 Problemas

Resolucion numerica de ecuaciones diferencialesordinarias

1. Cuestion de examen de junio de 2004:Aplicando el metodo de Taylor de grado 2 con dos pasos, calcula valoresaproximados de y(0,1) y y(0,2), sabiendo que

y′ + 2y = x2 + 2x y(0) = 1

2. Cuestion de examen de febrero de 2005:

Resuelve la ecuacion y′ = y2 − t2 + 1, y(0) = 0,5, siendo 0 ≤ t ≤ 1,2,mediante el metodo de Taylor de grado 2 con 5 pasos.

3. a) Deduce la formula de Adams-Bashforth de orden 4:

yi+1 = yi +h

24

[55f(xi, yi)− 59f(xi−1, yi−1)+

+ 37f(xi−2, yi−2)− 9f(xi−3, yi−3)]

y usala para resolver la ecuacion

y′ =yx− x2

x2y(1) = 2

con h = 1/128

b) Repite el apartado anterior para la formula de Adams-Moulton deorden 4:

yi+1 = yi +h

24

[9f(xi+1, y(i + 1)) + 19f(xi, yi)−

− 5f(xi−1, yi−1) + f(xi−2, yi−2)]

4. Usa la identidad

y(xi+2)− y(xi) =

∫ xi+2

xi

y′(x)dx

y la formula de cuadratura de Simpson 1/3 para obtener la formula delsiguiente metodo multipaso

yi+2 = yi +h

3

[f(xi, yi) + 4f(xi+1, yi+1) + f(xi+2, yi+2)

].

Problemas 173

5. Aplica el metodo clasico de Runge-Kutta para estimar 10 valores apro-ximados de la solucion de la ecuacion siguiente en el intervalo [0, 1]

y′ + 2y = x2 + 2x y(0) = 1

6. A partir de la identidad

y(xi+1)− y(xi) =

∫ xi+1

xi

y′(x)dx

deduce el metodo de Runge-Kutta de orden 2

yi+1 = yi + h[12f(xi, yi) +

1


]

Iniciacion a la Investigaci on en Matematica Aplicada

Documents

Transcript of Iniciacion a la Investigaci on en Matematica Aplicada