Indice

1. Introducci
on
1.1. Error. . . . . . . . . . .
1.2. Dgitos significativos. . .
1.3. Polinomios de Taylor . .
1.3.1. F
ormula de Euler

.
.
.
.

4
4
6
8
8

.
.
.
.
.
.
.
.

10
10
11
14
15
19
21
21
22

3. Ceros de las funciones Polin

omicas.
3.1. Sistemas de Ecuaciones No lineales. . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2. Metodo de Newton- Raphson Generalizado. . . . . . . . . . . . . . . . . . . . . . .

25
27
28

4. Sistemas de Ecuaciones Lineales.

4.1. Metodos directos . . . . . . . . . . . . . . . . . .
4.1.1. Metodo de Gauss. . . . . . . . . . . . . .
4.1.2. Eliminaci
on gaussiana con y sin pivoteo.
4.2. Factorizaci
on LU . . . . . . . . . . . . . . . . . .
4.2.1. Matrices definidas positivas. . . . . . . . .
4.2.2. Matrices diagonal dominante. . . . . . . .

.
.
.
.
.
.

33
33
33
35
37
39
39

.
.
.
.
.
.
.

41
41
41
42
43
44
45
46

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

2. Ecuaciones No lineales.
2.1. Metodo de Bisecci
on . . . . . . . . . . . . . . . . . . . . . .
2.2. Metodo de la secante. . . . . . . . . . . . . . . . . . . . . .
2.3. Estudio de la convergencia del metodo de Newton-Raphson.
2.4. Velocidad de convergencia del metodo de Newton- Raphson
2.5. Metodos de punto fijo. . . . . . . . . . . . . . . . . . . . . .
2.6. Convergencia . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7. Velocidad de Convergencia . . . . . . . . . . . . . . . . . . .
2.8. Convergencia del Metodo de Newton-Raphson . . . . . . . .

5. M
etodos iterativos.
5.1. Normas vectoriales y matriciales. . . . . . .
5.1.1. Preliminares. . . . . . . . . . . . . .
5.1.2. Normas matriciales . . . . . . . . . .
5.1.3. Otros resultados importantes para el
5.2. Metodos iterativos. . . . . . . . . . . . . .
5.2.1. Metodo de Jacobi. . . . . . . . . . .
5.3. Metodo de Gauss-Seidel. . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

. . . . . .
. . . . . .
. . . . . .
estudio de
. . . . . .
. . . . . .
. . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

. .
. .
. .
los
. .
. .
. .

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

. . . . .
. . . . .
. . . . .
metodos
. . . . .
. . . . .
. . . . .

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

. . . . . .
. . . . . .
. . . . . .
iterativos.
. . . . . .
. . . . . .
. . . . . .

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.

5.4. Metodo de sobrerelajacion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.5. Resultados importantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6. An
alisis de error y N
umero de Condicion . . . . . . . . . . . . . . . . . . . . . . . .

47
48
49

6. Resoluci
on Num
erica del Problema de Valor de Contorno de orden 2.
6.1. Metodo de diferencias finitas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53
53

7. Sistemas de Ecuaciones No lineales.

7.1. Metodo de Newton- Raphson Generalizado. . . . . . . . . . . . . . . . . . . . . . .
7.2. Metodo de Broyden. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56
57
62

8. Interpolaci
on Polin
omica.
8.1. Metodo de los coeficientes indeterminados. . . . . . .
8.2. Metodo de Lagrange. . . . . . . . . . . . . . . . . . .
8.3. Forma de Newton para el Polinomio de Interpolacion
8.3.1. Diferencias Divididas . . . . . . . . . . . . . .
8.4. El problema de Hermite . . . . . . . . . . . . . . . .
8.5. C
omo obtener el polinomio de Hermite ? . . . . . .

.
.
.
.
.
.

66
66
67
69
70
74
74

.
.
.
.
.
.
.
.
.
.
.
.
.
.

78
78
78
80
82
87
88
89
90
93
95
97
99
104
106

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

9. Integraci
on Num
erica.
9.1. F
ormulas de Newton-Cotes. . . . . . . . . . . . . . . . .
9.1.1. Metodo del Trapecio. . . . . . . . . . . . . . . . .
9.1.2. Metodo de Simpson. . . . . . . . . . . . . . . . . .
9.2. Cuadratura Gaussiana. . . . . . . . . . . . . . . . . . . .
9.3. Integrales Indefinidas . . . . . . . . . . . . . . . . . . . . .
9.4. Integrales Impropias . . . . . . . . . . . . . . . . . . . . .
9.5. Integrales M
ultiples . . . . . . . . . . . . . . . . . . . . .
9.6. Diferenciaci
on Numerica . . . . . . . . . . . . . . . . . . .
9.7. Mnimos cuadrados Discretos. . . . . . . . . . . . . . . . .
9.8. Algo sobre Normas de funciones. . . . . . . . . . . . . . .
9.9. Aproximaci
on de Mnimos Cuadrados Continuos. . . . . .
9.9.1. Metodo Generalizado . . . . . . . . . . . . . . . .
9.9.2. Sistemas Ortonormales. . . . . . . . . . . . . . . .
9.9.3. Proceso de ortonormalizacion de Gramm-Schmidt.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

10.Aproximaci
on Minimax.
108
10.0.4. Caso Discreto. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
10.1. Algoritmo de Remez. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

10.1.1. Caso Continuo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

10.2. Propiedad mnimax de los polinomios de Tchebyshev . . . . . . . . . . . . . . . . . 113
10.3. Potencias de xn en terminos de {Tn (x)} . . . . . . . . . . . . . . . . . . . . . . . . 116
11.M
etodos de un paso.
117
11.0.1. Metodos de solucion por serie de Taylor. . . . . . . . . . . . . . . . . . . . . 117
11.0.2. Metodos por serie de Taylor de orden superior. . . . . . . . . . . . . . . . . 120
11.0.3. Metodos de Runge-Kutta.(R-K) . . . . . . . . . . . . . . . . . . . . . . . . 121
12.METODOS MULTIPASOS.
122
12.1. Deducci
on de las f
ormulas predictor corrector . . . . . . . . . . . . . . . . . . . . . 123
13.M
etodos predictor - corrector.
124
13.0.1. F
ormulas de Adams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
14.Ecuaciones elpticas.
14.1. La ecuaci
on de Laplace . . . . .
14.2. Soluci
on Numerica . . . . . . .
14.2.1. La Ecuaci
on de Laplace.
14.2.2. El metodo de Liebmann.
14.3. Condiciones en la frontera. . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

15.La ecuaci
on de conducci
on de calor.
15.0.1. Metodos Explcitos . . . . . . . . . .
15.0.2. Convergencia y estabilidad . . . . .
15.0.3. Aproximaciones temporales de orden
15.1. Un metodo implcito simple . . . . . . . . .
15.2. El metodo de Crank-Nicholson . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

. . . . .
. . . . .
superior
. . . . .
. . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

128
128
131
131
133
134

.
.
.
.
.

136
137
139
139
139
141

16.La ecuaci
on de ondas.
143
16.1. Construcci
on de la ecuacion en diferencias. . . . . . . . . . . . . . . . . . . . . . . 143
16.2. Condiciones Iniciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

UNIVERSIDAD DE SANTIAGO DE CHILE

Topicos I.
1.

Introducci
on

Profesor : Mara Angelica Vega U.

El prop
osito de estos m
odulos es entregar los conceptos te
oricos necesarios y las aplicaciones en forma dirigida, de modo que el estudiante pueda participar en forma activa,
en el proceso ense
nanza - aprendizaje, elaborando su propio material de estudio. El
logro deseable es que el estudiante aprenda a aprender.

1.1.

Error.

En el an
alisis de un algoritmo, es necesario considerar el error en la aproximacion obtenida a
la soluci
on del problema. Este error puede producirse por varios factores, en el ejemplo 1.1 por
redondeo de la aritmetica y termino del proceso infinito. Estos errores son parte inevitable del
proceso de soluci
on y no deben confundirse con los errores cometidos por el usuario, como falla en
el algoritmo o error en la programacion.
Para cada algoritmo podemos conocer una cota de error, esto es una cota sobre el error en la
soluci
on calculada.
Una cota que no depende de la solucion y puede ser calculada previamente, se llama cota de error
a priori. Cuando una cota no puede ser calculada hasta despues de la completacion de la parte
principal de los c
alculos es llamada una cota a posteriori. La cota de error es mas grande que el
error, por tal raz
on debemos desarrollar metodos alternativos de estimacion de error para algunos
problemas. Las cotas de error y estimativos son a menudo calculables y nos dan informacion tanto
del comportamiento asint
otico del error como del mejoramiento de las aproximaciones.
Existen dos formas de medir o cuantificar la magnitud de un error.
Definici
on 1.1 . Sea x una aproximaci
on del valor exacto x. Un indicador natural es la distancia
entre el valor exacto y el aproximado.
Definimos,
El error absoluto como la distancia,
EA (x ) = |x x |
Observamos que el error absoluto no permite comparar la precisi
on entre aproximaciones de
diferentes
ordenes de magnitudes.
4

El error relativo como,

ER (x ) =

EA (x )
|x x |
=
|x|
|x|

Este error es m
as significativo para comparar aproximaciones.
El error porcentual es el error relativo multiplicado por 100,
EP (x ) = 100 ER (x )
Estos conceptos de error tienen s
olo un valor formal,puesto que no se conoce el valor exacto x. En
la pr
actica, se determina una cota o estimativo del error.
Otro concepto importante es el de intervalo precisi
on, diremos que [a, b] es un intervalo precisi
on
de x, que denotaremos por I(x), si
P (x [a, b]) = 1
es decir, la probabilidad que x este en el intervalo de precision es la cierta.
Relaci
on entre error absoluto e intervalo de precisi
on
Dada una cota del error absoluto de x , es posible determinar el intervalo de precision de x.
Sea  una cota del error absoluto,
EA (x )  |x x |  x  x x + 
es decir, x [x , x + ] , luego I(x) = [x , x + ].
Inversamente, dado I(x) = [a, b] un intervalo de precision de x, es posible determinar una
cota del error absoluto de cualquier aproximacion x en [a, b].
Sea x [a, b] |x x | max (x a, b x ) = 
de donde,
EA (x ) .
Si x =

a+b
2

=

ba
2 .

Observamos que el error absoluto no es significativo, a menos que tengamos alg

un conocimiento
de la magnitud de x. A menudo se usa el error relativo en termino de porcentajes.
Otro punto importante de considerar al seleccionar un algoritmo es la eficiencia, que se mide de
acuerdo al n
umero de operaciones aritmeticas basicas requeridas para dar una exactitud en el
resultado.

1.2.

Dgitos significativos.

El uso de la aritmetica de las computadoras lleva al concepto de dgitos significativos. El dgito

principal o m
as significativo de un n
umero real x 6= 0, es el dgito no nulo mas a la izquierda
de su expresi
on decimal. Todos los dgitos, incluyendo los ceros a la derecha del dgito significativo
principal, son significativos y el u
ltimo de la serie es el menos significativo. Sin embargo, los ceros
a la izquierda del dgito significativo principal no son significativos. Por ejemplo, si x = 0.0670 el
dgito significativo principal es 6 y el dgito menos significativo es el cuarto de izquierda a derecha
cero.
Definici
on 1.2
Diremos que x aproxima a x con dgitos significativos (o cifras) ssi es el entero m
as
grande no negativo tal que:
|x x |
< 5 10
|x|

|x x | < 5 10 |x|.

(1)

ii)
i) Diremos que x aproxima a x con dgitos significativos (o cifras) ssi es el entero m
as
grande no negativo tal que:
|x x | < 0.5 10
donde es el exponente de la potencia de 10 que resulta al transformar x a la notaci
on
cientfica normalizada.
Ejemplo 1.1 Para que x aproxime a 1000 con cuatro cifras significativas, debe satisfacer :
|x 1000|
< 5 104 .
|1000|

(2)

Es decir, 999.5 < x < 1000.5

Ejemplo 1.2
Si x = 5.98 la notaci
on normalizada de x = 0.598 101 , es decir, = 1.
ii)
i) Si x = 0.5 y x = 4.9, determine el n
umero de dgitos significativos de la aproximaci
on.
Ejemplo 1.3 Aproxime con 4 dgitos significativos de precisi
on el cero de la funci
on f (x) =
etodo de Newton.Raphson.
x cos(x) que se encuentra en el intervalo [0, pi
2 ], usando m
sol. Graficamos en el intervalo [0, pi
2]

Figura 1: f (x) = x cos(x) x [0, 2 ]

Algoritmo de Newton-Raphson: Dada una aproximaci

on inicial x0 , generamos una sucesi
on
{xn )}
,
tal
que
x
se
determina
por:
n
n=1
xk+1 = xk

(f (xk )
.
f 0 (xk )

(3)

Observando el gr
afico podemos seleccionar x0 = 0.7854, presentamos los resultados de cada iteraci
on en la siguiente tabla
n

xn

|xn+1 xn |
< 5 104
|xn+1 |

|xn+1 xn | < 0.5 104

0
1
2
3

0.7854
0.7395361684
0.7390851781
0.7390851332

0.06201702305 > 0.0005

0.0006102007094 > 0.0005
0.0000006075078226 < 0.0005

0.0458638 > 0.00005

0.00045099033 > 0.00005
0.0000000449 < 0.00005

Luego la aproximaci
on requerida es x3 = 0.7390851
Actividades Libres 1.1
1.

Use el metodo de Newton para aproximar con precisi

on 104 las races de

x 0.8 0.2 sin (x) = 0 en [0, 2 ].

2.

Aproxime con precisi

on 105 las races de
x
x
e + 2 + 2cos(x) 6 = 0.
7

1.3.

Polinomios de Taylor

Teorema 1.1 Sean

i) f C n [a, b] y f (n+1) existe en [a, b).
ii) x0 [a, b]
entonces para toda x [a, b] existe (x) entre x0 y x tal que :
f (x) = Pn (x) + Rn (x)
donde:
Pn (x)

= f (x0 ) + f 0 (x0 )(x x0 ) + f 00 (x0 )

+f ( n)(x0 )

y
Rn (x)

(x x0 )2
(x x0 )3
+ f 000 (x0 )
+
2!
3!

(x x0 )n
n!

= f (n+1) ((x))

(x x0 )(n+1)
(n + 1)!

Pn (x) se llama el polinomio de Taylor de grado n para f alrededor de x0 y Rn (x) se llama

el residuo o error de truncamiento asociado con Pn (x). La serie infinita que se obtiene al
tomar el lmite de Pn (x), cuando n se denomina serie de Taylor para f alrededor de x0 .
El termino de error de truncamiento se refiere generalmente al error involucrado al usar sumas
finitas o truncadas para aproximar la suma de una serie infinita.
1.3.1.

F
ormula de Euler

Un metodo numerico para de resolver el Problema de valor inicial :

(
P.V.I.

y 0 = f (x, y)
y(x0 ) = y0 .

es el metodo de Euler. Si escribimos la serie de Taylor, con respecto a alg

un punto xi de la solucion
y(x), entonces tenemos:
y(x) = yi + yi0 (x xi ) + yi00

(x xi )2
(x xi )3
+ yi000
+
2!
3!

Entonces, podemos aproximar la solucion en xi+1 por:

yi+1 = yi + yi0 h + yi00
8

h2
h3
+ yi000
+
2!
3!

donde, h = xi+1 xi . Es claro que mientras mas terminos se consideren en el desarrollo de Taylor
mejor ser
a la aproximaci
on. Si consideramos dos terminos del desarrollo, obtenemos:
yi+1 = yi + yi0 h yi+1 = yi + f (xi , yi )h

i = 0, 1, 2,

(4)

Observese que (197) es el metodo de Euler, con un error de orden O(h2 ), esto significa que los
terminos no considerados contienen potencias de h con exponente mayor o igual que 2.
Error de local del m
etodo de Euler.
Consideramos:
yi+1 = yi + hf (xi , yi )
y suponemos y(x) en serie de Taylor con segunda derivada continua en un intervalo que contiene
a xi
h2
(5)
y(xi+1 ) = y(xi ) + y 0 (xi )h + y 00 ()
2!
donde (xi , xi+1 ). Restando ambas ecuaciones :
y(xi+1 ) yi+1 = y(xi ) yi + h[f (xi , y(xi )) f (xi , yi )] +

(h2 ) 00
y ()
2

suponiendo yi = y(xi ), obtenemos la formula para el error local:

ei+1 = y(xi+1 ) yi+1 =

(h2 ) 00
y ().
2

(6)

Como generalmente no se conoce el valor de , aunque teoricamente exista, no se puede calcular

este error , pero es posible determinar un estimativo, acotando esta expresion por el maximo que
alcanza la segunda derivada en el intervalo, es decir
|ei+1 | M

h2
2!

(7)

con
M=

|y 00 (x)|

max
x(xi ,xi+1 )

Observemos que el m
aximo M no depende de h, luego si se reduce h en un factor de 21 , la cota
1
1
1
derror se reduce en 4 y una reduccion de h en un factor de 10
, la cota derror se reduce en 100
.
Actividades 1.1 Para el problema de valor inicial siguiente
y 0 = 0.2xy

y(1) = 1

Use el metodo de Euler para obtener una aproximaci

on a y(1.5) usando h = 0.1 y h = 0.05.
Resuelva esta ecuaci
on analticamente usando alguno de los metodos ya estudiados y elabore una
tabla que incluya en los puntos seleccionados (Use h = 0.1 ) el valor exacto, el valor aproximado
y el error relativo para cada aproximaci
on .

2.

Ecuaciones No lineales.
Modulo 2
Prof. M.Angelica Vega

2.1.

M
etodo de Bisecci
on

Problema 1. Dada una funcion f (x) no lineal y continua en el intervalo [a, b]. Determinar una
raz [a, b] de la ecuaci
on f (x) = 0.
Ante el problema planteado surgen las siguientes interrogantes :
1. Existe soluci
on?
2. Si existe soluci
on . Es u
nica ?
3. C
omo determinarla ?
1.- Existencia : Teorema de Bolzano.
La existencia de soluci
on est
a garantizada por el teorema de Bolzano, cuyo enunciado afirma que
bajo las hip
otesis :
i) f continua en [a, b].
ii) f (a) y f (b) de distinto signo. (es decir, f (a) f (b) < 0 )
entonces,
existe (a, b) tal que f () = 0.
2.- Unicidad :
La existencia de una u
nica solucion en el intervalo de precicion queda garantizada por el siguiente
teorema visto en C
alculo Diferencial .
Teorema 2.1 Si f es una funci
on que satisface las hip
otesis :
i) f continua en [a, b].
ii) f (a) f (b) < 0.
iii) f es diferenciable en (a, b)
iv) Signo de f 0 (x) es constante en (a, b) entonces,
existe un u
nico (a, b) tal que, f () = 0.
M
etodo de Bisecci
on. Se basa en el teorema de Bolzano, por lo tanto se supone :
i) f continua en I() = [a, b].
ii) f (a) f (b) < 0
Algoritmo de la Bisecci
on .
1.- Dividir el intervalo de precision I() por la mitad, obteniendose:
x1 =

(a + b)
2
8

2.- Analizar f (x1 )

Si f (x1 ) = 0 , entonces = x1 y detener el proceso.
Si f (x1 ) 6= 0 , entonces
Si f (x1 ) f (a) < 0 ,
Si

f (x)

tiene un cero en el intervalo

f (x1 ) f (a) > 0 ,

f (x)

[a, x1 ].

tiene un cero en[x1 , b].

3.- Repitiendo los pasos 1 y 2 se genera una sucesion de aproximaciones, hasta que se satisfaga una
condici
on de termino.
Ventajas y desventajas.
Es posible demostrar que si f (x) es continua en [a, b] y f (a)f (b) < 0, el procedimiento de biseccion
genera una sucesi
on que converge a la raz y
|xn |

(b a)
.
2n

(8)

Sin embargo este metodo tiene la desventaja que la convergencia es lenta .

Del resultado anterior, se puede deducir una formula para determinar el n
umero de iteraciones,
dada una cierta tolerancia .
ln ((ba)

n
.
(9)
ln(2)
Actividades 2.1
Demuestre que de (8) puede obtener (9).
Sea f (x) = x3 + 4x2 10, encuentre una raz aproximada.
a) C
omo podra verificar si es la u
nica?
b) Determine I(x).
c) C
uantas iteraciones se deben realizar para que la aproximaci
on sea correcta a lo menos
en 4 cifras significativas?

2.2.

M
etodo de la secante.

Se basa en encontrar una aproximacion a la raz mediante rectas secantes a la grafica de f (x), de
modo que, si (xk1 , f (xk1 )) y (xk , f (xk )) son dos puntos de la grafica de f (x) una aproximacion
a la raz es, xk+1 intersecci
on de la recta secante que une los puntos mencionados con el eje X.
De donde se deduce :
(f (xk ) (xk xk1 ))
xk+1 = xk
.
(10)
(f (xk ) f (xk1 ))

Figura 2:

Algoritmo de la Secante .
1.- Dadas dos aproximaciones x0 y x1 .
2.- Generamos una sucesi
on de aproximaciones {xk+1 } definida por:
xk+1 = xk

(f (xk ) (xk xk1 ))

.
(f (xk ) f (xk1 ))

M
etodo de la Regula Falsi.
Este algoritmo es una mezcla de los metodos de biseccion y secante.
1.- Si [a, b] es el intervalo de precision, unimos (a, f (a)) y (b, f (b)) para obtener x1 como interseccion
de la recta secante que pasa por los puntos precedentes y el eje X.
2.- Trazamos la recta secante por (x1 , f (x1 )) y (b, f (b)) cuya interseccion sobre el eje X determina
dos subintervalos, escogemos aquel que contiene a la raz, verificando :
Si f (x1 ) f (b) < 0 , la raz est
a en [x1 , b] , de lo contrario la raz pertenece al intervalo [a, x1 ] .
3.- Supongamos que f (x1 ) f (b) < 0 , repetimos el paso 1 considerando el intervalo [x1 , b] , de
modo que se obtiene la aproximacion x2 , como interseccion de la secante y el eje X.
4.- Procediendo as sucesivamente obtenemos la sucesion {xk+1 } definida por :
xk+1 = xk

(f (x[k]) (b xk ))
(f (b) f (xk ))

Gr
aficamente,
Actividades 2.2
10

(11)

Figura 3:

Usando el metodo de la secante, obtener el cero de

f = x 2 sen x
que est
a en el intervalo [1,5, 2]
Dada f = x3 + x2 3x 3 , determine en forma aproximada la raz que se encuentra en el
intervalo [0,3] .
M
etodo de Newton-Raphson.
El metodo de Newton-Raphson es uno de los metodos numericos mas conocidos en la resolucion
del problema de b
usqueda de races de la ecuacion :
f (x) = 0
Existen tres formas de derivar este metodo, un enfoque basado en el polinomio de Taylor, como
un procedimiento para obtener una convergencia mas rapida que la que se obtiene con los otros
metodos y el metodo gr
afico.
Discutiremos este u
ltimo. Geometricamente consiste en dada una funcion dos veces continuamente
diferenciable, nos aproximamos a la raz mediante rectas tangentes, de modo que si ((xk , f (xk )) es
un punto de la gr
afica de f (x), trazamos la recta tangente en este punto y donde la tangente corta
al eje X obtenemos la aproximacion xk+1 , es decir, geometricamente
de donde deducimos
(f (xk ))
xk+1 = xk 0
.
(12)
(f (xk ))

11

Figura 4:

Actividades 2.3 Use el metodo de Newton-Raphson para determinar una aproximaci

on de la
mayor raz negativa de la ecuaci
on f (x) = 0 con
f (x) = x5 + 0,85x4 + 0,70x3 3,45x2 1,10x + 1,265
use como aproximaci
on inicial x0 = 0,1.

2.3.

Estudio de la convergencia del m

etodo de Newton-Raphson.

Definici
on 2.1 Si {xn } es una sucesi
on de n
umeros reales que converge a la raz r, diremos que
la raz
on de convergencia es de orden si existen dos constantes C y y un entero N tal que,
|xn+1 r| C|xn r|

, n N

lm

|xn+1 r|
= C.
|xn r|

(13)

Observaci
on 2.1
i) Si usamos la notaci
on |en+1 | = |xn+1 r| (error absoluto cometido en la
(n+1)- esima aproximaci
on), el lmite puede abreviarse por
lm

|en+1 |
= C.
|en |

ii) Si C < 1 y N es un entero tal que

|xn+1 r| C|xn r|
diremos que la convergencia es a lo menos lineal .
12

n N

iii) Si existen C no necesariamente menor que 1 y un entero N tal que

|xn+1 r| C|xn r|2

n N

diremos que la convergencia es a lo menos cuadr

atica.

2.4.

Velocidad de convergencia del m

etodo de Newton- Raphson

Comentamos antes, que otra forma de obtener Newton-Raphson consiste en derivarlo a partir del
desarrollo en serie de Taylor de la funcion f (x).
En efecto, supongamos que desarrollamos en serie la funcion no lineal f (x) en torno del punto xk
= f (xk ) + f 0 (xk )(x xk ) +

f (x)

f 00 ()
2! (x

xk )2

evaluamos en xk+1
(14)

f (xk+1 )

= f (xk ) + f (xk )(xk+1 xk ) +

f 00 ()
2! (xk+1

xk ) .

donde, (xk , xk+1 ). Truncando el tercer termino, queda

f (xk+1 ) = f (xk ) + f 0 (xk )(xk+1 xk ).

(15)

En la intersecci
on con el eje X, f (xk+1 ) = 0, luego
0 = f (xk ) + f 0 (xk )(xk+1 xk ).

(16)

Ordenando, se obtiene
xk+1 = xk

f (xk )
.
f 0 (xk )

(17)

La f
ormula del metodo de Newton.
Adem
as este desarrollo permite estimar el error de la formula y la velocidad de convergencia del
metodo.
En efecto, supongamos que en el desarrollo anterior xk+1 es xr el valor exacto, entonces f (xr ) = 0.
y
f 00 ()
0 = f (xk ) + f 0 (xk )(xr xk ) +
(xr xk )2 .
(18)
2!
Restando (16) y (18) se obtiene :

ek+1 =

= f 0 (xk )(xr xk+1 ) +

= ek+1 f 0 (xk ) +

f 00 ()
2! (xr

xk )2

f 00 ()
2
2! (ek ) .

f 00 (xk ) 2
e
2f 0 (xk ) k

De esta u
ltima igualdad se deduce que el error es casi proporcional al cuadrado del error anterior
lo que significa que el n
umero de cifras decimales correctas se duplica en cada iteracion.
13

Tomando valor absoluto y lmite cunado k , se tiene que

|f 00 ()|
|ek+1 |
=
.
2
k |ek |
2|f 0 (xk )|
lm

(19)

De aqu deducimos que el metodo converge cuadraticamente, si f 0 (xk ) 6= 0, es decir si la raz es

simple.
Teorema 2.2 Si f C 2 [a, b] y r es una raz es simple en [a, b] (es decir, f (r) = 0 y f 0 (r) 6= 0 )
entonces, existe un > 0 tal que el metodo de Newton-Raphson generar
a una sucesi
on que converge
cuadr
aticamente a r para cualquier aproximaci
on inicial x0 [r , r + ] .
Deducimos que no ocurre necesariamente convergencia cuadr
atica si la
Observaci
on 2.2
raz no es simple.
El metodo de Newton-Raphson en general es muy eficiente, sin embargo hay situaciones en
que no es aplicable o se comporta en forma deficiente. Uno de estos casos es cuando la funci
on
posee races m
ultiples.
La siguiente figura muestra gr
aficamente, un ejemplo donde no es posible aplicar NewtonRaphson. Podra explicar porque ?

Figura 5: f (x) = x cos(x) x [0, 2 ]

Actividades 2.4 Los problemas que se presentan a continuaci

on son 2 casos especiales, donde no
es posible aplicar Newton-Raphson

14

Dada la funci
on
1
5
f (x) = x3 + x
2
2
aproxime la raz positiva, use como aproximaci
on inicial x0 = 1. y realice 6 iteraciones.
Que sucede ? Porque ocurre esto ?
Dada la funci
on
f (x) = (x 1)3 (x 2)
aproxime la menor raz, use como aproximaci
on inicial x0 = 1,75
Que sucede ? A que atribuye este hecho ?
Races m
ultiples.
Las races m
ultiples presentan algunas dificultades para los metodos que hemos analizado.
i) Si una funci
on no cambia de signo en races m
ultiples pares impide el uso de los metodos
cerrados que usan intervalos y los metodos abiertos tienen la limitacion que pueden ser
divergentes.
ii) Otra dificultad que puede presentarse, es el hecho que f (x) y f 0 (x) tiendan a cero. En estos
casos, hay dificultades con el metodo de Newton-Raphson y el de la secante (aproximacion
de f 0 (x)) que provocan una division por cero.
iii) Ralston y Rabinowitz demostraron que los metodos de Newton-Raphson y secante convergen linealmente cuando hay races m
ultiples y han propuesto algunas modificaciones que se
obtienen a partir del metodo de newton-Raphson para problemas con funciones que tienen
races m
ultiples.
iv) El siguiente metodo es una modificacion de Newton-Raphson, tiene convergencia cuadratica:
xk+1 = xk m

f (xk )
.
f 0 (xk )

(20)

donde, m es la multiplicidad de la raz, es decir m = 2 si la raz tiene multiplicidad 2.

v) Otra modificaci
on de Newton-Raphson, se obtiene definiendo una nueva funcion
u(x) =

f (x)
,
f 0 (x)

que tiene races en las mismas proporciones que la funcion original. Se sustituye u(x) en (17)
y se obtiene:
f (xk ) f 0 (xk )
xk+1 = xk 0
.
(21)
[f (xk )]2 f (xk ) f 00 (xk )

15

Actividades 2.5 Usar metodo de Newton-Raphson y el modificado (21)para aproximar la raz

m
ultiple de :
f (x) = x3 5x2 + 7x 3
use como aproximaci
on inicial x0 = 0. Realice a lo menos 5 iteraciones con el metodo de NewtonRaphson.
Nota 2.1 El material bibliogr
afico usado para la elaboraci
on de este taller :
An
alisis Numerico. R.L. Burden-J.D. Faires
Metodos numericos. S.C. Chapra- R.P. Canale

16

2.5.

M
etodos de punto fijo.

Modulo 3
Prof: Mara Ang
elica Vega U.
Sea f (x) = 0 una ecuaci
on no lineal. Los metodos para determinar una raz de esta ecuacion, que
se expresan en la forma
x = g(x)
se denominan M
etodos de Iteraci
on Funcional o M
etodos de punto fijo. Una solucion de
esta ecuaci
on se llama punto fijo de g y la funcion g(x) se llama funci
on de iteraci
on . Frente
al problema de buscar los puntos fijos de la ecuacion x = g(x), surgen las siguientes interrogantes :
1.- C
uando una funci
on g(x) tiene un punto fijo ?
2.- Si existe, es u
nico ?
3.- C
omo determinarlo ?
Las respuestas a la primera y segunda interrogante, es el resultado que veremos a continuacion
formalizado como Teorema de existencia y unicidad.
Teorema 2.3 Teorema de Existencia y unicidad. Si
i) g(x) es continua en [a, b].
ii) g(x) pertenece a [a, b], para todo x que pertenece al intervalo [a, b]. entonces, existe un punto
fijo para g(x) en (a, b).
Si adem
as,
iii) g(x) es diferenciable en (a, b)
iv) |g(x)| k < 1 , para todo x perteneciente al intervalo (a, b), entonces , g(x) tiene un u
nico
punto fijo en [a, b].
Demostraci
on.
1) Existencia. Si g(a) = a y g(b) = b la existencia de un punto fijo es trivial.
Suponemos entonces,
a < g(a) y b > g(b) (Por hip
otesis g(x) [a, b]).
Definimos:
h(x) = g(x) x
observamos que,
i) h(x) es continua en [a, b]
ii)
h(a) = g(a) a > 0
h(b) = g(b) b < 0

17

)
h(a) h(b) < 0.

Por teorema de Bolzano existe p (a, b) tal que

h(p) = 0 g(p) p = 0 p = g(p)
, luego p es punto fijo.
2) Unicidad. Supongamos que p y q son dos puntos fijos distintos. Por Teorema del Valor Medio,
( f continua y diferenciable, existe (p, q) tal que f 0 () =

f (b)f (a)
.)
ba

Tenemos,

|p q| = |g(p) g(q)| = |g 0 ()||p q| < k|p q| < |p q|.

Contradicci
on de suponer p 6= q, luego p = q en [a, b].
Algoritmo de punto fijo.
1.- Dada una aproximaci
on inicial x0 ,
2.- Generamos una sucesi
on de aproximaciones {xn } definida por :
xn = g(xn1 )
. Situaci
on gr
afica de algunas funciones.

Figura 7: x = x( 1/3)

Figura 6: x = (1/3)x + 0,2

Actividades 2.6 La ecuaci

on x3 + 4x2 10 = 0 tiene una raz en [1, 2], obtenga al menos 5
funciones de iteraci
on.
Convergencia
Teorema 2.4 Bajo las hip
otesis del teorema de existencia y unicidad, si p0 [a, b] la sucesi
on
definida por
pn = g(pn1 )
18

Figura 8: x = (3/2)x

converge al u
nico punto fijo. Adem
as se cumple
i)|pn p|

kn
|p0 p1 |
1k

n 1.

ii)|pn p| k n |p0 p| k n max {p0 a, b p0 },

puesto que, p [a, b].
Demostraci
on. Hacerla como ejercicio. Se utilizan las hip
otesis del Teorema de existencia y
unicidad y el Teorema del Valor Medio.
Actividades 2.7 Sea g(x) =
(1, 1).

(x2 1)
3

en [1, 1] . Demostrar que g(x) tiene un u

nico punto fijo en

Velocidad de Convergencia.
Teorema 2.5 Sea p un punto fijo de g(x), tal que p I y sea
xk = g(xk1 )

para

k = 1, 2, ...

i) Si g C 1 (I(p)) y g 0 (p) 6= 0, entonces

|xk+1 xk | g 0 (p)|xk xk1 |
de modo que, si x0 I(p) entonces,
{xk } converge linealmente a p con C = g 0 (p) si 0 < |g 0 (p)| < 1.
19

{xk } diverge linealmente a p con C = g 0 (p) si |g 0 (p)| > 1.

ii) Si g 00 (x) es continua y g 0 (p) = 0 entonces,
1
|xk+1 xk | g 00 (p)|xk xk1 |
2
si x0 I(p) entonces,
{xk } converge cuadr
aticamente a p con C = 12 g 00 (p).
Si g 0 (p) = 1 , {xk } puede converger o diverger muy lentamente.
Demostraci
on. Consultar An
alisis Numerico de M.J.Maron-R.J.L
opez.

2.6.

Convergencia del M
etodo de Newton-Raphson

Teorema 2.6 Si es una raz simple de f (x) = 0, entonces el algoritmo de Newton Raphson
tiene convergencia Cuadr
atica.
Sol. raz simple de f (x) = 0 implica que:
f (x) = (x ) h(x) con h() 6= 0
f 0 (x) = h(x) + (x ) h0 (x)
f 0 () = h() 6= 0
Desarrollando f (x) en serie de Taylor alrededor de x = se tiene :
f (x) = f () + f 0 ()(x ) +

f 000 ()(x )2
f 00 ()(x )2
+
2!
3!

Para alg
un en un entorno de . Consideremos el algoritmo de Newton-Raphson:
f (xn )
f 0 (xn )

xn+1

= xn

xn+1

= xn

xn+1

lmn

|xn+1 |
|xn |2

1
f 0 (xn )

f 0 ()(xn ) +

f 00 ()(xn )2
f 000 ()(xn )3
+
2!
3!



(xn )2 f 0 (xn ) f 0 () f 00 () f 000 ()(xn )

f 0 (xn )
(xn )
2!
3!

 

 1
f 00 ()   1 
=  0
(f 00 ()
) =  0
= cte
f ()
2
2f () 
(22)

Luego el metodo de Newton tiene orden de convergencia 2 para raz simple.

Teorema 2.7 Si es una raz multiplicidad p de f (x) = 0, entonces el algoritmo de Newton
Raphson tiene convergencia lineal.
Dem. raz simple de f (x) = 0 implica que:
20

f (x) = (x )p g(x) con g() 6= 0

f 0 (x) = p(x )p1 g(x) + (x )p g 0 ()
En efecto,
f (xn )
f 0 (xn )

xn+1

= xn

xn+1

= xn

xn+1

xn+1

lmn

|xn+1 |
|xn |2

(xn )p g(xn )
p(xn )p1 g(xn ) + (xn )p g 0 (xn )

(xn )p g(xn )
(xn )p1 (pg(xn ) + (xn ) g 0 (xn ))


g(xn )
= (xn ) 1
pg(xn ) + (xn )g 0 (xn )



g 0 () 

= lmn 1
g() 6= 0
pg() 
= (xn )

=1

(23)

1
p1
=
, con p > 1
p
p

Luego el metodo de Newton-Raphson tiene orden convergencia 1 o lineal para raz m

ultiple.
Actividades 2.8
1.

Considere la ecuaci
on:

e 2 2x = 0

(24)

i) Verifque la ecuaci
on (24) tiene una raz en el intervalo [4, 5].
ii) Demuestre que en el intervalo [4, 5] est
a la mayor raz positiva de la ecuaci
on (24).
iii) Para obtener la raz que est
a en [4, 5], se propone el siguiente algoritmo de punto fijo:
xk+1 =

xk
2

(25)
2
Determine justificadamente , si es posible garantizar la convergencia del algoritmo para
hallar la raz . Usando (25) y una aproximaci
on x0 , adecuada, trate de obtener una
aproximaci
on de la citada raz.
iv) Proponga un algoritmo de punto fijo convergente, y u
selo con x0 conveniente, para
determinar una aproximaci
on de la raz de (24) en [4, 5].
2.

Puede demostrarse que la forma que adopta un cable suspendido de sus extremos es catenaria.
Una ecuaci
on de una catenaria est
a dada por:
y = cosh
21

x
a

donde a es una constante por determinar. El eje Y equidista de sus puntos extremos. Un
cable telef
onico est
a suspendido, a la misma altura, de dos postes que est
an a 200 pies uno
del otro. El cable tiene una cada m
axima de 10 pies en el punto medio del cable. Determine
el valor de a. (Recuerde que cosh0 = 1.)
Nota 2.2 El material bibliogr
afico usado para la elaboraci
on de este taller :
An
alisis Numerico. R.L. Burden-J.D. Faires
An
alisis Numerico. M.J. Maron - R.J. L
opez.

22

3.

Ceros de las funciones Polin

omicas.

Una expresi
on de la forma:
Pn (x) = an xn + an1 xn1 + a1 x + a0
donde, los ai son constantes, an 6= 0, se llama polinomio de grado n.
Recordamos el metodo de Horner:
Teorema 3.1 Metodo de Horner.
Sea
Pn (x) =

n
X

ak xk

an neq 0

an = bn

k=0

i) Si bk = ak + bk+1 x0

k = n 1, 1, 0 entonces
P (x0 ) = b0

ii) Adem
as si,
Qn1 (x) = bn xn1 + bn1 xn2 + b2 x + b1
entonces
Pn (x) = (x x0 )Qn1 (x) + b0
Dem.
ii) Desarrollamos Q(x) en la expresi
on :
(x x0 )Qn1 (x) + b0

= (x x0 )(bn xn1 + bn1 xn2 + b2 x + b1 ) + b0

= bn xn + bn1 xn1 + b2 x2 + b1 x) (bn x0 xn1 + bn1 x0 xn2 + b2 x0 x + b1 x0 ) + b
= bn xn + (bn1 bn x0 )xn1 + + (b2 b3 x0 )x2 + (b1 b2 x0 )x + (b0 b1 x0 )
= an xn + an1 xn1 + a1 x + a0 = Pn (x).

i) Usando ii)
Pn (x) = (x x0 )Qn1 (x) + b0

P (x0 ) = b0

Corolario. P 0 (x0 ) = Q(x0 )

Dem.
= (x x0 )Q(x) + b0
= Q(x) + (x x0 )Q0 (x)
P 0 (x0 ) = Q(x0 )

P (x)

Observaci
on 3.1 Si se usa el metodo de Newton - Raphson para determinar las races de un
polinomio P (x), para evaluar P (xk ) y P 0 (xk ) puede usarse el metodo de Horner.

23

Ejemplo 3.1 Usando Newton - Raphson y Metodo de Horner, aproximar la mayor raz negativa
de
2x4 3x2 3x 4
Usando x0 = 2 como aproximaci
on inicial, relizar 3 iteraciones.
Soluci
on .
Primera iteraci
on.
i) C
alculo de P (x0 )
ai

a4
2

bi

a3
0
4
4

a2
3
8
5

a1
3
10
7

a0
4
14
10

Por teorema de Horner:

P (2) = 10

Q(x) = 2x3 4x2 + 5x 7

ii) C
alculo de P 0 (x0 )
Para determinar P 0 (2) = Q(2), usamos nuevamente metodo de Horner.
2
2

4
4
8

5
16
21

7
42
49

Luego P 0 (2) = 49 Usamos Metodo de Newton-Raphson

x1 = 2

10
1,796
49

Segunda y Tercera iteraci

on.
Actividades 3.1 Separar las races de los siguientes polinomios:
a) P3 (x) = 3x3 + 2x 1
c) P4 (x) = x4 2x2 + 1

24

UNIVERSIDAD DE SANTIAGO DE CHILE

CALCULO NUMERICO

Ecuaciones No lineales.
M
odulo 5
Prof: Mara Ang
elica Vega U.

3.1.

Sistemas de Ecuaciones No lineales.

Un sistema de n ecuaciones con n incognitas x1 , x2 , ..., xn se dice no lineal si una o mas ecuaciones
es no lineal.
Notaci
on.

f1 (x1 , ..., xn ) = 0

f2 (x1 , ..., xn ) = 0
..
.. o
. =.

fn (x1 , ..., xn ) = 0

f1 (x) = 0
x1

x2

f2 (x) = 0

.
..
.. donde x =
..
. =.

xn .
f (x) = 0
n

~ = ~0 (o F(x) = 0 .) Notese que la

Este sistema puede escribirse abreviadamente como, F~ (X)
soluci
on es un vector que debe satisfacer simultaneamente todas las ecuaciones del sistema.
Ejemplo 3.2 Analizar, si es posible determinar la soluci
on de :
i)
yex 2 = 0
2
x +y4 =0
ii)
xey x5 + y + 3 = 0
x + y + tan x sen y = 0
Soluci
on.
i) Seg
un la notaci
on anterior:
f1 (x, y) = yex 2 = 0
f2 (x, y) = x2 + y 4 = 0
25

Figura 9: x = (1/3)x + 0,2

Gr
aficamente, Una forma simple de resolver este sistema es, usar los metodos empleados en la
resoluci
on de sistemas lineales, es decir, despejar y en ambas ecuaciones, obteniendose:
x2 + 2ex 4 = 0
Usando cualquier metodo para aproximar races de ecuaciones no lineales se obtiene:
x1 (0,6, 3,7)T

x2 (1,9, 0,4)T

ii) Es posible usar el metodo anterior para determinar una soluci

on de:
f1 (x, y) = xey x5 + y + 3 = 0
f2 (x, y) = x + y + tan x sen y = 0
Justificar la respuesta.

3.2.

M
etodo de Newton- Raphson Generalizado.

La idea es ahora extender el metodo de Newton- Raphson, visto para funciones de una variable
real a funciones de dos o m
as variables, para resolver el sistema:
(
f (x, y) = 0
F(x) = 0
(26)
g(x, y) = 0
Analizaremos el caso m
as simple, es decir, funciones de dos variables. Para ello recordamos el
teorema de Taylor, para funciones de 2 variables .
Sea u = f (x, y) y todas sus derivadas parciales de orden menor o igual a (n + 1) continuas en
una regi
on D = {(x, y)/a x b c y d}. Sea (x0 , y0 ) D, para todo (x, y) D existen
26

 (x, x0 ) y (y, y0 ) tal que:

Pn (x, y)

= f (x0 , y0 ) + (x x0 )

f
f
(x0 , y0 ) + (y y0 ) (x0 , y0 )+
x
y

(x x0 )(y y0 ) 2 f
(x x0 )2 2 f
(x0 , y0 ) +
(x0 , y0 )+
2
2!
x
2!
xy

(y y0 )2 2 f
(x0 , y0 ) + +
2!
y 2

(27)

n  
1 X n
nf
=
(x x0 )nj (y y0 )j nj j (x0 , y0 )+
n! j=0 j
x
y


n+1
X n+1
1
n+1 f
(x x0 )n+1j (y y0 )j n+1j j (, ).
=
j
(n + 1)!
x
y
j=0

Pn (x, y) se llama Polinomio de Taylor de grado n en dos variables para la funcion f alrededor
de (x0 , y0 ). El u
ltimo termino corresponde al termino de error.
Deducci
on del m
etodo.
Aplicamos (27) para desarrollar las funciones no lineales f y g del sistema (26) en torno del punto
(xk , yk ).
f
f
f (x, y) = f (xk , yk ) + (x xk ) (xk , yk ) + (y yk ) (xk , yk )+
x
y
=

(x xk )2 2 f
(x xk )(y yk ) 2 f
(, ) +
(, )+
2
2
x
2
xy

(y yk )2 2 f
(, ) +
2
y 2
g
g
= g(xk , yk ) + (x xk ) (xk , yk ) + (y yk ) (xk , yk )+
x
y

=
g(x, y)

(x xk )2 2 g
(x xk )(y yk ) 2 g
(, )+
(,
)
+
2
x2
2
xy

(y yk )2 2 g
(, ) +
2
y 2

(28)

Evaluando en (xk+1 , yk+1 ) y trucando el termino de error se tiene:

f
f
(xk , yk ) + (yk+1 yk ) (xk , yk )
x
y

f (xk+1 , yk+1 )

= f (xk , yk ) + (xk+1 xk )

g(xk+1 , yk+1 )

g
g
= g(xk , yk ) + (xk+1 xk ) (xk , yk ) + (yk+1 yk ) (xk , yk ).
x
y

(29)

27

Como la aproximaci
on (xk+1 , yk+1 ) corresponde donde,
f (xk+1 , yk+1 ) = 0 = g(xk+1 , yk+1 ),
reordenando se obtiene:
(xk+1 xk )

f
f
(xk , yk ) + (yk+1 yk ) (xk , yk )
x
y

= f (xk , yk )
(30)

g
g
(xk+1 xk ) (xk , yk ) + (yk+1 yk ) (xk , yk )
x
y

= g(xk , yk ).

Con el fn de simplificar la notacion, escribimos el sistema en forma matricial

Jk = xk .

(31)

donde,

f
(xk )

g
(xk )
y

f
(x )
x k
Jk = g
(xk )
x

es la matriz jacobiana, en

xk =

xk
yk


para

xk =


 


f (xk )
x
xk+1 xk
y=
=
g(xk )
y
yk+1 yk

De (31) podemos escribir

xk+1

= J 1 F(xk )
= xk J 1 F(xk )

(32)

Esta u
ltima ecuaci
on es una generalizacion de la formula de Newton-Raphson para funciones no
lineales de una variable real.
Ejemplo 3.3 Resolver el ejemplo 1.1 i) usando metodo de Newton.Raphson generalizado.
Soluci
on. Consideremos:
f (x, y) = yex 2 = 0
g(x, y) = x2 + y 4 = 0
sabemos que tiene dos races, aproximaremos la raz positiva, usando

x0 =


0,6
3,7

Determinamos:

F(x) =


yex 2
x2 + y 4

28


,

J=

yex
2x

ex
1

(33)

Evaluando en el dato inicial x0 , tenemos



3,7e0,6 2
F(x0 ) =
(0,6)2 + 3,7 4


,

J0 =

3,7e0,6
2(0,6)

e0,6
1

de donde,


2,03060
1,2

  


0,548812
x
0,0306031

=
1
y
0,006

De aqu



 


1
1 0,548812
0,0306031
x

=
2,03060
0,006
y
2,68917 1,2

luego,


 

x
0,000865
=
y
0,05896

Como
x1 x0 = x0

x1 = x0 + x0

entonces,

x1 =



 

0,6
0,000865
0,599135
+
=
3,7
0,05896
3,64104

Si realizamos una segunda iteraci

on, se obtiene:

x2 = x1 + x1 =


0,599125
3,64105

Actividades 3.2 Use el metodo de Newton-Raphson para determinar la aproximaci

on x3 de la
soluci
on del sistema no lineal

x + 12 cosz = ln y1
(34)
y ln 4x + 1 = 2z + ez
4x + y sen z = 1
use como aproximaci
on inicial
1
3

x0 =
3

Ayuda. Programando el algoritmo de Newton Raphson generalizado, en MAPLE se obtiene:

,2399765202
0,2497429663
0,2499999723
x1 = ,3506602354 , x2 = 0,3672616918 , x3 = 0,3678789014
,551268419 101
0,26800671 103
0,1485685 106

0,2500000000
0,2500000000

0,3678794411
x4 = ,3678794411 , x5 =
9
10
,1098834 10
,7324581877 10
29

Nota 3.1 Observese que el metodo converge a

1
4

x=
e

0

30

Mtodo de Mller
Es un mtodo de interpolacin que aproxima la funcin f(x) cuya raz interesa por un
polinomio cuadrtico.
Para determinar dicho polinomio se consideran tres puntos cercanos a la raz de la ecuacin
f(x) = 0.
Sean
A( x0 , f 0 ) , B( x1 , f1 ) y C.( x2 , f 2 ) dichos puntos y sea P(u ) = au 2 + bu + c el polinomio
que interpola a f ( x) en dichos puntos.
Para simplificar la obtencin del polinomio cuadrtico se supone que el origen pasa por el
punto que esta en medio de los otros dos.
Sea
Sea h1 = x1 x0 , h2 = x0 x2

Los coeficientes a, b, c del polinomio P (u ) se determinan resolviendo el sistema

P(0) = f 0

a(0) 2 + b(0) + c = f 0

P(h1 ) = f1 ah12 + bh1 + c = f1

P(h2 ) = f 2

ah2 2 bh2 + c = f 2

( en este desarrollo se considera x2 < x0 < x1 , es decir, x0 es el punto del medio)

Se obtiene,con d =

h2
:
h1

c = f0 ,

a=

d f1 (1 + d ) f 0 + f 2
,
(1 + d )d h12

b=

f1 f 0 ah12
h1

Con los valores determinados de a, b, c se resuelve P (u ) = au 2 + bu + c = 0

eligiendo la raz ms cercana a x0 . Este valor es:

r = x0

2c

b b 2 4ac
El signo del denominador se toma de modo que coincida con el de b : si b > 0 se considera el signo + y
si b < 0 se elige el signo - en el denominador.
De ese modo el denominador toma su mayor valor absoluto.
El valor r calculado es uno de los tres puntos a considerar en el siguiente paso,es decir,
si r queda a la derecha de x 0 se toman los puntos x0 , x1 , r y si r queda a la izquierda de x0 ,
los puntos a considerar son x0 , x2 , r.
En ambos casos los puntos se redefinen de modo que x0 sea el punto del medio de los tres escogidos.

Ejemplo:

y := sin( x )

1
x
2

Su grfica es:

La raz est entre x2 =1.8 y x1 = 2.2. Tomando x0 =2.0 entre x2 y x1

x0 := 2.0
x1 := 2.2
x2 := 1.8

los correspondientes valores de f(x), h1, h2 y d son:

y0 := -.0907025732
y1 := -.2915035962

y2 := .0738476309
h1 := .2

h2 := .2
d := 1.000000000

Resolviendo el sistema se obtienen los coeficientes del polinomio cuadrtico:

a := -.4531352362
b := -.9133780680
c := -.0907025732

Por ser b < 0 se considera el signo (-) en el denominador de la expresin para r, raz del
polinomio ( r es una primera aproximacin a la raz de f(x) )
r := 1.895252107

La raz r del polinomio se encuentra a la izquierda de x0. Los nuevos valores son
x00 := 1.895252107
x11 := 2.0
x22 := 1.8

h11 := .104747893
h22 := .095252107
d1 := .9093462816

y los nuevos valores de la funcin son:

y00 := .0001983067
y11 := -.0907025732

y22 := .0738476309

Los valores de los nuevos coeficientes del polinomio son:

a1 := -.4730106876
b1 := -.8182594122

c1 := .0001983067

Con dichos valores se obtiene una nueva aproximacin a la raz resolviendo el polinomio
cuadrtico:
r1 := 1.895494425

El valor de la funcin en r1 es:

y3 := -.1294 10 -6

y la solucin exacta (obtenida con Maple) de la ecuacin propuesta es:

rexac := 1.895494267

La magnitud del error absoluto en r1 es:

er := .158 10 -6

Ejemplo 2: Calcular la raz en el segundo cuadrante de la ecuacin:

y := ( 1 + x ) sin( x ) 1 = 0

x2 := 2.8
x0 := 2.9
x1 := 3.0

h1 := .1
h2 := .1

d := 1.000000000
y2 := 1.726773367
y0 := -.0669276161

y1 := -.4355199676
a := 71.25543160

b := -10.81146668
c := -.0669276161
r := 2.894043416

y3 := -.0458478202

x22 := 2.894043416
x00 := 2.9
x11 := 3.0

h11 := .1
h22 := .005956584

d1 := .05956584000
y22 := -.0458478202
y00 := -.0669276161
y11 := -.4355199676
a1 := -1.387518017

b1 := -3.547171713
c1 := -.0669276161

Una segunda aproximacin a la raz es r1 y la raz "exacta" es rexac:

r1 := 2.880990772

rexac := 2.880986321

El valor de la funcin en x = r1 es:

y33 := -.0000155441

M
odulo 6
Prof. M.Ang
elica Vega U.

4.

Sistemas de Ecuaciones Lineales.

4.1.

M
etodos directos

Estos metodos consisten en sustituir la matriz A por una sucesion de matrices elementales, de
manera que se obtiene un sistema sea mas facil de resolver.
4.1.1.

M
etodo de Gauss.

El metodo de eliminaci
on gaussiana es uno de los metodos directos clasico para resolver el sistema
Ax = b, consiste de dos etapas.
Etapa 1. Triangularizaci
on.
Para efectos de la operatoria, resulta mas simple, efectuar las operaciones sobre la matriz aumentada A en lugar de A.
El metodo consiste en obtener una sucesion de sistemas equivalentes
A(k) x = b(k) ,

(k = 1, 2, ... , n)

mediante operaciones elementales de fila, tal que para cada k = 2, 3 , . . . , n.

..
(1)
(1)
(1)
(1)
a
a
.
.
.
.
.
.
.
.
.
.
.
.
a
.
a
11
12
1n
1 n+1

..
(2)
(2)
(2)
0
a22 . . .
...
...
...
a2n
.
a2 n+1

..
..
..
..

0
.
0
.
.
.

.
(k)

(k)
(k)
(k)
(k)
.
A = 0
0
...
akk
ak k+1
...
ak n
.
ak n+1

(k)
(k)
(k)
(k)
0
0
. . . ak+1 k ak+1 k+1 . . . ak+1 n .. ak+1 n+1

..
..
..
..
..
..

.
.
...
.
.
...
.
.

..
(k)
(k)
(k)
(k)
0
0
. . . an k
an k+1
. . . an n
. an n+1 .

(35)

donde,

(1)
a1 n+1
.
..

= a(k)

n
n+1

..
.
(k)
an n+1 .

b(k)

con
A(1) = A ,
31

(36)

(k)
A(k) = (aij ) se obtiene de la siguiente manera:

(k)
aij

(k1)
a

ij
= 0

a(k1)

ij

si, i = 1, 2, . . . k 1 ;
si, i = k, k + 1, . . . n ;
(k1)

ai k1
(k1)

ak1 k1

(k1)
ak1 j

sii = k, k + 1 . . . n ;

j = 1, 2, . . . n + 1
j = 1...k 1 ,

(37)

j = k, k + 1 . . . n + 1 ..

El proceso termina cuando A(n) es triangular superior, teniendose el sistema A(n) x = b(n) .
N
otese que en esta etapa s
olo se obtiene una matriz triangular superior. En la etapa siguiente, se
resuelve el sistema.
Etapa 2. Retrosustituci
on.
Resolvemos la n-esima ecuaci
on
(n)
a n+1
.
xn = n(n)
an n
Usando xn se puede despejar xn1 de la n 1-ecuacion
(n1)

xn1 =

(n1)

an1 ,n+1 an1 n xn

;
an1 n1

continuando con este proceso, se tiene que, para cada

i = n 1, n 2, . . . , 2, 1 ,
Pn
aii ,n+1 j=i+1 aiij xj
aii ,n+1 aii n xn aii n1 xn1 aii i+1 xi+1
xi =
=
aii
aii

(38)

Nota 4.1 Denotaremos por :

mi k =

ai k
ak k

al n
umero que produce un cero en la fila i columna k, denominado multiplicador y a los elementos
ai i le llamaremos elementos pivotes.
(k)

Observemos que, ai j , se obtiene de :

(k)

(k1)

ai j = ai j

(k1)

mk1
i k1 ak1 j

Actividades 4.1
i) Ejemplo 4.1 Use eliminaci
on de Gauss para resolver:
3x1 0,1x2 0,2x3
0,1x1 + 7x2 0,3x3
0,3x1 0,2x2 + 10x3
efectuar los c
alculos con 6 cifras significativas.
32

= 7,85 ,
= 19,4 ,
= 71,4 ,

Soluci
on : la matriz aumentada es

3 0,1 0,2 7,85

7
0,3 19,4
A = 0,1
0,3 0,2
10
71,4

luego,los multiplicadores que producen ceros en la primera columna son :

m2 1 =

0,1
a2 1
=
= 0,03333333
a1 1
3

m3 1 =

a3 1
0,3
=
= 0,09999999
a1 1
3

ii) Ejemplo 4.2 Use eliminaci

on de Gauss para resolver los siguientes sistemas: a)
0,0003x1 + 3,0000x2
1,0000x1 + 1,0000x2

=
=

2,0001 ,
1,0000 ,

1,0000x1 + 1,0000x2
0,0003x1 + 3,0000x2

=
=

1,0000 ,
2,0001 ,

b)

c) Compare los resultados de a) y b) al considerar 3, 4 y 5 cifras significativas. A que conclusi

on llega ?
4.1.2.

Eliminaci
on gaussiana con y sin pivoteo.

De la Actividad I.ii), se observa que pueden surgir dificultades en algunos casos, como cuando el
pivote es peque
no en relaci
on a los otros elementos de la columna a la que pertenece el elemento
pivote. Una forma de mejorar las soluciones se usa como estrategia pivotear. El llamado pivoteo
parcial consiste en seleccionar el elemento de mayor valor absoluto ubicado en la misma columna
que el pivote, bajo la diagonal e intercambiar filas. Al procedimiento de buscar el elemento de
mayor valor absoluto tanto en las columnas como en las filas y luego intercambiar se conoce como
pivoteo total.
El pivoteo completo o total no es muy usado puesto que el intercambio de columnas implica un
intercambio de las inc
ognitas, lo que agrega complejidad injustificada.
Actividades 4.2

1.

Use eliminaci
on de Gauss para resolver el sistema :
x1 + x2 x3
6x1 + 2x2 + 2x3
3x1 + 4x2 + x3

usando,
33

= 3 ,
= 2,
= 1,

Gauss simple.
Eliminaci
on gaussiana con pivoteo parcial.
Eliminaci
on gaussiana con pivoteo total.
2.

Para k = 0, 1, 2, ..., sea x(k) = [xk , yk , zk ]T generada por :

4xk+1
6yk+1
4xk+1

= xk yk + 2 ,
= +2xk + yk zk 1 ,
= xk + yk zk + 4 ,

a) Demuestre que la sucesi

on de vectores {x(k) }
k=1 es convergente para cualquier vector
(0)
inicial x .
b) Determine el vector x al cu
al converge usando Gauss con pivote parcial.
c) Encuentre el vector x iterando el sistema dado.
3.

Suponga que trabaja con una m

aquina que redondea a dos decimales para resolver el siguiente
sistema :
0,10 101 x + y
xy

= 1
= 0

Use el metodo de eliminaci

on de Gauss para encontrar la soluci
on.
Efect
ue un pivoteo parcial y compare con la soluci
on anterior. La soluci
on exacta del
sistema es x = y = 100
.
101

34

4.2.

Factorizaci
on LU
A = LU

Una primera interrogante es Existe tal factorizacion? . Un primer resultado que garantiza la
existencia de esta factorizaci
on dice,
Teorema 4.1 Si puede aplicarse la eliminaci
on gaussiana al sistema Ax = b sin intercambio de
filas, entonces la matriz A puede factorizarse como
A = LU ,
donde

U =

a11

(1)

a12

(1)

0
..
.

(2)
a22

0
y

...

..

...
..
.
..
.
0

(1)

a1n
..
.

m21

L =
(n1)
.
an1 n1
..

(n)

ann
mn1

0
1
..

.
...

...
..
.
..
.
mn n1

0
..
.

0
1

(39)

N
otese que al factorizar A en LU ,
Ax = b

L(U x) = b

entonces, es posible resolver el sistema utilizando dos sistemas que involucran matrices triangulares,
se resuelve primero Ly = b y luego U x = y. De modo que si L = (lij ), entonces se obtiene el
vector y por sustituci
on hacia adelante,
b1
y1 =
l11
y los yi para cada i = 2, 3, . . . n, mediante


i1
X
1
yi =
bi
lij yj .
lii
j=1
As conociendo el vector y podemos obtener en el sistema U x = y el vector x por retrosustitucion.
El Teorema 1 proporciona una factorizacion de la matriz A, observe que los elementos de la matriz
L son los elementos multiplicadores mi j que se obtuvieron al realizar la eliminacion gaussiana,
excepto los elementos de la diagonal principal que son 1. Esta factorizacion se conoce como Factorizaci
on de Doolitle.
Queda analizar una segunda interrogante . si A admite una factorizacion LU , esta sera u
nica?

35

Analicemos el problema general, el factorizar A como LU , significa que hay

ciones

u11 u12
l11 0
0 ... 0
a11 a12 . . . a1n
l21 l22 0 . . . 0 0 u22

a21 a22 . . . a2n
l31 l32 l33 . . . 0 0

0
=

..

..
.
.
..
.. .
.
..
..
..
.. ..
.
.
.
.
.
.
.
.
.
. .
an1 an2 . . . ann
0
0
ln1 ln2 ln3 . . . lnn

que resolver las ecuau13

u23
u33
..
.

...
...
...
..
.

u1n
u2n
u3n
..
.

...

unn .

(40)
Notemos que las ecuaciones anteriores, no determinan a L y a U en forma u
nica. Para cada i se
puede asignar un valor no nulo lii o uii (pero no a ambos). Porque ? Justifique . Por ejemplo, si
lii = 1 para i = 1, 2, . . . , n, es la descomposicion de Doolitle que se obtiene como vimos al aplicar
eliminaci
on gaussiana. Otra eleccion simple es uii = 1 para cada i en este caso la factorizacion es
la de Crout, o bien asumir que lii = uii en cual caso tenemos la factorizacion de Choleski
Actividades 4.3
i) Descomponer la matriz A usando la factorizaci
on de Crout, en el sistema :
x1 + x2 x3
6x1 + 2x2 + 2x3
3x1 + 4x2 + x3

= 3 ,
= 2,
= 1,

ii) Problema N 3 de la PEP1 (Sem II-2000) Para determinar la temperatura en cuatro puntos
interiores equidistantes de una barra, hay que resolver la ecuaci
on :
Atk+1 = tk ,

k = 0, 1, 2, ...

donde,

10
3 1 0
0
1 3 1 0 0 12

A=
0 1 3 1 , t = 12
10
0
0 1 3

(41)

(tk es la distribuci
on de la temperatura en los puntos interiores de la barra en el instante k.
a) Mediante el metodo de Gauss (sin pivote) determine t1 .
b) Obtenga t2 y t3 , usando una descomposici
on A = LU .
iii) Resolver el sistema

x1
2x1

3x1
x1

+ 2x2
4x2
+ 8x2
+ 2x2

3x3
5x3
8x3
6x3

+
+

x4
x4
x4
4x4

= 1
= 0

= 2
= 1

(42)

iv) Investigue c
omo a partir de la factorizaci
on de Doolitle puede obtener la factorizaci
on de
Crout ? Ilustre el metodo con un ejemplo.
36

4.2.1.

Matrices definidas positivas.

Definici
on 4.1 Una matriz A simetrica es definida positiva si y s
olo si
xT Ax > 0

, x

Rn {0}

Ejemplo 4.3 Verifique si la matriz dada es definida positiva.

2 1 0
1 2 1
0 1 2

(43)

(44)

En la pr
actica puede resultar difcil determinar si una matriz es definida positiva a partir de la
definici
on, el siguiente resultado constituye un criterio de gran utilidad para este proposito.
Teorema 4.2 A es simetrica y definida positiva,
el determinante:

 a11


Dk =  ...

 ak1

si y s
olo si Dk > 0 , k = 1, ..., n, donde Dk es

..
.

a1k
..
.
akk

Dem. Consultar, Introduction to Numerical Analysis de Stoer - Bulirsch.

Teorema 4.3 Si A es simetrica y definida positiva, entonces :
i) A es no singular
ii) aii > 0, para cada i = 1, 2, ..., n.
Dem.
i) Supongamos A no invertible, entonces exite x 6= 0 tal que el sistema Ax = 0, de donde
xT Ax = 0, lo que contradice el hecho de A definida positiva. Por tanto, A es no singular.
ii) Ejercicio
Teorema 4.4 Si A es una matriz real, simetrica y definida positiva, entonces A = LLT , donde L
es una matriz triangular inferior con diagonal positiva.
Dem. Consultar An
alisis Numerico de Burden - Faires o Kincaid - Cheney
4.2.2.

Matrices diagonal dominante.

Definici
on 4.2 Sea A una matriz cuadrada de orden n, diremos que A es diagonal dominante
si
n
X
|aii |
|aij |, i, 1 i n
(45)
j=1
j6=i

Veamos algunos resultados importantes,

37

Teorema 4.5 Si A es diagonal dominante, entonces A es invertible.

Teorema 4.6 Si A es diagonal dominante y todos los elementos de la diagonal son reales positivos

, entonces la parte real de cada valor propio es positiva. Dem. Algebra

Lineal. Grossman. Mc Graw
Hill.
Teorema 4.7 Si A es simetrica, diagonal dominante, con los coeficientes diagonales positivos,
entonces A es definida positiva.
Dem. Usando el teorema anterior, si i > 0 deduciremos que xt Ax > 0. Como A es simetrica,
tiene un conjunto de autovectores ortonormales (teorema espectral), entonces cualquier vector x
puede escribirse como una combinaci
on lineal de ellos,
x = c1 x1 + c2 x2 + ... + cn xn ,
entonces
Ax = c1 Ax1 + c2 Ax2 + ... + cn Axn = c1 1 x1 + c2 2 x2 + ... + cn n xn .
De la ortogonalidad y la normalizaci
on xti xi = 1. Luego,
xt Ax = (c1 xt1 + c2 xt2 + ... + cn xtn )(c1 1 x1 + c2 2 x2 + ... + cn n xn )
= c21 1 x1 + c22 2 x2 + ... + c2n n xn
como i > 0, se tiene que i xt Ax > 0.
Observaci
on 4.1 El recproco del teorema anterior es falsa, es decir existen matrices simetricas
definidas positivas que no son diagonal dominante. De un contraejemplo.

38

Modulo 7
Prof : Mara Angelica Vega U.

5.

M
etodos iterativos.

5.1.
5.1.1.

Normas vectoriales y matriciales.

Preliminares.

Actividades 5.1
i) Si V es un espacio vectorial, Que es una norma sobre V ? , C
omo se define la distancia
entre dos vectores ? Cu
ando una sucesi
on de vectores es convergente ?
ii) Sea V = IRn , las normas de uso frecuente en IRn son,
Pn
kxk1 = i=1 |xi | , llamada norma 1
kxk2

1/2

llamada norma eucldea o norma 2

kxk

= m
ax1in |xi |

llamada norma infinito o del m

aximo.

Pn

i=1

x2i

Ejemplo 5.1 Compare la longitud de los siguientes tres vectores en IR4 , usando cada una de las
normas anteriores.
x = (4, 4, 4, 4)T

y = (0, 5, 5, 5)T

w = (6, 0, 0, 0)T

Definici
on 5.1 Diremos que dos normas son equivalentes si existen dos n
umeros y positivos
tales que
kvkq kvkp kvkq .
Demuestre las siguientes afirmaciones.(Haga una en clases ).
iv) Las normas kvk y kvk1 son equivalentes, en efecto
kvk kvk1 nkvk
v) Las normas kvk y kvk2 son equivalentes
kvk kvk2

nkvk

1
kvk2 kvk kvk2 .
n

vi) Las normas kvk1 y kvk2 son equivalentes

kvk2 kvk1

nkvk2 .

Observaci
on 5.1 Por que es importante que las normas estudiadas sean equivalentes ?
Porque se puede utilizar cualquiera de ellas. As si, por ejemplo, se tiene una sucesion de vectores
que tiende a v para una cierta norma, entonces la sucesion tambien converge a v para las otras
normas (equivalentes).
40

5.1.2.

Normas matriciales

Definici
on 5.2 Sea An el conjunto de las matrices de orden n, con las operaciones habituales,
definidas sobre IR (o C
I ). Definiremos una norma matricial como una aplicaci
on
k k : An IR
que verifica las siguientes propiedades :

i)
ii)
iii)
iv)

kAk 0 , A An y kAk = 0 A = 0
kAk = ||kAk , IR y A An
kA + Bk kAk + kBk , A, B An
kABk kAkkBk , A, B An .

Comentarios:
i) An puede ser considerado como un espacio vectorial de dimension n2 , As las cosas, las tres
primeras propiedades son similares a las correspondientes de una norma vectorial; mientras
que la u
ltima propiedad es propia de las normas matriciales.
ii) Un concepto importante para nuetros propositos es norma matricial inducida, puesto que
permite relacionar las normas vectoriales con las normas matriciales, pilar fundamental para
el estudio de convergencia y analisis de error en el calculo numerico.
Definici
on 5.3 Sean A una matriz n n y k kv una norma vectorial definida en IRn , llamaremos
norma matricial inducida o norma matricial subordinada, con respecto a la norma vectorial dada
al n
umero:
kAvkv
= sup kAvk.
kAk = sup
v6=0 kvkv
kvk=1
Actividades 5.2 Investigue, como demostrar los siguientes resultados importantes.
i) Para toda norma matricial inducida, se cumple
kAxkv kAkkxkv ,

x IRn .

En este caso se dice que la norma matricial y la norma vectorial son compatibles.
ii) Sea A = (aij ) una matriz cuadrada real, entonces
kAk1
kAk2
kAk

n
X
kAvk1
= max
|aij |
1jn
kvk1
i=1
q
q
kAvk2
= sup
= (AT A) = (AAT ) = kAT k2 (norma espectral)
kvk2
n
X
kAvk
= sup
= max
|aij | (norma del m
aximo.)
1in
kvk
j=1

sup

41

(46)

Nota 5.1 Si A es una matriz de orden n, (A) denota el radio espectral de A, es decir
(A) = max |i |
i (A)

Observaci
on 5.2 Existen normas matriciales que no son subordinadas a ninguna norma vectorial,
por ejemplo, la norma

1/2
X
1/2
kAkE =
|aij |2
= (tr(A A))
i,j

es una norma matricial no subordinada.

Nota 5.2 La matriz aumentada de un sistema lineal Ax = b en general se almacena con un error
de redondeo, si la soluci
on calculada se obtiene usando una estrategia de pivoteo que limita el error,
entonces la soluci
on calculada y la soluci
on exacta son iguales hasta alrededor de la exactitud de
la m
aquina usada. Por lo tanto es u
til saber si la matriz de los coeficientes est
a mal condicionada.
Se entiende por matriz mal condicionada si peque
nos cambios en los elementos de la matriz
pueden producir grandes cambios en el vector soluci
A,
on.
Un indicador de mal condicionamiento es el llamado n
umero de condici
on de la matriz A que
se define como
K(A) = kAk kA1 k
Este n
umero tiene la propiedad que, K(A) 1, luego un n
umero de condicion de una matriz
cercano a 1 es indicador de buen condicionamiento y un n
umero de condicion muy grande, es
indicador de mal condicionamiento.
Actividades 5.3
i) Determine kxk , kAk y kA1 k para las k k1 y k k , si

2
x = 5
3

2
A=
3

1 4
2 0
3 2

ii) Cu
al es el n
umero de condici
on de la matriz A.
iii) Determine la norma espectral de A. (Para entretenerse fuera de la clase ).
5.1.3.

Otros resultados importantes para el estudio de los m

etodos iterativos.

Si A es una matriz real de orden n, se tiene

i) (A) kAk, para cualquier norma k k.
42

(47)

ii) Diremos que A es una matriz convergente, si

lm (ak )ij = 0

para cada i = 1, 2, ..., n

j = 1, 2, ..., n.

Las siguientes afirmaciones son equivalentes:

i)
ii)

A es una matriz convergente.

lm k(An )k = 0, para alguna norma natural k k.
n

iii) (A) < 1.

iv) lm An x = 0.
n

Actividades 5.4 ! Entretengase ! , demuestre que A es una matriz convergente si

 1

0
2
A=
1
1
4

(48)

Teorema 5.1 Las siguientes afirmaciones son equivalentes:

i)
ii)

A es una matriz convergente.

lm k(An )k = 0, para alguna norma natural k k.
n

iii) (A) < 1.

iv) lm An x = 0.
n

5.2.

M
etodos iterativos.

Los metodos directos requieren alrededor de 2n3 /3 operaciones aritmeticas para resolver un sistema
lineal de n n, lo que limita el tama
no de los sistemas que pueden ser resueltos de manera directa.
Por ejemplo cuando se resuelven numericamente ecuaciones diferenciales pueden surgir sistemas
de 20000 variable.
Una alternativa que puede lograr gran exactitud para n grande son los metodos iterativos.
Para obtener un metodo iterativo procedemos como antes, re-escribimos el sistema de ecuaciones
como una iteraci
on de punto fijo. Es decir,
Ax = b.

(49)

lo escribimos, en la forma
x(k+1) = M x(k) + c

para cadak 1.

(50)

donde la matriz M , se llama matriz de iteraci

on. Relacionando nuestros conocimientos sobre
normas con estos metodos, existen algunos resultados importantes que garantizan la convergencia
de estos metodos. ( Se sugiere cosultar la demostracion en los textos de la bibliografa entregada.)
43

1) El radio espectral de la matriz M , (M ) cumple

(M ) kM k
para cualquier norma matricial inducida.
2) Sea M matriz n n, el proceso iterativo (61) converge para cada k 1 y c 6= 0 a la u
nica
soluci
on si y s
olo si (M ) < 1.
Para transformar (49) en un problema de punto fijo, se descompone la matriz A en la suma,
A=L+D+U
donde, D es la diagonal principal de A, L es la parte diagonal inferior de la matriz A y U es la
parte triangular superior de A con ceros en las diagonales.
5.2.1.

M
etodo de Jacobi.

i) Forma matricial
Transformamos el sistema,
Ax = b (D + L + U )x

Dx

=b
= b (L + U )x
= D1 [b (L + U )x]
= D1 (L + U )x + D1 b
= MJAC x + c.

(51)

De aqu que la matriz de de iteracion de Jacobi es,

MJAC = D1 (L + U ).
y c = D1 b.
ii) Forma en terminos de las componentes
La i-esima componente del vector correspondiente a (k + 1)-esima iteracion es,

(k+1)

xi

1
bi
aii

44

n
X
j=1
j6=i

(k)
aij xj .

(52)

5.3.

M
etodo de Gauss-Seidel.

i) Forma matricial
En este caso se hace la transformacion,
Ax = b (D + L + U )x

(D + L)x

=b
= b Ux
= (D + L)1 [b U x]
= (D + L)1 U x + (D + L)1 b
= MGS x + c.

(53)

De aqu que la matriz de de iteracion de Gauss-Seidel es,

MGS = (D + L)1 U.
y el vector c = (D + L)1 b.
ii) Forma en terminos de las componentes
Este metodo es un mejoramiento del metodo de Jacobi, que consiste en calcular la compo(k)

nente xi

a partir de las componentes de x(k1) en la forma siguiente, para i > 1, han

sido calculadas las

(k1)

(k)

(k)

x1 , ..., xi1 que supuestamente son mejores aproximaciones que las

(k1)

x1
, ..., xi1 a las componentes x1 , ..., xi1 , de la solucion real , por lo tanto se usa esta
idea al resolver la i-esima ecuacion,

(k)

xi

Pi1 
j=1

(k)

aij xj

Pn

j=i+1

(k1)

aij xj

+ bi
para

aii

i = 1, 2, ..., n.

Actividades 5.5
i) Dado el sistema lineal:

x1
5x1
x1

0x2
+ x2
+ 6x2

+ 3x3
+ 2x3
+ 2x3

=
2
= 5
= 11

(54)

a) Es Definida positiva la matriz A. Justifique.

b) Determine la soluci
on exacta del sistema.
c) Comenzando con x(0) = 0 determine cuatro iteraciones usando el metodo de Jacobi y
cuatro iteraciones usando Gauss-Seidel
ii) El metodo iterativo de Jacobi puede expresarse como :
x(k+1) = MJAC x(k) + c.
donde = MJAC es la correspondiente matriz de iteraci
on.
45

(55)

a) Verifique que (55) puede escribirse como:

x(k+1) = x(k) + z(k) .

(56)

donde,
z(k) = c BJAC x(k)
(z

(k)

con

BJAC = I MJAC

se llama vector de correcci

on)

b) Para obtener un metodo de relajaci

on, se introduce en (56) el coeficiente para obtener
:
x(k+1) = x(k) + z(k) .
Use el metodo (57) con

3
1

0
0

= 0,5 para estimar la aproximaci

on x(2) del sistema :

x1
10
1 0
0

12

3 1 0

x2

=
12
1 3 1 x3
10
0 1 3
x4

(57)

(58)

con x(0) = 0.

5.4.

M
etodo de sobrerelajaci
on

Estos metodos se obtienen modificando el metodo de Jacobi o el de Gauss-Seidel, introduciendo

un par
ametro llamado par
ametro de relajaci
on.
i) Forma matricial
En efecto, la forma matricial del esquema, se obtiene considerando la siguiente descomposicion
para la matriz A.


(1 )
D
A=
+L
D + U.

(59)

donde :
es un par
ametro
D,L,U son las matrices diagonal y triangulares inferior y superior definidas anteriormente.
sustituyendo (59) en Ax = b, se tiene :


D
(1 )
x + Lx +
D + U x = b.

multiplicando (60) por :

Dx + Lx + [(1 )D + U ]x = b
46

(60)

y de aqu :
x = (D + L)1 [(1 )D U ]x + (D + L)1 b
que en forma recursiva puede expresarse en la forma :
x(k) = M x(k1) + c

(1)

con :
M = (D + L)1 [(1 )D U ]
c = (D + L)1 b
ii) En termino de componentes
Para fines computacionales la ecuacion matricial, puede escribirse:

i1
n
X
X

(k)
(k1)
(k)
(k1)
bi
,
aij xj
aij xj
xi = (1 )xi
+
aii
j=1
j=i+1

(2)

para i = 1, 2, ..., n.
Si = 1 , tenemos el metodo de Gauss-Seidel.
Si 0 < < 1, los metodos de relajacion se denominan metodos de sub-relajaci
on y se pueden usar
en sistemas que no son convergentes por el metodo de Gauss-Seidel.
Si > 1 los metodos se denominan metodos de sobre-relajaci
on y se pueden usar para acelerar la
convergencia en sistemas que son convergentes mediante el metodo de Gauss-Seidel.

5.5.

Resultados importantes

i) Para cualquier x(0) Rn , la sucesion definida por

x(k+1) = M x(k) + c para cadak 1.

(61)

converge a la soluci
on u
nica si y solo si (M ) < 1.
ii) Si kM k < 1 para cualquier norma matricial subordinada, la sucesion generada por (61)
converge para cualquier x(0) Rn y se satisafacen las cotas de error:
kx x(k) k kM kk kx(0) xk
kx x(k) k

kM kk
kx(1) x(0) k
1 kM k

iii) Si A es diagonal dominante, entonces para cualquier eleccion de x(0) los metodos de Jacobi
y Gauss-Seidel convergen a la solucion exacta del sistema.
iv) Si A es definida positiva y 0 < < 2, entonces el metodo SOR converge para cualquier
aproximaci
on inicial.
v) Si A es definida positiva y tridiagonal, entonces (MGS ) = (MJAC )2 < 1 y la eleccion
optima de , es

2
p
=
.
1 + 1 (MJAC )2
47

5.6.

An
alisis de error y N
umero de Condici
on

Consideremos el sistema A
x = b y supongamos que A es una matriz invertible, analicemos las
siguientes situaciones:
1.

Si se perturba la matriz A1 para obtener una nueva matriz B la solucion x = A1 b resulta

e = Bb. Nos preguntamos Que tan grande es la solucion
perturbada y la escribiremos x
perturbada en terminos absolutos y relativos ? sol. Consideremos || ||p una norma vectorial
y su correspondiente norma matricial subordinada, determinamos
i) el error absoluto
e|| = ||x Bb|| = ||x BAx|| = ||(I BA)x|| ||I BA|| ||x||
||x x
ii) el error relativo
e||
||x x
||I BA||
||x||
representan las medidas de los errores en terminos absolutos y relativos.

2.

e Si x y x
e satisfacen
Consideremos ahora que se perturba el vector b para obtener el vector b.

e
e en terminos absolutos y
respectivamente Ax = b y Ae
x = b. En cuanto difiere x y x
relativos ?
i) En terminos absolutos
e k=k A1 (b b)
e kk A1 kk b b
ek
e|| =k A1 b A1 b
||x x
ii) En terminos relativos
e||
||x x

e = ||A1 || ||b b||

e
||A1 || ||b b||

||A1 || ||Ax||
e||
||x x
||x||

||b||
||b||

e
||b b||
||b||

||A1 || ||A||

e
||b b||
||b||

es decir, el error relativo de la solucion x esta acotado por (A) veces el error relativo de b .
Observaci
on 5.3
i) El n
umero de condici
on depende de la norma vectorial seleccionada.

48

ii) Si (A) es peque

no, entonces peque
nas perturbaciones en b conducen a peque
nas perturbaciones en x.


1
1+
Ejemplo 5.2 Analizar el condicionamiento de la matriz A =
con > 0.
1
1
Sol. Determinamos A1
A1 =

1
1 +

1
1

Usando norma infinito, tenemos,

||A|| = 2 +

(A) =

(2 + )

||A1 || = 2 (2 + )
2
=

4 + 4 + 2
4
> 2
2

Si 0,01, entonces (A) 40000. Esto significa que una peque

na perturbaci
on sobre b puede
inducir una perturbaci
on o error relativo de 40000 veces mayor en la soluci
on del sistema Ax = b.
Observaci
on 5.4 En el caso de los metodos iterativos, podemos medir que tan buena es la soluci
on
e, calculando
aproximada x
i) Ae
x
ii) y r = b Ae
x. Este vector se denomina vector residual.
e se llama vector residual.
Observamos adem
as que la diferencia e = x x
Podemos resumir una relaci
on importante entre vector error y vector residual.
e donde b
e = b r implica r = b.
e
e es la soluci
Teorema 5.2 Si x
on exacta de Ae
x = b,
1 ||r||
||e||
||r||

(A)
A) ||b||
||x||
||b||
Criterios de parada.
Se define el vector el residual r(k) = b Ax(k) y el vector error e(k) = x x(k) .
Criterio 1. Dada una tolerancia , iterar hasta que,
kr(k) k
< .
kr(0) k
Criterio 2. Sean, b , x , x(0) Rn y A una matriz no singular, entonces:
kr(k) k
ke(k) k

(A)
.
ke(0) k
kr(0) k
donde, (A) es el n
umero de condicion de la matriz A.
49

Criterio 3. El criterio 1, depende del valor inicial x(0) , de modo que si la aproximacion inicial
no es buena, los resultados pueden ser grotescos. Se recomienda usar en su lugar,
kr(k) k
< .
kbk
Actividades 5.6
i) Demuestre el teorema anterior
ii) Demuestre que si ||A|| < 1 entonces
I A es invertible.
||(I + A)1 || (1 ||A||)1
(I + A)1 = I A + A2 A3 +
1
||(I + A)1 ||
1 + ||A||
Teorema 5.3 El radio espectral
(A) = inf ||A||
||.||

es el nfimo sobre todas las normas matriciales subordinadas.

Actividades 5.7
1) Contra los deseos de su corredor Jeanette compra x ttulos de acciones de ENDESA, y ttulos
de acciones de CTC Chile y z ttulos de acciones de Chilgener . A inicios de Enero, Febrero
y Marzo, las acciones de ENDESA valan $5,18 , $3,14 y $3,01 por ttulo, las acciones de
CTC Chile $2,10 , $6,03 y $2,18 por ttulo y las acciones de Chilgener, $2,90 , $2,80 y $5,53
por ttulo, respectivamente.
a) Determine el vector (x(2) , y (2) , z (2) ) mediante el metodo de Gauss-Seidel, si Jeanette
gast
o en total $2945,30 , $4113,85 y $3417,00 en Enero, Febrero y Marzo, respectivamente. (x(0) = (150, 150, 150)T ).
(Justifique).
b) Determine una estimaci
on del error absoluto de la aproximaci
on anterior.
c) De acuerdo a lo obtenido en a), Cu
antos ttulos de cada empresa adquiri
o Jeanette ?
2) El sistema de ecuaciones Ax = b, donde






3 2
x1
1
A=
,x =
,b =
1 2
x2
2

(62)

puede resolverse mediante el siguiente metodo iterativo :

x(k+1) = x(k) + (Ax(k) b),


1
x(0) =
,
1
50

(63)

a) Determinar todos los valores de para los cuales se puede garantizar la convergencia
del algoritmo.
b) Determinar si para = 0,4 el algoritmo es convergente.

51

2
Aproximaci
on de funciones.
M
odulo 10
Prof. Mara Angelica Vega U.

El prop
osito de estos m
odulos es entregar los conceptos te
oricos b
asicos
necesarios y algunas aplicaciones en forma dirigida, de modo que el estudiante pueda participar en forma activa, en el proceso ense
nanza - aprendizaje,
elaborando su propio material de estudio. El logro deseable es que el estudiante aprenda a aprender. Es claro que este material debe ser complementado
con el apoyo de la bibliografa sugerida.

2.1

Interpolaci
on Polin
omica.

El Problema :
Dado el conjunto de (n + 1) puntos, del grafico de una funcion f , denotados por
(x0 , f (x0 )) , (x1 , f (x1 )) . . . (xx , f (xn )),
se requiere encontrar una funcion polinomica P (x) al que llamaremos polinomio interpolante para f , de a lo mas grado n que satisfaga la condicion
de interpolacion, expresada por
P (xi ) = f (xi ) ,

i = 0, 1, ...n.

(2.1)

Soluciones al Problema :
Inspirandonos en el criterio de Taylor, de acuerdo a (2.1), observemos
que el polinomio de Taylor es combinacion lineal de los elementos de la base
{1, (x x0 ), (x x0 )2 , (x x0 )3 ...(x x0 )n ...}
para el espacio vectorial de los polinomios de una variable real sobre IR. Es
obvio que los polinomios de Taylor son u
tiles para aproximar una funcion
sobre intervalos peque
nos alrededor del punto x0 , pero en la practica este no
siempre es el caso y es necesario usar un metodo mas eficiente que incluya
informacion en otros puntos.

 POLINOMICA.

2.1. INTERPOLACION

2.1.1

M
etodo de los coeficientes indeterminados.

Tomando en cuenta lo anterior, podemos proponer que el polinomio interpolante sea una combinacion lineal de los elementos de la base canonica del
espacio de los polinomios por la sencillez en los calculos, es decir
Pn (x) =

n
X

ai xi = a0 + a1 x + a2 x2 + ... + an xn .

(2.2)

i=0

As el problema se reduce a encontrar los coeficientes de las potencias de

x, de donde se obtiene un sistema de ecuaciones cuyas incognitas son estos
coeficientes. Veamos estas ideas aplicadas a un problema practico.
Actividad I.
Suponga que una empresa contrata servicios y entrega la siguiente informacion exacta,

Piezas Vendidas (en millones)

Ganancia Neta (en millones)

1
11.2

2
15.3

3
17.1

4
16.9

5
15.0

Tabla 1.
se pide determinar cuantas piezas deben venderse para maximizar la
ganancia neta, dejando estipulado en las bases del contrato que no habra
futuras contrataciones si la respuesta es equvoca . Se debe hacer una recomendacion y justifcarla.

Consideraciones sobre el m
etodo.

Este metodo que resulta muy simple en este caso, presenta la dificultad
que la matriz de los coeficientes asociada al sistema de ecuaciones que se
obtiene al aplicar la condicion de interpolacion, es la matriz de Vandermonde,
que se caracteriza por volverse rapidamente mal condicionada cuando n, el
orden de la matriz crece en funcion del grado del polinomio, por tanto este
metodo no es conveniente para hallar polinomios interpolantes.
De aqu que es importante estudiar otra forma de escribir este polinomio,
de modo de obviar el que tengamos que resolver un sistema de ecuaciones
lineales mal condicionado. Los metodos de Lagrange y de Newton son los
que tradicionalmente se estudian en cualquier curso de calculo numerico,
ambos se basan en escribir el polinomio de interpolacion como combinacion
lineal de elementos de alguna base del espacio vectorial de los polinomios
distinta de la canonica.

2.1.2

M
etodo de Lagrange.

Sea Pn (x) el polinomio que interpola a la funcion f en los (n + 1) puntos

(x0 , f (x0 )) , (x1 , f (x1 )) , (x2 , f (x2 )) , , (xn , f (xn ))
y supongamos que seg
un una cierta base
{l0 (x), l1 (x), ..., ln (x)},
para el espacio vectorial de los polinomios, el polinomio interpolante se escribe :
Pn (x) =

n
X

f (xi )li (x) = f (x0 )l0 (x) + ... + f (xn )ln (x).

(2.3)

i=0

Si resolvemos de esta forma el problema, nuestro trabajo consistira en encontrar los polinomios li (x) de modo que se verifique
f (xi ) = Pn (xi ) = f (x0 )l0 (xi ) + + f (xi )li (xi ) + + f (xn )ln (xi ). (2.4)
Observemos que para que se verifique la condicion (2.4) para cada i, es
suficiente que el polinomio li (x) verifique que:
li (xj ) = 0

0 j n

i 6= j

li (xi ) = 1

(2.5)
(2.6)

De la condicion (2.5) se deduce que las races del polinomio li (x) son los
x = xj

0jn

i 6= j ,

por lo tanto podramos escribir este polinomio como :

li (x) = (x x0 )...((x xi1 )(x xi+1 )...(x xn ).

(2.7)

Pero este polinomio no cumple con la condicion (2.6), puesto que si lo evaluamos en el nodo xi vale:
li (xi ) = (xi x0 )...((xi xi1 )(xi xi+1 )...(xi xn ).

(2.8)

Luego si dividimos el polinomio li (x) por el valor li (xi ), obtenemos el polinomio :

li (x) =

n
Y

(x xj )
(x x0 )...((x xi1 )(x xi+1 )...(x xn )
=
.
(xi xj )
(xi x0 )...(xi xi1 )(xi xi+1 )...(xi xn )
j=0 j6=i

(2.9)
Que por construccion cumple las condiciones (2.5)y (2.6).
Los polinomios li (x) se conocen con el nombre polinomios de Lagrange,
mientras que Pn (x) recibe el nombre de polinomio interpolante de Lagrange
para f (x) en los puntos (x0 , f0 ) , (x1 , f1 ) , ..., (xn , fn ), que suele denotarse
generalmente Ln (x).
Aplicaremos la idea de Lagrange al problema del ejemplo anterior con el
fn de obtener alg
un tipo comparacion entre ambos metodos.
Actividad II.
Solucione el problema planteado en la actividad I, usando un polinomio
de Lagrange para interpolar la funcion representada por los datos dados.

 POLINOMICA.

2.1. INTERPOLACION

n 2.1
Observacio
Sabemos que si el polinomio interpolante existe, es
u
nico, por lo tanto es obvio pensar que el polinomio es el mismo usando
cualesquiera de los dos metodos anteriores, m
as a
un usando cualquier
otro metodo.
El metodo de Lagrange para determinar el polinomio interpolante comparado con el metodo de coeficientes indeterminados, es un tanto laborioso por la cantidad de c
alculos involucrados al calcular los elementos
de la base, pero tiene un valor te
orico importante en la deducci
on de
otros metodos numericos, como por ejemplo f
ormulas de integraci
on
numerica, aproximaci
on de soluciones de ecuaciones diferenciales, etc.

Si en la formula (2.3), n = 1 el polinomio interpolante adopta la forma,

P1 (x) = f (x0 )

(x x1 )
(x x0 )
+ f (x1 )
.
(x0 x1 )
(x1 x0 )

(2.10)

o (2.10) en forma simplificada,

P1 (x) = f (x0 )l0 (x) + f (x1 )l1 (x)

con,
l0 (x) =

(2.11)

x x0
x x1
, l1 (x) =
.
(x0 x1 )
(x1 x0 )

cuya grafica es la recta que pasa por los puntos dados por lo que este problema se conoce como Interpolaci
on lineal.
Otra situacion particular es el caso en que n = 2, es decir se tienen los
puntos,
(x0 , f (x0 )) , (x1 , f (x1 )) , (x2 , f (x2 )),
el polinomio P2 (x) que interpola a la funcion f (x) en los tres puntos es

P2 (x) = f (x0 )

(x x1 )(x x2 )
(x x0 )(x x2 )
+ f (x1 )
(x0 x1 )(x0 x2 )
(x1 x0 )(x1 x2 )
(2.12)

(x x0 )(x x1 )
+f (x2 )
.
(x2 x0 )(x2 x1 )
.

Simplificando la notacion igual que antes,

P2 (x) = f (x0 )l0 (x) + f (x1 )l1 (x) + f (x2 )l2 (x)

(2.13)

6
con,
l0 (x) =

(x x1 )(x x2 )
(x x0 )(x x2 )
, l1 (x) =
.
(x0 x1 )(x0 x2 )
(x1 x0 )(x1 x2 )

La grafica corresponde a la parabola que pasa por los puntos dados, en

este caso hablamos de interpolaci
on cuadr
atica.

Consideraciones sobre el m
etodo.
Uno de los inconvenientes que tiene la forma de Lagrange para el polinomio de interpolacion, es que al agregar otro punto de interpolacion o
nodo , los elementos de la base polinomios de Lagrange,
li (x) , para i = 0, 1, ..., n + 1.
deben ser calculados nuevamente uno a uno. Surge entonces la inquietud
de como construir el polinomio de interpolacion en forma mas eficiente, de
modo de aprovechar la informacion ya procesada al agregar nuevos nodos.
Este hecho amerita que busquemos otra forma de escribir el polinomio de
interpolacion, es decir expresarlo en otra base del espacio vectorial.

 POLINOMICA.

2.1. INTERPOLACION

2.1.3

Forma de Newton para el Polinomio de Interpolaci

on

En la forma de Newton el polinomio interpolante Pn (x) es escrito usando una

base parecida a la usada para expresar el polinomio de Taylor, observemos
ambas bases :
{1 , (x x0 ) , (x x0 )2 , , (x x0 )n }.

(2.14)

base usada para definir el polinomio de Taylor y

{1 , (xx0 ) , (xx0 )(xx1 ) , , (xx0 )(xx1 ) (xxn1 )}. (2.15)
base usada para definir el polinomio de Newton. De aqu podemos expresar
el polinomio en la forma,
Pn (x) = a0 +a1 (xx0 )+a2 (xx0 )(xx1 )+ an (xx0 )(xx1 ) (xxn1 ).
(2.16)
Notemos que por construccion este polinomio es de grado a lo mas n y
es logico pensar los coeficientes ai , para i = 0, 1, ..., n deben determinarse
de modo que se verifique la condicion de interpolacion. Para ello evaluamos
Pn (x) en cada xi .
Para x0 obtenemos:
Pn (x0 ) = f (x0 ) = a0 .
Analogamente si evaluamos Pn (x) en x1 tenemos:
a1 =

f (x1 ) f (x0 )
x1 x0

Pn (x) evaluado en x2
a2 =

f (x1 ) f (x0 )
f (x2 ) f (x0 )

.
(x2 x0 )(x2 x1 ) (x2 x1 )(x1 x0 )

As obtenemos sucesivamente los otros valores de las constantes a3 , , an .

Para simplificar la escritura de los coeficientes ai , introducimos el concepto de diferencias divididas lo que permite escribir el polinomio de Newton
en una forma mas simple.
Diferencias Divididas
Consideremos (n+1) puntos cualesquiera Pi (xi , f (xi )) , i = 0...n, no
necesariamente equiespaciados, definimos :
1. diferencias divididas de orden cero por:
f [xi ] = f (xi ).
2. Si consideramos un par de puntos indexados cualesquiera
Pi (xi , f (xi ))

y Pi+1 (xi+1 , f (xi+1 ))

8
definimos la primera diferencia dividida o diferencia dividida de orden
1 como el n
umero :
f [xi , xi+1 ] =

f (xi+1 ) f (xi )
.
xi+1 xi

Al relacionar f [xi , xi+1 ] con su significado geometrico, sabemos que

corresponde a la pendiente del segmento lineal que une Pi y Pi+1 ,
entonces si f [xi , xi+1 ] > 0 implica que f (x) es creciente en [xi , xi+1 ]
y si f [xi , xi+1 ] < 0 implica que f (x) es decreciente en [xi , xi+1 ]. Es
decir, este concepto nos proporciona informacion adicional acerca de
la grafica de la funcion.
3. Consideremos ahora un tro de puntos consecutivos
Pi (xi , f (xi )) , Pi+1 (xi+1 , f (xi+1 ))

y Pi+2 (xi+2 , f (xi+2 ))

definimos la segunda diferencia dividida o diferencia dividida de orden

2 como el n
umero :
f [xi , xi+1 , xi+2 ] =

f [xi+1 , xi+2 ] f [xi , xi+1 ]

.
xi+2 xi

Como f [xi , xi+1 , xi+2 ] es el cuociente entre el cambio de pendiente y

el cambio en x, de modo que si f [xi , xi+1 , xi+2 ] > 0 geometricamente
la grafica de f (x) es concava hacia arriba y si f [xi , xi+1 , xi+2 ] < 0 la
grafica de f (x) es concava hacia abajo en el intervalo [xi , xi+2 ].
O sea, tenemos ademas informacion acerca de la concavidad de la
funcion.
4. En general considerando los n + 1 puntos consecutivos
(xi , fi ) , (xi+1 , fi+1 ) , ..., (xi+n , fi+n )
podemos definir recursivamente la n-esima diferencia dividida o diferencia dividida de orden n por :
f [xi , ..., xi+n ] =

f [xi+1 , xi+n ] f [xi , xi+n1 ]

.
xi+n xi

Vemos que hay una relacion entre las diferencias divididas y el valor de
los coeficientes del polinomio de Newton (2.16).
Los calculos para determinar las diferencias divididas pueden facilitarse
mucho, si estas se representan abreviadamente en una tabla.

 POLINOMICA.

2.1. INTERPOLACION
xi
x0

f [xi ]
f [x0 ]

f [xi , xi+1 ]

f [xi , xi+1 , xi+2 ]

9
f [xi , xi+1 , xi+2 , xi+3 ]

f [x0 , x1 ]
x1

f [x1 ]

f [x0 , x1 , x2 ]
f [x0 , x1 , x2 , x3 ]

f [x1 , x2 ]
x2

f [x2 ]

x3

f [x3 ]

x4

f [x4 ]

f [x1 , x2 , x3 ]
f [x1 , x2 , x3 , x4 ]

f [x2 , x3 ]
f [x2 , x3 , x4 ]

f [x2 , x3 , x4 , x5 ]

f [x3 , x4 ]
f [x3 , x4 , x5 ]
f [x4 , x5 ]
x5

f [x5 ]

x6

f [x6 ]

f [x5 , x6 ]

Tabla 2. Diferencias Divididas

Al observar la tabla 2 podemos darnos cuenta que las diferencias divididas
que aparecen en la diagonal que va desde el extremo superior izquierdo al
extremo inferior derecho son precisamente los coeficientes de (2.16), usando
este hecho el polinomio interpolante puede expresarse en forma simplificada
como:
Pn (x) = f [x0 ] + f [x0 , x1 ](x x0 ) + f [x0 , x1 , x2 ](x x0 )(x x1 ) +
+f [x0 , x1 , ..., xn ](x x0 )(x x1 ) (x xn1 ).

Actividad III.
Consideremos nuestro ejemplo modelo y construyamos el polinomio interpolante en la forma de Newton.
C
omo medimos qu
e tan buena es la aproximaci
on ?
Cuando en calculo diferencial se estudia el polinomio de Taylor de grado
n para aproximar una funcion f C (n+1) [a, b] en un punto x0 [a, b] su desarrollo involucra el termino Rn (x) llamado residuo o error de truncamiento
y que se define como :
Rn (x) =

f (n+1 )()
(x x0 )(n+1)
(n + 1)!

, (x0 , x).

(2.17)

Cabe pensar entonces que al aproximar una funcion mediante un polinomio

de interpolacion usando mas de un punto, tambien existe un termino de error
involucrado y de hecho es una extension de (2.17). En el siguiente resultado
formalizamos este hecho.
Teorema 2.1 Si x0 , x1 , ..., xn son (n + 1) puntos distintos en [a, b] y f
C (n+1) [a, b] entonces, para cada x [a, b] existe un (x) en (a, b) tal que :

10

E(x) = f (x) Pn (x) =

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

(2.18)

donde Pn (x) es el polinomio interpolante de Lagrange.

Demostracion. Se deja como inquietud, sin embargo podemos comentar
que se usa como parte del argumento el teorema de Rolle.

n 2.2
Observacio
Es posible obtener una relaci
on entre diferencia dividida y el termino
f (n+1 )()
,
(n + 1)!
de modo que podemos expresar la f
ormula de error en terminos de
diferencias divididas. Para ello consideremos :
Pn (x) = f [x0 ] + f [x0 , x1 ](x x0 ) + f [x0 , x1 , x2 ](x x0 )(x x1 ) +
+f [x0 , x1 , ..., xn ](x x0 )(x x1 ) (x xn1 )
si agregamos un punto m
as (xn+1 , f (xn+1 )) de modo que ahora tenemos
(n + 2) puntos y usamos la notaci
on :
xn+1 = x

f (xn+1 ) = f (
x)

x
6= xi

i = 0, ..., n

entonces de la f
ormula de Newton
f (
x) = Pn+1 (
x) = Pn (
x) + f [x0 , x1 , ..., xn , x
]w(
x)
o
f (
x) Pn (
x) = f [x0 , x1 , ..., xn , x
]w(
x)
de aqu y (2.18), obtenemos
f [x0 , x1 , ..., xn , x
] =

f (n+1) ()
para alg
un I
(n + 1)!

luego reemplazando f [x0 , x1 , ..., xn , x

] en (2.18) tenemos,
E(
x) = f (
x) Pn (
x) = f [x0 , x1 , ..., xn , x
]w(
x).

(2.19)

que es otra forma de expresar el error de interpolaci

on.
El polinomio de Newton puede expresarse para datos equiespaciados,
usando el concepto de diferencias progresivas y regresivas, para mayores detalles consultar la bibliografa recomendada al final.

 POLINOMICA.

2.1. INTERPOLACION

2.1.4

11

El problema de Hermite

El problema mas simple de Hermite es encontrar el polinomio H(x) del

menor grado posible, tal que interpole a la funcion f (x) y a su derivada
f 0 (x) en (n + 1) puntos x0 , , xn , es decir, el polinomio debe satisfacer
H(xi )

= f (xi ) , i = 0, ..., n

H 0 (xi ) = f 0 (xi ) , i = 0, ..., n

y el Problema General consiste en suponer que :
1. En el nodo i se conocen los valores
p(j) (xi ) , p(j1) (xi ) , p0 (xi ) , p(xi ).

2. Sea ki el n
umero de condiciones de interpolacion sobre el nodo xi .
3. Si x0 , ..., xn son los nodos entonces
p(j) (xi ) = f (j) (xi ) , j = 0, 1, ki1 , i = 0, ..., n

De los dos problemas presentados estudiaremos el caso simple, aquellos

polinomios que geometricamente se parecen a la funcion en los puntos de
interpolacion (xi , f (xi )) en el sentido que las rectas tangentes al polinomio
y a la funcion tienen el mismo valor, es decir poseen misma derivada en una
vecindad del punto de interpolacion.

2.1.5

C
omo obtener el polinomio de Hermite ?

Podemos determinar el polinomio de Hermite basandonos en la forma de

Lagrange o extendiendo la forma de Newton.
1. Polinomio de Hermite basado en la forma de Lagrange.
Consideremos los nodos x0 , ..., xn y supongamos que en cada nodo tenemos:
p(xi ) = f (xi ) , i = 0, ..., n
p0 (xi ) = f 0 (xi ) , i = 0, ..., n
Observemos que tenemos (n+1) condiciones sobre f (x) y (n+1) condiciones
sobre f 0 (x) en total (2n + 2) condiciones, por lo que el polinomio es de
grado (2n + 1). Procediendo en forma similar a la usada para determinar la
forma del polinomio interpolante de Lagrange, se propone que la forma del
polinomio de Hermite es :

H2n+1 (x) =

Pn

i=0 f (xi )hi (x)

Pn

i=0 f

0 (x )h

i i (x).

(2.20)

12
i (x) elementos de una base para el espacio vectorial de los
Con hi (x) y h
polinomios, tales que,
hi (xj ) = ij

i (xj ) = 0.
, h

h0i (xj )

, h0 i (xj ) = ij .

(2.21)
=0

Al aplicar las condiciones de interpolacion y las relaciones (2.21) el lector

puede comprobar que se obtiene que :
hi (x) = [1 2(x xi ) li0 (xi )]li2 (x) , i = 0, ..., n
(2.22)
i (x) = (x xi ) l2 (x) , i = 0, ..., n.
h
i
i (x) son de grado (2n + 1) y como antes li (x)
Notese aqu, que hi (x) y h
es el polinomio de Lagrange de grado n definido anteriormente, por:
li (x) =

n
Y

(x xj )
.
(xi xj )
j=0 j6=i

Resumiendo estas ideas el El polinomio de Hermite de grado (2n + 1) que

interpola a una funcion f (x) y a su derivada en (n + 1) puntos es :

H2n+1 (x) =

Pn

i=0 f (xi )hi (x)

Pn

i=0 f

0 (x )h

i i (x).

(2.23)

con,
hi (x) = [1 2(x xi ) li0 (xi )]li2 (x) , i = 0, ..., n.
(2.24)
i (x) = (x xi ) l2 (x) , i = 0, ..., n.
h
i
Que podemos decir del error cometido ?
Debera deducirse en forma analoga a como se obtiene el error del polinomio interpolante de Lagrange. En efecto,
Si x0 , ..., xn son (n + 1) puntos distintos en [a, b] y f C (2n+2) [a, b] y
H(x) es el polinomio de grado (2n + 1) tal que :
H(xi )

= f (xi ) , i = 0, ..., n

H 0 (xi ) = f 0 (xi ) , i = 0, ..., n

entonces, a cada x [a, b] le corresponde (x) [a, b] tal que,

E(x) = f (x) H(x) =

n
f (2n+2) ((x)) Y
(x xi )2 .
(2n + 2)! i=0

(2.25)

 POLINOMICA.

2.1. INTERPOLACION

13

n 2.3 Es conveniente se
Observacio
nalar que en el caso particular que existe un s
olo nodo x0 , se conoce la funci
on f (x) y sus n derivadas el polinomio
de Hermite es el polinomio de Taylor, es decir
H(x) = P (x)
= P (x0 ) + P 0 (x0 )(x x0 ) + P 00 (x0 )

(x x0 )2
(x x0 )n
+ + P (n) (x0 )
.
2!
n!
(2.26)

2. Polinomio de Hermite basado en la forma de Newton.

La idea es ahora extender el metodo de las diferencias divididas de Newton para resolver el problema de interpolacion de Hermite. Ello se logra
construyendo la tabla como antes, con la variante de incluir un nodo en la
tabla tantas veces como condiciones dadas existan sobre el. Para ilustrar
esta situacion consideraremos a modo de ejemplo el caso en que se tienen
dos puntos de interpolacion x0 y x1 y se conocen los valores de la funcion y
su derivada. Inspirandonos en el polinomio de Taylor podemos intuir que el
polinomio que interpola la funcion y su derivada en los nodos dados es de la
forma :

H3 (x) = a0 + a1 (x x0 ) + a2 (x x0 )2 + a3 (x x0 )2 (x x1 ).

(2.27)

donde obviamente los coeficientes a0 , a1 , a2 y a3 se pueden obtener de la

siguiente tabla :
xi
x0

f [xi ]
f [x0 ]

f [xi , xi+1 ]

f [xi , xi+1 , xi+2 ]

f [xi , xi+1 , xi+2 , xi+3 ]

f [x0 , x0 ] = f 0 (x0 )
x0

f [x0 ]

f [x0 , x0 , x1 ]
f [x0 , x1 ]

x1

f [x0 , x0 , x1 , x1 ]

f [x1 ]

f [x0 , x1 , x1 ]
f [x1 , x1 ] = f 0 (x1 )

x1

f [x1 ]

Tabla 4. Diferencias Divididas para el polinomio de Hermite

Veremos ahora porque,

f [x0 , x0 ] = f 0 (x0 )

y f [x1 , x1 ] = f 0 (x1 ).

En efecto,

f [x0 , x0 ] = lim f [x0 , x] = lim

xx0

xx0

f (x) f (x0 )
= f 0 (x0 ).
x x0

(2.28)

14
Un razonamiento similar es valido para:

f [x1 , x1 ] = lim f [x1 , x] = lim

xx1

xx1

f (x) f (x1 )
= f 0 (x1 ).
x x1

(2.29)

De donde el polinomio de Hermite es,

H3 (x) = f (x0 )+f 0 (x0 )(xx0 )+f [x0 , x0 , x1 ](xx0 )2 +f [x0 , x0 , x1 , x1 ](xx0 )2 (xx1 ).
(2.30)
n 2.4 En general se puede demostrar que, si los puntos xi son
Observacio
distintos y f (n) (x) es continua en el nodo xk , entonces :
lim

xk+n xk

f [xk , , xk+n ] =

f (n) (xk )
, xk = xk+1 = = xk+n .
n!

(2.31)

Actividad I.
1) Encuentre un polinomio p(x) de grado menor o igual que 4 que cumpla
con las condiciones :
f (0) = 0

f 0 (0) = 0

f (2) = 16

f 0 (2) = 32

f 00 (2) = 48

Usando
el primer nodo como pivote
el u
ltimo nodo como pivote.

2) PEP 2 (Problema 3). Segundo Semestre del 2000.

Un automovil recorrio 160 km en un camino que une dos ciudades, y
empleo, en este trayecto, 2 hora y 20 minutos. Use los datos de la
tabla adjunta para determinar :
a) La distancia recorrida en los primeros 45 minutos de viaje. (Use
un polinomio c
ubico.)
b) El tiempo empleado en llegar a la mitad del camino recorrido.
Tiempo (en min)
0
10
30
60
90
120
140

Distancia (en Km)

0.00
8.00
27.00
58.00
100.00
145.00
160.00

3 ) PEP 2 (Problema 1). Segundo Semestre del 2001.

Considere la funcion f (x) = ex + senx + 2 definida en [0, ].

 POLINOMICA.

2.1. INTERPOLACION

15

a) Aproxime la funcion dada, utilizando interpolacion lineal por

i
tramos y los puntos xi =
, i=0,1,2,...,4.
4
b) Mediante la aproximacion encontrada en el punto anterior, estime
el valor de la funcion en el punto x = 1.25 y determine su error
absoluto comparando el valor aproximado con el exacto.
4) PEP 2 (Problema 2). Segundo Semestre del 2001. Usando mnimos
cuadrados, se desea determinar una funcion que represente los datos
de la siguiente tabla:
k
xk
f (xk )

0
1
5.12

1
3
3.00

2
6
2.48

3
9
2.34

4
15
2.18

a) Encuentre una aproximacion de mnimos cuadrados usando la

funcion

h(x) = +
x
b) Utilizando los datos de la tabla y la aproximacion anterior, calcule
:
max{|h(xk ) f (xk )|/k = 0, 1, ..., 4}
Que interpretacion puede darle a este n
umero ?
5) Determine un polinomio P (x) tal que
p(1) = 10

p0 (1) = 11

p00 (1) = 22

p(3) = 18

p0 (3) = 83.

M
odulo 11
Prof. Mara Angelica Vega U.

9.2.

Polinomios de Tchebyshev y economizaci

on de series.

El estudio de lo polinomios de Tchebyshev conduce a los resultados siguientes:

i) Permiten una colocaci
on optima de los puntos de interpolacion para minimizar el error en la
interpolaci
on de lagrange.
ii) Permiten reducir el grado del interpolante con una perdida mnima de precision.
Definici
on 9.1 Se definen los polinomios de Tchebyshev por
Tn (x) = cos(n arc cos (x))

1x1

n = 0, 1, 2,

(142)

la sustituci
on = arc cos (x) implica que:
1x1

Tn (x) = cos(n )

(143)

Usando identidades trigonometricas la ecuacion 143 se transforma en

Tn+1 (x) = 2xTn (x) Tn1) (x) para cada

n = 1, 2,

y x [1, 1]

con T0 (x) = 1 y T1 (x) = x. Esta relacion de recurrencia implica que para cada n 1,
polinomio de grado n con coeficiente principal 2n1 .
9.2.1.

(144)
Tn es un

Propiedades o Teoremas importantes

Propiedad 9.1 El polinomio de Tchebyshev Tn d egrado n 1 tiene n ceros simplesn en [1, 1]



(2k 1)
x
ek = cos
(2n)

para

k = 1, 2, , n.

(145)

Adem
as tiene extremos absolutos en


(k)
xk = cos
(n)

para

k = 0, 1, , n.

(146)

con


(k)
Tn (xk ) = cos
(n)

para

k = 0, 1, , n.

(147)

Las gr
aficas de algunos polinomios de Tchebyshev se presentan en la siguiente figura: En estos
gr
aficos se puede observar las siguientes 3 propiedades de los polinomios de Chebyshev:
1. |Tn (x)| 1 1 x 1.
81

Figura 11: Graficas de T1 , ..., T5

2. Tn tiene n ceros reales distintos en el interior de [1, 1].

3. |Tn (x)| alcanza su m
aximo modulo de 1 en [1, 1] en n+1 puntos, incluyendo ambos extremos
1 y Tn (x) toma los valores 1 alternadamente en estos puntos.
Observaci
on 9.1
i) Los polinomios m
onicos (coeficiente inicial 1) Ten de Tchebyshev, se
pueden deducir de los polinomio de Tchebyshev Tn , dividiendo por 2n1 ,
Te0 (x) = 1

Ten (x) = 21n Tn (x)

para cada

n 1.

ii) Estos polinomios satisfacen la relaci

on de recurrencia,
1
Te2 (x) = xTe1 (x) Te0 (x)
2

1
Ten+1 (x) = xTen (x) Ten1 (x)
4

para cada n 2.
iii) Debido a la relaci
on lineal entre Ten y Tn la propiedad 8.1, implica que los ceros de Ten ocurren
tambien en,


(2k 1)
para k = 1, 2, , n.
(148)
x
ek = cos
(2n)
y que los valores extremos de los Ten ocurren en


k
xk = cos
(n)

para

82

k = 0, 1, , n.

(149)

iv) En estos valores extremos para n 1

(1)k
Ten (xk ) = n1
2

para cada k = 0, 1, 2, , n

De esto u
ltimo, se deduce la siguiente propiedad
Propiedad 9.2 Los polinomios de la forma Ten , cuando n 1, tienen la propiedad
1
2n1

= m
ax |Ten (x)| max |Pn (x)|
x[1,1]

para toda

x[1,1]

Pn

Y
f
n

En esta u
ltima expresi
on, se produce la igualdad si y solo s Pn = Ten
Observaci
on 9.2 1. La respuesta a la interrogante d
onde colocar los nodos interpolantes para
minimizar el error en la interpolaci
on de Lagrange? Recordemos que: Si f C n+1 [1, 1] y
x 0, , xn son n
umeros distintos en el intervalo [1, 1], entonces para cada x [1, 1]
existe (x) en (1, 1) tal que,
f (x) = P (x) +

f (n+1) ((x))
(x x0 )(x x1 ) (x xn ),
(n + 1)!

(150)

donde P es el poinomio interpolante de Lagrange. Vimos que una alternativa de minimizar

el error es ubicar los nodos de una forma
optima, lo equivale a encontrar los nodos que
minimizan la expresi
on.
|(x x0 )(x x1 ) (x xn )|
sobre todo el intervalo [1, 1]. Como
(x x0 )(x x1 ) (x xn ).

(151)

es un polinomio m
onico de grado (n + 1), la propiedad anterior implica que este mnimo se
obtiene ss
(x x0 )(x x1 ) (x xn ) = Ten+1 (x).
(152)
Luego si se escoge xk como el (k + 1) cero de Ten+1 para cada k = 0, 1, n, es decir cuando
xk se escoja como


(2k + 1)
x
ek+1 = cos
para k = 1, 2, , n.
(153)
(2(n + 1))
el valor m
aximo de |(x x0 )(x x1 ) (x xn )| ser
a minimizado. Como,
m
axx[1,1] |Ten+1 (x)| = 2n

1
= max |(x x
f0 )(x x
f1 ) (x x
en+1 )|
2n
x[1,1]

maxx[1,1] |(x x0 )(x x1 ) (x xn )|.

para cualquier elecci
on de x0 , x1 , , xn en el intervalo [1, 1]
83

(154)

2.

Esta tecnica de elecci

on de puntos que minimizan el error interpolante puede extenderse
f
acilmente a un intervalo cerrado general [a, b] usando el cambio de variables
x
e=

1
[(b a)x + a + b]
2

umeros correspondientes x
ek
para transformar los n
umeros xk , en el intervalo [1, 1] en los n
en el intervalo [a, b].
Ejemplo 9.1 Sea f (x) = xex en [0, 1,5]. Construir dos polinomios interpolantes de grado a lo m
as
tres.
i) Usando los nodos igualmente espaciados x0 = 0, x1 = 0,5, x2 = 1, x3 = 1,5
ii) Usando los nodos de Tchebyshev.
iii) Hacer una tabla comparativa entre los valores exactos de la funci
on, el polinomio usando
e
nodos igualmente espaciados y el polinomio de Tchebyshev T4 y entre los errores absolutos
de los dos polinomios.
Sol.
i) Usando los nodos igualmente espaciados se obtiene:
l0 (x) = 1,3333x3 + 4,0000x2 3,6667x + 1,
l1 (x) = 4,000x3 10,000x2 + 6,0000x,
l2 (x) = 4,000x3 + 8,0000x2 3,0000x,
l3 (x) = 1,3333x3 2,0000x2 + 0,66667x.
De aqu el polinomio interpolante es
P3 (x) = 1,3875x3 + 0,057570x2 + 1,2730x.
ii) Usando los nodos de Tchebyshev. Transformamos los ceros x
ek = cos( (2k+1)
, k = 0, 1, 2, 3
8
de Te4 de [1, 1] a [0, 1,5], usando la transformaci
on lineal
x
ek = 0,75 + 0,75e
xk
obteniendose,
x
e0 = 1,44291 , x
e1 = 1,03701 , x
e2 = 0,46299 , x
e3 = 0,05709
Los polinomios de Lagrange para estos nodos son:
e
l0 (x) = 1,8142x3 2,8249x2 + 1,0264x 0,049728,
e
l1 (x) = 4,3799x3 + 8,5977x2 3,4026x + 0,16705,
e
l2 (x) = 4,3799x3 11,112x2 + 7,1738x 0,37415,
e
l3 (x) = 1,8142x3 + 5,3390x2 4,7976x + 1,2568.
y el polinomio interpolante de grado a lo m
as 3 es:
Pe3 (x) = 1,3811x3 + 0,044445x2 + 1,3030x 0,014357.

84

7.

Integraci
on Num
erica.
M
odulo 7
Prof. Mara Angelica Vega U.

7.1.

F
ormulas de Newton-Cotes.

Problema. Se trata de encontrar el valor de

b

f (x)dx
a

especialmente cuando no puede ser resuelta por metodos analticos.

Se obtiene aproximando f (x) en [a, b] mediante un polinomio de interpolaci
on, es decir,
f (x)
Rb
a

= Pn (x) + E(x) ,

f (x)dx =

Rb
a

Pn (x)dx +

Rb
a

entonces
E(x)dx

de modo que si
Pn (x) =

n
X

f (xi )li (x)

i=0

es el polinomio interpolante de Lagrange,

Z

f (x)dx =
a

n
X

Z
f (xi )

i=0

pi (x)dx.

(106)

En general, si el integrando es:

Z

w(x)f (x)dx =
a

n
X

Hi (x)f (xi ) + Eint .

(107)

i=0

donde w(x) es una funci

on peso, los Hi (x) se determinan imponiendo que Ec , el error de cuadratura
sea cero, si f es un polinomio de grado menor o igual que n.
7.1.1.

M
etodo del Trapecio.

Si en la f
ormula, (??), n=1 entonces los puntos de interpolacion son x0 y x1
Z

x1

f (x)dx =
x0

1
X
i=0

Hi (x)f (xi ) = H0 f (x0 ) + H1 f (x1 ).

luego hay que calcular, H0 y H1 , en efecto

Z

x1

H0 =

x1

l0 (x)dx =
x0

x0

x x1
1
dx =
x0 x1
h

x1

(x x1 )dx =
x0

h
.
2

Analogamente resulta:
H1 =

h
.
2

Luego,
R x1

f (x)dx =

x0

h
[f (x0 ) + f (x1 )].
2

(108)

denominada Regla del trapecio.

Que representa gr
aficamente esta formula ?
Error de cuadratura
Que podemos decir del error cometido ? Para ello integramos el error de la aproximacion polinomial, en efecto
f 00 ()
E = f (x) P1 (x) =
(x a)(x b)
2
utilizando el teorema del valor medio para integrales, se tiene que
Z

E(x)dx =

ET =
a

1
(b a)3 f 00 ().
12

(109)

Notese aqu que ba = h es la longitud del intervalo y que la formula (??) proporciona un resultado
exacto para polinomios de grado a lo mas 1. Esta formula se conoce con elnombre de m
etodo del
trapecio simple.
Una forma de disminuir el error de cuadratura, es considerar mas puntos de interpolacion por
ejemplo (n + 1), en este caso el intervalo [a, b]lo dividimos en n subintervalos de longitud h y en
cada subintervalo aplicamos el metodo anterior , obteniendose,

Rb
a

n1
X
ba
f (x)dx =
[f (x0 ) + f (xn ) + 2
f (x0 + ih)].
2n
i=1

denominada m
etodo del trapecio generalizada y el termino de error
ET =

1 h3 00
f ()
12 n2

de donde un estimativo del error es :

|ET |

h3
max |f 00 (x)|
12n2 axb

(110)

7.1.2.

M
etodo de Simpson.

Si en la f
ormula, (??), n=2 entonces los puntos de interpolacion son x0 , x1 y x2 y
Z

x2

f (x)dx =
x0

2
X

Hi (x)f (xi ) = H0 f (x0 ) + H1 f (x1 ) + H2 f (x2 ).

i=0

luego hay que calcular, H0 , H1 y H2 , en efecto

R x1
x0

f (x)dx =

R x2
x0

(x x1 )(x x2 )
(x x0 )(x x2 )
f (x0 ) +
f (x1 ))dx
(x0 x1 )(x0 x2 )
(x1 x0 )(x1 x2 )

R x2 (x x0 )(x x1 )
f (x2 )dx +
x0 (x x )(x x )
2
0
2
1

x2

x0

(x x0 )(x x1 )(x x2 ) (3)

f ((x))dx.
6

evaluando las integrales (Consultar Analisis Numerico de Burden-Faires.) se obtiene :

Z

x1

f (x)dx =
x0

h5
h
[f (x0 ) + 4f (x1 ) + f (x2 )] f (4) .
3
90

denominada Regla de Simpson simple. Como la derivada involucrada en el error es de cuarto

orden, entonces la f
ormula ser
a exacta para polinomios de grado menor que 4.
Observemos que si queremos generalizar este metodo a n subintervalos, el metodo requiere de dos
ba
intervalos consecutivos, por lo que n debe ser par. Generalizamos definiendo: n = 2m , h =
2m
tenemos
a = x0 < x1 < ... < x2m = b , xi = x0 + ih , i = 0, ..., 2m
entonces
Rb
a

#
"
m
m1
X
X
h
f (x2i1 ) + f (x2m )
f (x2i ) + 4
+f (x0 ) + 2
f (x)dx =
3
i=1
i=1
(111)
m
h X (4)

f (i ).
90 i=1
5

Un concepto importante, es saber determinar el grado de exactitud o de precision de una formula

de cuadratura. El grado de exactitud o precisi
on de una formula de cuadratura es el entero
positivo n tal que E(Pk ) = 0, para todos los polinomios Pk de grado menor o igual que n, pero
E(Pn+1 ) 6= 0 para alg
un polinomio de grado (n + 1).
Actividades 7.1 Problema 1. (PEP No 3. Sem I - 2000) Una partcula con masa m = 20 kg se
mueve a traves de un fluido y est
a sometida a una resistencia viscosa dada por:
R(v) = v 3/2

donde v es su velocidad. La relaci

on entre el tiempo t, la velocidad v, y la resistencia R est
a dada
por:
Z v(t)
m
dv (en segundos)
t=
R(v)
v0
donde v0 es la velocidad inicial.
a) Suponiendo que v0 = 15 m/seg y usando cuatro intervalos de igual longitud, determinar el
tiempo T que se requiere para que la partcula llegue a v(T ) = 7,0 m/seg.
b) Estime el error cometido en (a)
Problema 2. Considere la siguiente f
ormula de integraci
on numerica:
2

f (x)dx =
2

p
p
2
[5f (2 0,6) + 8f (0) + 5f (2 0,6)]
9

(1)

a) Usando la f
ormula (1), estime el valor de la siguiente integral:
Z

2[0,2969 x 0,126x 0,3516x2 + 0,2843x3 0,1015x4 ]dx

b) Determine el polinomio de mayor grado para el que la f

ormula (1) es exacta.

7.2.

Cuadratura Gaussiana.

Recordamos que las f

ormulas de Newton-Cotes se obtuvieron al integrar los polinomios interpolantes de Lagrange. Como el termino de error en el polinomio interpolante involucra la derivada
de orden (n + 1) de la funci
on que se esta aproximando, la formula es exacta cuando se aproxima
cualquier polinomio de grado menor o igual a n. Luego las formulas de Newton-Cotes tienen grado
de precisi
on de a lo menos n.
Todas las f
ormulas de Newton-Cotes requieren que se conozcan los valores de la funcion que se
va a integrar, en puntos equiespaciados ( datos tabulados), sin embargo si la funcion esta dada
explcitamente, los puntos para evaluar la funcion pueden escogerse de modo que la aproximacion
tenga una mayor precisi
on.
La cuadratura Gaussiana consiste en escoger los puntos de evaluacion de la funcion de modo
que la aproximaci
on sea
optima.
Problema. Se trata de escoger los puntos x0 , x1 , ..., xn en [a, b] y las constantes H0 , H1 , ..., Hn de
modo que el error obtenido en:
Z b
n
X
Hi f (xi )
f (x)dx
a

i=0

sea mnimo, para cualquier f . La mejor eleccion es la que maximice el grado de precision de la
formula.
Como los valores de los Hi son arbitrarios y los xi solo a que la funcion cuya integral se esta
aproximando, debe estar definida en estos puntos , hay (2n + 2) parametros involucrados. Si los
coeficientes de un polinomio se consideran tambien como parametros, entonces la clase mas grande
de polinomios para la cu
al es exacta la formula sera 2n + 1.
Actividades 7.2 Se considera la f
ormula de cuadratura:
Z

f (x)dx = A[f (x1 ) + f (x2 ) + f (x3 )) + E

1

i) Determinar las ecuaciones que permiten calcular A , x1 , x2 , x3 de modo que la f

ormula
sea exacta para polinomios del mayor grado posible.
ii) Si en la f
ormula anterior se impone x2 = 0. Determinar A , x2 , x3 . Cu
al es el mayor
grado del polinomio para el que la f
ormula obtenida es exacta?
Sol.

2
2
2
, x3 =
, A= .
2
2
3

x2 = 0 , x1 =

La f
ormula es exacta para polinomios de grado 3.

De la actividad anterior , es claro que para aplicar cuadratura Gaussiana, necesitaremos algunos
conceptos relacionados con conjuntos ortogonales de funciones.
Definici
on 7.1 Diremos que un conjunto de funciones {0 , 1 , , n } es ortogonal en [a, b],
con respecto a la funci
on peso w(x) 0, si
Rb
a

Rb
a

0 kj w(x)dx = 0 , cuando j 6= k y
0 kj w(x)dx

> 0 , cuando j = k

Algunas Propiedades importantes.

1) Si {0 , 1 , , n } es un conjunto ortogonal de polinomios definidos en [a, b] y i es de
grado i para cada i = 0, 1, , n, entonces para cualquier polinomio Q(x) de grado a lo mas
n, existen constantes u
nicas 0 , 1 , , n }, tales que:
Q(x) =

n
X

i i (x)

i=0

2) Si {0 , 1 , , n } es un conjunto ortogonal de polinomios definidos en [a, b] con respecto a

una funci
on de peso continua w y k es un polinomio de grado k para cada k = 0, 1, , n,
entonces k tiene k races en el intervalo (a, b).
Resumen Polinomios Ortogonales.
Polinomios de Legendre.
i) Son ortogonales con respecto a la funcion peso w(x) = 1 en [a, b] = [1, 1]
ii) F
ormula Recursiva.

Pn+1 =

2n + 1
n+1


xPn (x)

n
n+1


Pn1 (x)

para todo n 1, P0 (x) = 1, P1 (x) = x.

iii) F
ormula para los pesos.
Wj =

2
, j = 1, ..., n
(n + 1)Pn+1 (xj )Pn0 (xj )

En esta f
ormula n representa el grado del polinomio.
iv) F
ormula para el error.
E(f ) =

22n+1 (n!)4
f (2n) ().
(2n + 1)[(2n)!]3

con [1, 1] y n grado del polinomio.

Polinomios de Laguerre.
i) Son ortogonales con respecto a la funcion peso w(x) = ex en [a, b] = [0, )

ii) F
ormula Recursiva.

Ln+1 =

1
n+1


(2n + 1 x)Ln(x)

n
n+1


Ln1 (x)

para todo n 1, L0 (x) = 1, L1 (x) = 1 x.

iii) F
ormula para los pesos.
Wj =

[(n)!]2
.
Ln+1 (xj )L0n (xj )

iv) F
ormula para el error.
E(f ) =

(n!)2 f (2n) ()
.
(2n)!

con [0, ).
Polinomios de Hermite.
2

i) Son ortogonales con respecto a la funcion peso w(x) = ex en [a, b] = (, )

ii) F
ormula Recursiva.
Hn+1 = 2xHn (x) 2nHn1 (x)
para todo n 1, H0 (x) = 1, H1 (x) = 2x.
iii) F
ormula para los pesos.
Wj =
iv) F
ormula para el error.

2(n+1) (n)!
.
[Hn0 (xj )]2

n! f (2n) ()
E(f ) =
.
2n [(2n)!]

con (, ).
Polinomios de Tchebyshev.
i) Son ortogonales con respecto a la funcion peso w(x) = (1 x2 )1/2 en [a, b] = [1, 1]
ii) F
ormula Recursiva.
Tn+1 (x) = 2xTn (x) Tn1 (x)
para todo n 1, T0 (x) = 1, T1 (x) = x.
iii) F
ormula para los pesos.

.
n
donde n representa el n
umero de ceros del polinomio.
Wj =

iv) F
ormula para el error.
E(f ) =

2
f (2n) ().
22n (2n)!

con (1, 1).

Actividades 7.3

1.

Estimar
3

Z
1

dx
x

usando un polinomio de Legendre de grado 3. Estimar adem

as el error.
Sol.
r
3
3
x0 = ( ) , x1 = 0 , x2 = ( )
5
5
5
8
5
H0 = , H 1 = , H 2 =
9
9
9
1
, I 1,098039 , |E| < 0,045714
g(y) =
y+2
r

2.

Estimar

ex x2 dx ,

usando un polinomio de Hermite de grado 2.

Sol.

3.

2
2

sqrt
, x2 =
, H 1 = H2 =
, I=
2
2
2
2

x1 =

Calcular :
Z

cosx
dx
1 x2

uasndo un polinomio de Tchebyshev de grado 2.

Sol. I 2,38884
4.

La f
ormula de Gauss - Tchebyshev,
Z 1
n
f (x)
X

f (xi )
dx approx
n i=1
1 x2
1
es exacta para polinomios de grado 2n-1.
a) Demuestrelo para el caso n=3.
b) Aproxime usando lo anterior:
Z
0

cosx
dx.
2 x2

7.3.

Integrales Indefinidas

Supongamos que desamos evaluar la integral:

Z

F (x) =

f (t)dt.

(112)

para a x b. Podemos escoger un adecuado paso de h para aproximar el lado derecho de (112),
para x = a + h , x = a + 2h , ... cada integral puede estimarse por alguna regla de cuadratura por
ejemplo trapecio. Los valores intermedios de F (x) pueden estimarse por interpolacion.
Otro metodo es reemplazar f por alg
un polinomio de aproximacion p y aproximar F integrando p.
Si |f (x) p(x)| <  para a x b entonces,

 Z
Z x


F (x)
p(t)dt <

a

dx (b a)

para a x b. Tomando [a, b] como [1, 1] una adecuada eleccion de p es la serie de Tchebyshev
b
X

r Tr (x).

(113)

r=0

donde,
r =

N
2 X
f (xj )Tr (xj ).
N j=0

(114)

y xj = cos( j
on (113), es necesario integrar Tr (x), cuya
N ), con N > n. Para integrar la expresi
integral indefinida es un polinomio de grado r+1. Por lo tanto, queda una suma de T0 , T1 , ..., Tr+1 .
Para evaluar en general la integral de Tr (x) es mas facil usar t = cos(), de modo que:
R
R
Tr (x)dx = cosr sin d
1
2

1
=
2

[sin (r 1) sin (r + 1)]d

cos(r + 1) cos(r 1)

r+1
r1


+C

Si r 6= 1, escogemos C de modo que se cumpla la condicion Tr (1) = (1)r , de modo que:



Z x
1 Tr+1 (x) Tr1 (x)
(1)r+1
Tr (t)dt =

+ 2
.
(115)
2
r+1
r1
r 1
1
donde r > 1. Adem
as,
Rx
1

Rx
1

T0 (t)dt

1
T1 (x) + 1
4

T1 (t)dt

1
1
T1 (x)
4
4

Luego integrando (113), tenemos:

n
X

!
r Tr (t) dt =

r=0

n+1
X

j Tj (x).

(116)

1 j n 1.

(117)

j=0

De (115) se deduce que :

j =

1
(j+1 j1 ) ,
2j

y
0 =

n
X
1
1
(1)r+1 r
0 1 +
2
4
r2 1
r=2

Si se define n+1 = n+2 = 0, la expresion (117) se cumple para j = n y n + 1.

7.4.

Integrales Impropias

Record
ando c
omo se definen las integrales impropias de primera y segunda clase, tenemos los dos
problemas siguientes:
Rb
i) Estimar a f (x)dx, donde f tiene una singularidad en x = b, entonces definimos sobre el
intervalo [a, b ], con 0 < < b a. Suponemos que
b

Z
lm

f (x)dx
a

R1
1
existe. Por ejemplo, 0 (1 x) 2 dx
R
ii) Estimar a f (x)dx, donde f esta definida sobre cualquier intervalo [a, b] y
Z
lm

existe. Por ejemplo,

R
1

f (x)dx
a

x2 dx

Una tecnica que es u

til algunas veces en ambos casos es hacer un cambio de variables adecuado.
R1
1
1
As la sustituci
on t = (1x2 ) 2 transforma la integral impropia i) 0 (1x) 2 sin xdx en la integral
R1
propia 2 0 sin 1 t2 dt que puede estimarse usando alguna regla de cuadratura estandard .
Otra tecnica usada es remoci
on de la singularidad. Consideremos la integral impropia
Z
I=
0

xp cosxdx ,

0<p<1

con singularidad en x = 0. Podemos aproximar cosx por el primer termino de la serie de Taylor
alrededor de x = 0, nos da
1

I=

xp (cosx 1)dx.

dx +

La primera integral tiene el valor (1 p)1 . La segunda integral puede calcularse mediante un
metodo numerico, puesto que su integrando no tiene singularidad, ya que cosx 1 se comporta
como 21 x2 en x = 0. Sin embargo, xp (cosx 1) se comporta como 12 x2p , cuya segunda
derivada y las derivadas de orden superior son singulares en x = 0. Por lo tanto no es posible
garantizar que la integraci
on numerica de resultados exactos.
Ejemplo 7.1

7.5.

Integrales M
ultiples

Dscutiremos integrales sobre una region rectangular. Mediante un cambio de variables, el rectangulo
puede transformarse en cuadrado. Usando cualquier cuadratura en una dimension, tenemos :
R b R b
a


f (x, y)dx dy

R b PN

f
(x
,
y)
dy
i
i
i=0

PN

PN

i=0

i=0 i

Rb
a

f (xi , y)dy

P

N
j=0


j f (xi , yj ) ,

donde, en realidad xi = yi . Los n

umeros xi i representan los puntos y los pesos respectivamente,
en una dimensi
on. Luego para dos dimensiones tenemos:
Z

f (x, y)dx dy =
a

N
X

i j f (xi , yj ).

(118)

i,j=0

denominada unaregla de integraci

on del producto. Por ejemplo, para un producto con la regla
de Simpso compuesta en dos dimensiones con 4 subintervalos en cada direccion, el valor relativo a
los pesos 1, 4, 2, 4, 2, 1, son multiplicados, los resultados de muestran en la siguiente tabla:
Tabla relativa a los pesos de un producto usando Simpson
1
4
2
4
1

4
16
8
16
4

2
8
4
8
2

4
16
8
16
4

1
4
2
4
1

La suma de los n
umeros en la tabla es (1 + 4 + 2 + 4 + 1)2 Consideraremos la funcion f (x, y) 1,
por lo tanto la regla integrar
a la funcion sobre la region con area (b a)2 , cada n
umero en la tabla
2
(b a)
debe ser multiplicado por
para obtener los pesos i j de (118).
144
El estimativo de error de las reglas del producto, da un error estimado para la regla unidimensional.
Se puede obtener una cota de error basados en la regla de Simpson , con paso h. Si R denota la
regi
on cuadrada a x b , a y b, entonces:
Z

f (x, y)dx =
a

7.6.

N
X

i f (xi , y) (b a)

i,j=0

h4 (iv)
f (y , y).
180 x

(119)

Diferenciaci
on Num
erica

Se trata de encontrar metodos para aproximar f 0 dados los valores de x en ciertos puntos. Si p
es una aproximaci
on polinomial a f , podemos usar p0 como una aproximacion a f 0 . Sin embargo,
puede ocurrir que el m
odulo m
aximo de f 0 (x) p0 (x) sobre el intervalo [a, b] puede ser mucho mas
grande que el m
odulo m
aximo de f (x) p(x). Por lo tanto se buscaran otros metodos. La formula
del error es:
Y
f (n+1) (x )
.
(120)
f (x) pn (x) =
(x)
(n + 1)!
n+1
donde,

n+1 (x)

= (x x0 ) (x xn ). Derivando,
f 0 (x) p0n (x) = 0n+1 (x)

f (n+1) (x ) (n+1) (x) d (n+1)

+
f
(x ).
(n + 1)!
(n + 1)! dx

Recordamos la f
ormula de diferencias avanzadas:




s
s
pn (x) = pn (x0 + sh) = f0 +
4f0 + +
4n f0
1
n
De x = x0 + sh se tiene

(121)

(122)

dx
= h, por lo tanto,
ds
ds d
1 d
pn (x0 + sh) =
pn (x0 + sh).
dx ds
h ds

(123)





1
1
d
s
4f0 + (2s 1)42 f0 + +
4n f0 .
n
h
2
ds

(124)

p0n (x) =
As
p0n (x) =

Para calcular 0n+1 (xr ) escribimos:

n+1 (x) = (x xr )j6=r (x xj ).

(125)

De aqu obtenemos:
0n+1 (x)


=


d
d
(x xr ) j6=r (x xj ) + (x xr ) j6=r (x xj ).
dx
dx

(126)

En x = xr , el segundo termino sera cero y obtenemos:

0n+1 (xr ) =

(xr xj ) = (1)nr hn r!(n r)!.

(127)

j6=r

ya que xr xj = (r j)h. Por lo tanto, (121) puede expresarse como:

f 0 (xr ) p0n (xr ) = (1)nr hn

r!(n r)! (n+1)

f
(x )r.
(n + 1)!

(128)

De aqu podemos deducir algunas reglas:

Con n = 1 en (124) obtenemos:
f1 f0
.
h

(129)

1
f1 f0
1
4f0 =
hf 00 (0 ).
h
h
2

(130)

p01 (x) =
y con r = 0 en (128) tenemos:
f 0 (x0 ) =

Con n = 2 en (124) y r = 1 obtenemos la derivada centrada en x1 , interpolando los puntos x0 , x1

y x2 , si x = x1 entonces s = 1 y la formula (124) se transforma en :
p02 (x1 ) =

1
f2 f0
1
[4f0 + 42 f0 ] =
.
h
2
2h

(131)

Esta es una aproximaci

on a f 0 (x1 ), con error (128).
f 0 (x1 )

f2 f0
1
= h2 f 000 (1 ).
2h
6

(132)

Se pueden obtener las f

ormulas de diferencia progresiva , centrada y regresiva respectivamente:
f (x + h) f (x)
.
h

(133)

f 0 (x)

f (x + h) f (x h)
.
2h

(134)

f 0 (x)

f (x + h) f (x h)
.
2h

(135)

f 0 (x)
y

Las derivadas de orden superior de f son aproximadas por las derivadas de orden superior del
polinomio interpolante pn . Si procedemos como antes, tenemos:
p00n (x) =

1 d2
pn (x0 + sh)
h2 ds2

y de (124) tenemos:
p00n (x)


 



1
d2
d2
s
s
2
3
n
= 2 4 f0 + 2
4 f0 + + 2
4 f0 .
3
n
h
ds
ds

(136)

donde x = x0 + sh y n 2. Para n = 2 encontramos que:

p002 (x) =

1 2
4 f0
h2

Consideramos esto como una aproximacion a f 00 (x) en x = x1 . As

f 00 (x1 )

f2 2f1 + f0
.
h2

(137)

Generalizando esta estrategia para cualquier k n obtenemos:

p(k)
n (x) =

1
hk

dk
dsk

s
k

4k f0 + +

dk
dsk

s
n


4n f0 .

(138)

De aqu se obtiene,
f (k) (x)

4k f0
.
hk

(139)

Si consideramos la f
ormula (139) para aproximar f (k) (x) para un x centrado en el intervalo [x0 , xk ]
y si k = 2m, entonces tenemos:
42m f0
f (2m) (xm )
.
(140)
h2m
Si k = 2m1, el intervalo [x0 , x2m1 ] no tiene un punto central tabulado, en este caso se generaliza
el procedimiento adaptado para la estimacion de la primera derivada. Usamos (140) con n = k+1 =
2m, encontrando
1
1
f (2m1) (xm ) 2m1 (42m1 f0 + 42m f0 ).
(141)
h
2
Ya que,
42m f0 = 4(42m1 f0 ) = 42m1 f1 42m1 f0 .

(142)

obtenemos:
f (2m1) (xm )

42m1 f0 + 42m1 f1
.
2h2m1

(143)

De (143) y (140) se pueden obtener aproximaciones para la tercera y cuarta derivada:

1
(f4 2f3 + 2f1 f0 ).
h3

(144)

1
(f4 4f3 + 6f2 4f1 + f0 ).
h4

(145)

f (3) (x2 )
y
f (4) (x2 )

Se pueden obtener otras reglas de derivacion considerando mas terminos en (138)

y de (154) tenemos:
p00n (x)







1
d2
d2
s
s
2
3
n
= 2 4 f0 + 2
4 f0 + + 2
4 f0 .
3
n
h
ds
ds

(166)

donde x = x0 + sh y n 2. Para n = 2 encontramos que:

p002 (x) =

1 2
4 f0
h2

Consideramos esto como una aproximacion a f 00 (x) en x = x1 . As

f 00 (x1 )

f2 2f1 + f0
.
h2

Generalizando esta estrategia para cualquier k n obtenemos:







1 dk
dk
s
s
k
n
p(k)
(x)
=
4
f
+

+
4
f
0
0 .
n
k
n
hk dsk
dsk

(167)

(168)

De aqu se obtiene,
f (k) (x)

4k f0
.
hk

(169)

Si consideramos la f
ormula (169) para aproximar f (k) (x) para un x centrado en el intervalo [x0 , xk ]
y si k = 2m, entonces tenemos:
42m f0
.
(170)
f (2m) (xm )
h2m
Si k = 2m1, el intervalo [x0 , x2m1 ] no tiene un punto central tabulado, en este caso se generaliza
el procedimiento adaptado para la estimacion de la primera derivada. Usamos (170) con n = k+1 =
2m, encontrando
1
1
f (2m1) (xm ) 2m1 (42m1 f0 + 42m f0 ).
(171)
h
2
Ya que,
42m f0 = 4(42m1 f0 ) = 42m1 f1 42m1 f0 .

(172)

obtenemos:
f (2m1) (xm )

42m1 f0 + 42m1 f1
.
2h2m1

(173)

De (173) y (170) se pueden obtener aproximaciones para la tercera y cuarta derivada:

f (3) (x2 )

1
(f4 2f3 + 2f1 f0 ).
h3

(174)

y
1
(f4 4f3 + 6f2 4f1 + f0 ).
h4
Se pueden obtener otras reglas de derivacion considerando mas terminos en (168)
f (4) (x2 )

90

(175)

M
odulo de clases
Prof. Mara Angelica Vega U.

9.7.

Mnimos cuadrados Discretos.

Consiste en buscar una una funcion aproximada que ajuste la forma de la tendencia que siguen los
datos sin pasar necesariamente por los puntos dados. Un modo de hacer esto es determinar una
curva que minimice la diferencia entre los puntos dados y la curva. Este metodo se conoce como
regresi
on por mnimos cuadrados. El caso mas simple de una aproximacion por mnimos cuadrados
es el ajuste de las observaciones mediante una lnea recta. La figura siguiente muestra la mejor
lnea recta.
Regresi
on lineal. Este problema clasico puede enunciarse : Dada una funcion continua f sobre
un intervalo [a, b] y un conjunto de n observaciones (x1 , y1 ), (x2 , y2 ), , (xn , yn ), buscamos una
lnea recta de ecuaci
on
f (x) = a0 + a1 x
tal que la diferencia entre el valor real de y y el valor aproximado obtenido por f (x) sea el menor
posible. Recordemos que se definio el error o residuo entre la funcion y) y su aproximacion como :
e = y a0 a1 x.
C
omo obtener la recta de regresi
on ? La regresion lineal es un metodo poderoso en el sentido
que ajusta los datos a la mejor lnea. Sabemos que este es un problema de optimizacion y debemos
encontrar los valores a0 y a1 de f (x) = a0 + a1 x de modo que la suma
Er =

n
X

(ei )2 =

i=1

n
X

(yi f (xi ))2 .

(176)

i=1

sea mnima. Para ello derivamos la ecuacion (176) con respecto a cada coeficiente e igualamos a
cero, de donde tenemos :
Er
a0

= 2

Pn

a0 a1 xi ) = 0

Er
a1

= 2

Pn

a0 a1 xi )xi = 0

i=1 (yi

i=1 (yi

De aqu se tiene el sistema de ecuaciones lineales, escrito en la forma Aa = b:



  Pn
Pn
a0
Pni=1 x2i
Pni=1 yi
=
a1
i=1 xi yi
i=1 xi

 Pn
Pni=1 1
i=1 xi


Pnn

i=1

xi

o equivalentemente


  Pn

Pn
a0
yi
i=1
P
Pni=1 x2i
=
n
a1
i=1 xi
i=1 xi yi

Resolviendo el sistema matricial obtenemos los valores de a0 y a1 .

91

Ejemplo 9.2 Ajuste a una lnea recta los siguientes cinco puntos. Que puede decir del error que
se obtiene
P1 (1, 5,12) , P2 (3, 3) , P3 (6, 2,48) , P4 (9, 2,34) , P5 (15, 2,18).
Soluci
on. Determinamos cada uno de los terminos que aparecen en (9.7)
C
alculo de las sumas.
n=5
P5
xi = 34
Pi=1
5
y = 15,12
P5i=1 i2
x
i=1 i = 352

,
,
,
,

P5
xi = 34
Pi=1
5
yi = 15,12
Pi=1
5
x
i yi = 82,76
Pi=1
5
2
y
i=1 i = 51,5928

a0
a1

de donde tenemos el sistema de ecuaciones :

Sistema de ecuaciones


5
34

34
352




=

15,12
82,76

Usando eliminaci
on gaussiana encontramos los valores del vector a
a0 = 4,15298

a1 = 0,166027.

de donde la ecuaci
on de la recta es:
L(x) = 4,15 0,17 x
Actividades 9.4
i) Encuentre el sistema matricial que se obtiene al ajustar, usando:
f (x) = P (x) = a0 + a1 x + a2 x2 .
ii) Dado un conjunto de m observaciones (x1 , y1 ), ..., (xm , ym ). aproximar los datos con un polinomio de grado n, usando el procedimiento de mnimos cuadrados. Observe que el polinomio
es ahora
n
X
f (x) = P (x) =
ai xi .
i=1

Problema 2. Supongamos que queremos ajustar n datos, P1 (x1 , y1 ), , Pn (xn , yn ) con una funci
on g(x) de dos par
ametros 1 y 2 , es decir
g(x) = 1 1 (x) + 2 2 (x).
Notemos que 1 (x) y ( x) son funciones conocidas elementales de x que no son necesariamente
elementos de una base para los polinomios. Obviamente para encontrar la g(x) que minimiza el
cuadrado de los errores hay que determinar los parametros 1 y 2 . En efecto,
E(g) =

n
X

(yi g(xi ))2 .

i=1

92

(177)

Derivando parcialmente con respecto a 1 y 2 , igualando a cero y reordenando se obtiene el

siguiente sistema de ecuaciones en forma matricial A = b
 Pn

  Pn

Pn
1
1 (xi )yi
i=1
Pni=1 1 (xi )1 (xi ) Pni=1 2 (xi )1 (xi )
P
=
n
2
i=1 1 (xi )2 (xi )
i=1 2 (xi )2 (xi )
i=1 2 (xi )yi
Las ecuaciones del sistema se denominan ecuaciones normales y cualquier modelo lineal (suma
ponderada de dos funciones) de dos parametros pueden ajustarse a un sistema de ecuaciones
lineales. El error en este caso puede determinarse luego de encontrar 1 y 2 , usando:
E(g) =

n
X

(yi )

i=1

n
X

1 (xi )yi + 2

i=1

n
X

!
2 (xi )yi

(178)

i=1

Para ilustrar este modelo, tomemos el ejemplo 0.2 y ajustemos los datos usando una funcion de
dos par
ametros como sigue.
Ejemplo 9.3 Ajustar mediante el metodo de mnimos cuadrados, los datos del ejemplo 1.15 con
2
una hiperbola desplazada del tipo g(x) = 1 + . Determine el error cometido al ajustar los datos
x
con la recta y con la par
abola.

9.8.

Algo sobre Normas de funciones.

Definici
on 9.2 La norma de una funci
on f que pertenece a alguna clase de funciones continuas
C es una funci
on de C al conjunto de los n
umeros reales no negativos :
|| ||

: C R+
f 7 ||f || 0

tal que verifica las siguientes propiedades :

i) ||f || > 0

ii) R y f C
iii) f, g C

||f || = 0 f = 0
||f || = |||f ||

||f + g|| |||f || + ||g||

Algunas normas conocidas con f C

i)
||f || = max |f (x)|.
axb

(179)

Esta norma definida sobre [a, b] es llamada norma del m

aximo, norma infinito o norma
de Tchebyshev. Verificar como ejercicio que satisface las propiedades de una norma.

93

ii)
(Z

) p1

|f (x)|p dx

||f ||p =

(180)

conocida como p norma, donde f C[a, b] y p 1. Obviamente cada valor de p 1 es

una norma, aunque en la practica se utiliza p = 1 y p = 2.
Normas Discretas (es decir para f C(X))
i)
||f || = max |f (xi )|.
0iN

(181)

llamada norma infinito discreta.

ii)

p1
p
.
||f ||p = N
i=0 |f (xi )| dx

(182)

conocida como p norma discreta, siendo p 1.

La desigualdad triangular es difcil de demostrar para esta u
ltima norma, excepto para p = 1 y
p = 2.
Nota 9.1
i) Observemos que se ha usado ||f || para denotar dos normas diferentes (179) sobre
C[a, b] y (181) an
alogamente sucede con ||f ||p .
ii) Existe la siguiente relaci
on entre la norma (181) y las p normas (182) dada por :
lm

p1
p
N
= max |f (xi )|.
i=0 |f (xi )| dx
0ixN

(183)

iii) La mejor aproximaci

on respecto a la norma infinito se llama aproximaci
on minimax,
puesto que esta aproximaci
on minimiza el error m
aximo ||f q|| sobre todo q Pn (x)
iv) La mejor aproximaci
on respecto a la norma 2 se llama aproximaci
on mnimo cuadrado.

94

M
odulo N 9
Profesores : Edgard Shipley - Mara Angelica Vega U.
Erkl
arungen: Zusammengestohlen aus Verschiedenem diesem und jenem
Ludwig van Beethoven.

9.9.

Aproximaci
on de Mnimos Cuadrados Continuos.

Problema 2. Sea f C[a, b] , queremos determinar un polinomio Pn de grado n que minimice

el error :
Z b
(f (x) Pn (x))2 dx
a

Sean Pn (x) =

n
X

ak xk

Z
y E(a0 , ..., an ) =

(f (x) Pn (x))2 dx.

(184)

k=0

es decir, la funci
on error es,
b

(f (x)

E(a0 , ..., an ) =
a

n
X

ak xk )2 dx.

(185)

k=0

Para determinar los coeficientes ak , k = 0, 1, 2, ..., n, debemos resolver el siguiente sistema:

E
= 0 para j = 0, 1, 2, ..., n.
aj
De (??), obtenemos
Z

b
2

f (x) dx 2

E(a0 , ..., an ) =
a

n
X

ak

f (x)x dx +
a

k=0

n
X

!2
k

ak x

dx.

(186)

k=0

Luego,
E
aj

= 2

Rb

= 2

Rb

f (x)xj dx +

Rb
a

f (x)xj dx + 2

Pn

k=0

Pn

k=0

ak

ak xk xj dx

Rb
a

xj+k dx = 0 , j = 0, 1, ..., n.

De aqu el sistema es:

n
X
k=0

Z
ak
a

xj+k dx =

f (x)xj dx , j = 0, 1, ..., n.

(187)

Las ecuaciones (??) son las ecuaciones normales y tienen solucion u

nica si f C[a, b]. La matriz
de los coeficientes de este sistema es:
 j+k+1

b
aj+k+1
H = (hi,j ) =
j+k+1
95

que es una matriz de Hilbert mal condicionada que produce dificultades con el error de redondeo.
Observemos que al dar valores a j en (??), se obtiene:
j = 0,

R

j = 1,

R

..
.

..
.

j = n, 2

R

b
a

f (x)dx +

b
a

f (x)xdx +

b
a

f (x)xn dx +


k
a
x
dx = 0
k
k=0

R b Pn
a

R b Pn

k=0

R b Pn
a




ak xk xdx = 0

k=0




ak xk xn dx = 0

Ejemplo 9.4 (Burden) Determinar el polinomio de grado dos que aproxima por mnimos cuadrados a la funci
on f (x) = sin (x) en el intervalo [0, 1].
Sol. El sistema (??) queda:

1
2

1
3

1
2

1
3

1
4

1
3

1
4

1
5

a0

R1
0

sin xdx

R1
a1 =
x sin xdx

R1 2
a2
x sin xdx
0

2 4
3

El polinomio aproximante de grado dos para f (x) = sin (x) en el intervalo [0, 1] es
P2 (x) = 0,050465 + 4,12251x 4,12251x2
9.9.1.

M
etodo Generalizado

Podemos usar en lugar de la base canonica para polinomios,{1, x, x2 , ..., xn } , cualquier base de
familia de polinomios ortogonales en el intervalo [a, b] respecto de una funcion peso w(x). Para ello
debemos tener presente los siguientes conceptos:
Definici
on 9.3 Un espacio con producto interno (E.P.I.) es un (E, h , i) donde:
hf, hi = hh, f i
hf, ah + bgi = ahf, hi + bhf, gi
hf, f i = ||f ||2 > 0

96

Propiedad 9.1 Algunas propiedades importantes en un E.P.I.

Pn
Pn
h 1 ai fi , gi
= 1 ai hfi , gi
||f + g||2

= ||f ||2 + 2hf, gi + ||g||2

Si f g

||f + g||2 = ||f ||2 + ||g||2

|hf, gi|

= ||f ||||g||

||f + g||2 + ||f g||2

= 2(||f ||2 + ||g||2 )

Teorema 9.1 Sea G un subespacio de un espacio con producto interno E. Para toda f E , g G
son equivalentes las proposiciones:
i) g es una mejor aproximaci
on a f G
ii) f g G

si

hf g, hi = 0

Definici
on 9.4 Si B = {0 , 1 , 2 , ..., n } es una base de los polinomios ortogonales en [a, b]
respecto de w(x) se cumple:
(
Z b
0,
si j 6= k
hj , k i =
w(x)j (x)k (x)dx =
k , si j = k
a
Observaci
on 9.1 Recordemos que,
Un conjunto de funciones 0 , 1 , ..., n es linealmente independiente en [a, b] si
n
X

cj j = 0 cj = 0 para cada j = 0, 1, ..., n

j=0

De lo contrario decimos que el conjunto es linealmente dependiente.

Sea j un polinomio de grado j para cada j = 0, 1, ..., n y sea 0 , 1 , ..., n un conjunto linealmente independiente. Si n es el conjunto de todos los polinomios de grado n, entonces
cualquier polinomio de n puede expresarse como
Pn (x) = nj=0 cj j (x)
.
Pn
De modo que el problema de determinar un polinomio Pn (x) =
j=0 cj j (x) de grado n que
minimice el error

2
Z b
Z b
n
X
E(c0 , ..., cn ) =
w(x)[f (x) Pn (x)]2 dx. =
w(x) f (x)
cj j dx.
(188)
a

97

j=0

se reduce a determinar la solucion del sistema de ecuaciones lineales obtenido mediante

h
i
Rb
Pn
E
= 0 a w(x) f (x) j=0 cj j (x) i (x)dx = 0
ci
E
=0
ci

Pn

j=0 cj hj , i i

= hf, i i ci i = hf, i i , i = 0, 1, ..., n

Por tanto los coeficientes de Pn (x) estan dados por:

hf, j i
cj =
=
j

Rb
a

(x)f (x)j (x)dx

j = 0, 1, ..., n..
Rb
(x)2j (x)dx
a

El producto interior definido por hf, gi =

que est
a dada por

Rb
a

p
||f ||2 = hf, f i =

(189)

(x)f (x)g(x)dx permite definir la llamada Norma-2

Z

! 21

b
2

(x)f (x) dx

.[a, b]

(190)

Con lo anterior,el problema de minimizar el error (??), puede expresarse como:

Problema 2. Determinar
mn ||f Pn ||2 , con Pn n
donde n es el espacio de los polinomios generados por B.
Ejemplo 9.5 Sean 0 (x) = 3 , 1 (x) = x + 2 , 2 (x) = x2 + 4x 6. Expresar
P (x) = a + bx + cx2 como nj=0 cj j (x) En efecto,
1
0 (x) = 1
3

2
x = 1 (x) 2 = 1 (x) 0 (x)
3

y
2
x2 = 2 (x) 4(1 (x) 0 (x)) + 20 (x)
3
Luego,

P (x) = a






1
2
14
0 (x) + b 1 (x) 0 (x) + c 2 (x) 41 (x) + 0 (x)
3
3
3

y
P (x) =

1
(a 2b + 14)0 (x) + (b 4c)1 (x) c2 (x)
3

Proceso de Gramm-Schmidt
Existe un metodo para generar una familia de polinomios ortogonales en un intervalo dado, respecto
de una funci
on peso (x). Este metodo es el llamado metodo de Gramm-Schmidt que consiste en
determinar los polinomios de la siguiente manera:
98

1. 0 (x) = 1

1 (x) = (x B1 )0 (x)

2. B1 se determina de modo que h1 , 0 i = 0. Se obtiene ;

Rb
(x)x0 (x)2 dx
B1 = aR b
(x)20 (x)dx
a
3.

Para k 2, k (x) = (x Bk )k1 Ck k2

Rb
a

Bk = R b
a

Rb

(x)x2k1 (x)dx

(x)2k1 (x)dx

Ck =

(x)xk1 (x)k2 dx
Rb
(x)2k2 (x)dx
a

Ejemplo 9.6 Los polinomios de Legendre son ortogonales en [1, 1] con funci
on peso (x) = 1 y
se define 0 (x) = 1. Aplicando la f
ormula anterior recursivamente, se obtiene
1 (x) = x , 2 (x) = x2

1
3
, 3 (x) = x2 x
3
5

.
Ejemplo 9.7 Encontrar la aproximaci
on de mnimos cuadrados de grado dos, de la funci
on f (x) =
ex con x [1, 1] usando polinomios de Legendre.
Sol.: Sea P2 (x) = a0 (x) + b1 (x) + c2 (x) el polinomio que minimiza ||f P2 ||22 . En este caso
las ecuaciones que determinan los coeficientes son:
R1

ex 0 (x)dx
1
R1 2
(x)dx
1 0

a =

R1
b

ex 1 (x)dx
1
R1 2
(x)dx
1 1

R1
= 1
R1

ex xdx

R1

R1

ex 2 (x)dx
1
R1 2
(x)dx
1 2

1
1
[ex ]11 = (e e1 )
2
2

x2 dx

3(2)
3
= = 3e1
2e
e

ex (x2 31 )dx
15
dx =
(e 7e1 )
R1
1
2
4
(x

)
3
1

El polinomio cuadr
atico que minimiza la aproximaci
on por mnimos cuadrados es:
P2 (x)

= a0 (x) + b1 (x) + c2 (x)

= 21 (e e1 )0 (x) 3e1 1 (x) +

15
4 (e

7e1 )2 (x)

En el gr
afico, se muestra la funci
on y su aproximaci
on : El error es 0,001440569656, luego ||ex
1
P2 (x)||2 = 0,3795483706e

99

Figura 11: f (x) = ex

Ejemplo 9.8 Encontrar las aproximaciones por mnimos cuadrados a

y(t) = t2 en (0, 1), usando polinomios de Legendre de grado 1.
Sol.
1) El polinomio de Legendre de grado 0 es P0 (x) = 1 y P1 (x) = x
2) Transformaci
on del intervalo (0, 1) a (1, 1).
Consideremos la recta x = at + b, si t = 0 x = 1 b = 1
luego x = at 1 y si t = 1 x = 1 a = 2
luego, la recta de transformaci
on es x = 2t 1
y=

(x + 1)2
4

t=

adem
as

x+1
2 ,

es decir (0, 1) (1, 1) y

P0 (x) = 1

P1 (x) = x

Usando (4)
ak

a0

a1

2k + 1
2

1
=
2

3
2

1
1

y(x)Pk (x)dx
1

(x + 1)2
1
dx =
4
3
(x + 1)2
1
xdx = .
4
2

Luego,
y(x) =

1
1
1 x
P0 (x) + P1 (x) = +
3
2
3 2
100

y
y(t) = t

1
6

Ejemplo 9.9 Dada f (x) = e2x en [1, 3]. Determinar la mejor aproximaci
on a
g(x) = a + bx2 + cx4 .
Sol. La base {1, x2 , x4 } genera el subespacio G E, consideremos g G , f E y recordemos
que
hf g, hi = 0 hf, hi = hg, hi.
Si h = 1, x2 , x4 entonces,
h=1

, he2x , 1i

= ha + bx2 + cx4 , 1i

R3

= ha, 1i + hbx2 dx, 1i + hcx4 , 1i

e2 dx

0,0664
, he2x , x2 i
R3
1

R3
1

bx2 dx +

R3
1

x2 dx + b

R3
1

cx4 dx

x4 dx + c

, he2x , x4 i

= ha + bx2 + cx4 , x4 i

e2 x4 dx = a

0,4967

R3
1

cx6 dx

26
a + 48,4b + 312,92c
3

e2 x2 dx = ha, x2 i + hbx2 dx, x2 i + hcx4 , x2 i

0,153

R3

R3

= ha + bx2 + cx4 , x2 i

=a

h = x4

adx +


3
x3 3
x5 
= 2a + b 1 + c  .
3
5 1
26
= 2a + b + 48,4c
3

e2x 3
2 1

h = x2

R3

R3
1

x4 dx + b

R3
1

x6 dx + c

R3
1

cx8 dx

= 48,4a + 312,24b + 2186,8c

Resolviendo el sistema
a = 0,144
2

b = 0,045

c = 0,0034

Luego g(x) = 0,144 0,045x + 0,0034x .

9.9.2.

Sistemas Ortonormales.

Existe otro enfoque de los problemas de aproximaci

on en un E.P.I, que consiste en construir
un sistema ortonormal.
101

Definici
on 9.5 Un conjunto {1 , , n } es ortonormal ssi hi , j i = ji
Teorema 9.2 Sea {1 , 2 , , n } un sistema ortonormal
en un E.P.I. que denotaremos (E, h, i).
Pn
La mejor aproximaci
on a f mediante un elemento i=1 ci i se obtiene ssi ci = hf, i i.
hf

n
X

ci i , j i = hf, j i

i=1

n
X

ci hi , j i = hf, i i ci = 0

i=1

Es decir, la mejor aproximaci

on de f (x) se puede escribir en la forma
f (x) =

n
X

hf (x), i (x)ii (x)

con

{i (x)}

base ortonormal.

i=1

El procedimiento sugerido de estos conceptos para aproximar elementos de un espacio E mediante

elementos de un subespacio G es :
i) Obtener una base ortonormal {1 , 2 , ..., n }. para G
ii) Determinar la mejor aproximacion a partir de
n
X

ci hi , j i =

i=1

n
X

ci i

i=1

Ejemplo 9.10 Obtener la aproximaci

on de la forma g(x) = c1 + c2 x2 + c3 x3 , a f (x) = sen (x),
usando la base ortonormal:

x 5x3 3x 63x5 70x3 + 15x

q , q
q
,
2

2
2
3

11

Sol.
i) Determinamos cada ci , i = 1, 2, 3
q R
q
1
c1 = 32 1 x sen (x)dx = 2 32 (sen (1) cos(1))
c2

c3

q R
7 1
2

11
2

(5x3 3x) sen (x)dx = 2

1

7
2 (18 sen (1)

28cos(1))

q
11
5
3
(63x

70x
+
15x)
sen
(x)dx
=
2
2 (4320 sen (1) 6728cos(1))
1

R1

ii) La mejor aproximaci

on es
r
r
r
3
7
11
2
g(x) = 2
(sen (1)cos(1))x+2
(18 sen (1)+28cos(1))x +2
(4320 sen (1)6728cos(1)).
2
2
2
Actividades 9.5 Determina la mejor aproximaci
on de la funci
on f (x) = sen (x) mediante un
polinomio g(x) = c1 x + c2 x3 + c3 x5 en [1, 1]. (Usando || || ).
102

9.9.3.
1.

Proceso de ortonormalizaci
on de Gramm-Schmidt.

Definimos
u1 =

2.

v1
,
||v1 ||

luego ||u1 || = 1 , w1 = u|

Se define w2 = v2 hv2 , u1 iu1

hw2 , u1 i = hv2 hv2 , u1 iu1 , u1 i
= hv2 , u1 i hv2 , u1 ihu1 , u1 i1 = 0 w2 u1 .
w2 u1 , pero w2 no tiene necesariamente norma 1. Por tanto ,
u2 =

3.

v2 hv2 , u1 iu1
w2
=
||w2 ||
v2 hv2 , u1 iu1

Se define w3 = v3 hv3 , u1 iu1 hv3 , u2 iu2 ortogonal a u1 y u2 .

hw3 , u1 i = hv3 hv3 , u1 iu1 hv3 , u2 iu2 , u1 i
= hv3 , u1 i hv3 , u1 ihu1 , u1 i hv3 , u2 ihu2 , u1 i
hw3 , u1 i = hv3 , u1 i hv3 , u1 i = 0.
Luego w3 u1 y w3 u2 se comprueba en forma analoga. Como w3 u1 y w3 u2 , definimos
u3 =

w3
||w3 ||

Por tanto {u1 , u2 , u3 } son ortonormales.

4.

Generalizamos a un =

wn
con
||wn ||

wn = vn hvn , u1 iu1 hvn , u2 iu2 hvn , un1 iun1 .

Ejemplo 9.11 Ilustrativo. Sea f (x) = e2x con x [1, 3] encontrar la mejor aproximaci
on a f (x)
i) respecto de la base = {x, 3x2 + 1} con w(x) = 1 en || ||2 .
ii) respecto de la base = {x + 1, x2 } con w(x) = ex en || || .
Sol. i) Nuestro prop
osito es encontrar g(x) =

P
hf (x), gi (x)igi (x). Para ello

Se debe ortonormalizar, para obtener g1 (x) y g2 (x) respecto a = {v1 (x) = x, v2 (x) =
3x2 + 1}

103

g1 (x) =

v1 (x)
,y
||v1 (x)||2
3

||v1 (x)||22

= hv1 (x), v1 (x)i =

1

de donde, ||v1 (x)||2 =

26
3

x2 dx =

x3 3
26
| =
3 1
3

y
r

x
g1 (x) = q

26
3

g2 (x) =

1 (v1 (x)) dx =

3
x
26

w2 (x)
, donde w2 (x) = v2 (x) hv2 (x), g1 (x)ig1 (x)
||w2 (x)||2
r

hv2 (x), g1 (x)i =

(3x + 1)
1

Luego, w2 (x) = (3x2 +1)21,74

3
26 x

3
xdx = 21,74
26

es ortonormal a g1 (x), perono necesariamente unitario.

Falta determinar
||w2 (x)||22

Z
=

Luego, ||w2 (x)||2 =

[(3x + 1) 21,74

(w2 (x)) dx =
1

3 2
x] dx = 16,98
26

16,98 y

g2 (x) =

q
3
(3x2 + 1) 21,74 26
x
4,12

Luego, = {g1 (x), g2 (x)} es ortonormal seg

un la norma || ||2 y la mejor aproximaci
on es
g(x) = hf (x), g1 (x)ig1 (x) + hf (x), g2 (x)ig2 (x)
y
hf (x), g1 (x)i =

R3
1

e2x

3
26 xdx

= 0,033

q
3
(3x + 1) 21,74 26
x

hf (x), g2 (x)i =

R3
1

2x

4,12

xdx = 0,045

entonces la mejor aproximaci

on es:
r
g(x) = 0,033

3
(3x2 + 1) 21,74 26
x
3

x 0,045
26
4,12

104

M
odulo N 10
Profesores : Mara Angelica Vega U.
Erkl
arungen: Zusammengestohlen aus Verschiedenem diesem und jenem
Ludwig van Beethoven.

10.
10.0.4.

Aproximaci
on Minimax.
Caso Discreto.

Problema Este problema no es directo como el metodo de mnimos cuadrados. Por ejemplo si
1
queremos encontrar el polinomio minimax P1 (x) = a0 + a1 x que aproxima a f (x) = x 2 en un
intervalo [ 14 , 1]
1

E = mn 1max |x 2 (a0 + a1 x)|

[a0 ,a1 ]

4 x1

La mejor lnea recta es aquella para la cual el error maximo ocurre en x =

punto interior x = . Entonces
1
a0 + a41 ( 14 ) 2 = E
1

a0 + a1 2
a0 + a1 1

1
4

, x = 1 y en alg
un

= E
= E.

Observamos que hay 3 ecuaciones y 4 incognitas, a0 , a1 , E y . La cuarta ecuacion se obtiene del

hecho que la derivada es cero en x = para el error E = f (x) (a0 + a1 x). Es decir,
1 1
2 b=0
2
Ya que E en x =
a0 +

a1
4

( 14 ) 2

a0 + a1 1

=E
= E.

De aqu,
9
17
1
2
, =
, a=
, E=
3
16
48
48
17 2
1
Luego la mejor aproximaci
on es y(x) =
+ x, que aproxima con un error no mejor que
.
48 3
48
b=

10.1.

Algoritmo de Remez.

Permite construir aproximaciones minimax, debido al matematico ucraniano E.Ya. Remez.

Algoritmo de Remez Calcula iterativamente el polinomio minimax p Pn , para f C[a, b] en
un conjunto de (n + 2) puntos. La sucesion de polinomios p converge al polinomio minimax p Pn
y la sucesi
on de conjuntos X converge a X sobre el cual f p es equioscilante.
105

Paso 1. Escoger un conjunto X = {x1 , x2 , ..., xn+2 }.

Paso 2. Resolver el sistema de ecuaciones lineales,
f (xi ) p(xi = (1)i e

1in+2

(191)

para determinar un n
umero real e y un polinomio p Pn .
Paso 3. Cambiar el conjunto X, de la siguiente manera:
Cuando [a, b] es [1, 1], una forma com
un de elegir el conjunto X en paso 1 es el conjunto de puntos extremos de Tn+1 . Para este conjunto inicial, el paso 2 da el polinomio
minimax inmediatamente cuando f (x) = xn+1 . En el caso general las ecuaciones lineales
del paso 2 siempre tienen una u
nica solucion.
Para llevar a cabo el paso 3, necesitamos primero estimar ||f p|| , donde p Pn es
el polinomio que ha sido calculado en el paso 2.
Si ||f p|| es suficientemente cercana a |e| entonces por teorema 4.3 p es proximo al
polinomio minimax y terminamos el proceso iterativo.
De otra forma, sea un punto tal que
|f () p()| = ||f p|| .
Incluimos el punto en el conjunto X y borramos uno de los puntos existentes de modo
que f p tenga signos alternados sobre los (n + 2) puntos de X.
Si est
a entre a y x1 , descartamos x1 o xn+2 y anaalogamente si esta entre xn+2 y b,
descartamos x1 o xn+2 . De otra forma, esta entre dos xi consecutivos y en este caso
descartamos uno de ellos apropiadamente.
1

Ejemplo 10.1 Encontrar el polinomio minimax de grado 2, para f (x) = ( 83 x + 58 ) 2 en [1, 1].
Sol. Escogemos X = {1, 0,5, 0,5, 1} inicialmente, los puntos extremos de T3 y resolvemos el
sistema de ecuaciones lineales (??) de donde obtenemos:
p(x) = 0,79188 + 0,24665x 0,04188x2

|e| = 0,00335

(192)

con 5 cifras decimales.

Evaluando f p en intervalos intervalos de 0,01, encontramos que hay m
aximo y mnimo local en
1 , 0,61 , 0,40 y 1 y f p tiene signos alternados en estos puntos, los que tomamos como el
nuevo conjunto X. Resolviendo el sistema(??) nos d
a otra vez
p(x) = 0,79196 + 0,24650x 0,04196x2

|e| = 0,00350

(193)

Encontramos que ||f p|| = 0,00350 con 5 cifras decimales y por lo tanto aceptamos este p como
polinomio minimax.

106

10.1.1.

Caso Continuo.

En los espacios con Producto Interno (E.P.I.)se pueden definir las siguientes normas:
1.

Norma uniforme o norma minimax o norma de Tchebyshev.

||f || = ||f || = max |f (x)|
axb

.
2.

Norma Euclideana o norma 2.

s
Z

w(x)f 2 (x)dx.

||f ||2 =
a

3.

Norma 1.
b

Z
||f ||1 =

w(x)|f (x)|dx.
a

4.

Norma minimax ponderada. ||f || = maxaxb w(x)|f (x)|

Recordemos que, p Pn es una mejor aproximaci

on de f C con respecto a una norma
dada si
||f g|| = inf ||f g||
qPn

Las aproximaciones con respecto a la norma infinito se denominan aproximaciones

minimax, porque es la mejor aproximaci
on que minimiza el error m
aximo ||f
q|| q Pn

Problema 10.1
i) Determinar || x||2 .
ii) Hallar un polinomio g(x) = a0 + a1 x ortogonal a f (x) = 1.

iii) Encontrar el polinomio minimax para f (x) = 1 + x2 en [0, 1].

Sol.
i)

|| x||2 = h x, xi =

x ln |x|dx =
1

e2
e2
1

+ .
2
4
4

ii)
hg(x), f (x)i = 0

g(x)f (x)

ha0 + a1 x, 1)i = 0 ha0 , 1i + a1 hx, 1i = 0

Re
1

a0 1 ln |x|dx + a1

a0

Re
1

ln |x|dx + a1


a0 + a1
107

Re
1

Re
1

x 1 ln |x|dx = 0
x ln |x|dx = 0

e2
e2
1

+
2
4
4


=0

Tenemos infinitas soluciones. En t

erminos de a0 ,

a0 = a1

luego g(x) = a1

e2
e2
1

+
2
4
4

e2
e2
1

4
2
4


a0 = a1

e2
e2
1

4
2
4


,


+ a1 x es ortogonal a 1.

iii) Polinomio minimax. De (??).

a1

f (1) f (0)
= 21
1

a0

c

1 + f (c)
( 2 1)
2
2

donce c es tal que, f 0 (c) = a1

f 0 (x) =
luego,

x
c
f 0 (c) =
2
1+x
1 + c2

c
c
= 21
= 0,42 c = 0,455 f (c) = 1,0986
1 + c2
1 + c2

luego,

a0 =

1 + 1,0986
2

( 2 1)

1,0986
2


= 0,955

Por tanto, la aproximaci

on minimax lineal a f (x) en [0, 1] es:

P1 (x) = 0,955 + ( 2 1)x.

Observaci
on 10.1 Observamos que puede haber m
as de una mejor aproximaci
on a f .
Determinar la mejor aproximaci
on polinomial a f L [a, b] equivale a resolver el
problema : minimizar m
ax |f (x) pn (x)| si x [a, b] y pn n , es decir, se debe
determinar el polinomio de grado n que minimiza el valor m
aximo de |f (x) pn (x)|.
Proposici
on 10.1 Existe una u
nica mejor aproximaci
on que resuelve el problema mnimax.
Teorema 10.1 (Teorema de la alternancia para polinomios) Para cualquier funci
on
f (x) C[a, b] existe una u
nica mejor aproximaci
on polinomial mnimax pn (x) que
est
a caracterizada por la propiedad alternante de que hay (a lo menos) n+2 puntos en
[a, b] en los cuales f (x) pn (x) obtiene sus valores m
aximos absolutos (nominalmente
||f pn || ) con signos alternados.

108

Este teorema se atribuye a Tchebyshev. Pero Borel (1905) estableci

o que la condici
on
de alternancia es una condici
on necesaria y suficiente y que existe s
olo un polinomio
que cumple con dicha caracterizaci
on
Ejemplo 10.2 La funci
on f (x) = x2 se aproxima por f = p1 (x) = x 0,125 en [0, 1].
Entonces, el error f (x) p1 (x) = x2 x + 0,125 tiene un valor m
aximo absoluto en
los puntos 0, 0,5 y 1 y alcanza los valores 0,125, 0,125, 0,125 que tienen signos alternados. Por lo tanto, es la u
nica aproximaci
on mnimax. La gr
afica de la funci
on y la
aproximaci
on es:

Figura 12: f (x) = x2

109

10.2.

Propiedad mnimax de los polinomios de Tchebyshev

Los polinomios de Tchebyshev,

Tn (x) con x = cos tiene (n + 1) valores extremos en los


puntos x = yk = cos k
con k = 0, 1, 2, ..., n
n
Proposici
on 10.2 Los polinomios de Tchebyshev alcanzan su valor
m
aximo en mag

nitud igual a 1 en [1, 1] en los (n + 1) puntos x = yk = cos k
con
k
= 0, 1, 2, ..., n y
n
presenta signos alternados en dichos puntos.
La proposici
on evoca el teorema de la alternancia y se puede aplicar de la siguiente
forma:
Sea f (x) = xn y sea pn1 (x) su aproximaci
on polinomial mnimax de grado n 1 en
el intervalo [1, 1]. Por el teorema de la alternancia f (x) pn1 (x) = xn pn1 (x),
debe tener dicha propiedad en n + 1 puntos.
El coeficiente principal de Tn (x) es 2n1 de modo que 21n Tn (x) es de la misma
forma que f (x)pn1 (x) = xn pn1 (x) y con la misma propiedad de la alternancia.
Se sigue de lo anterior que
xn pn1 (x) = 21n Tn (x)
Claramente 21n Tn (x) es un polinomio m
onico de grado n. Del Teorema de la alternancia y de lo anterior es el siguiente corolario.
Corolario 10.1 La aproximaci
on polinomial mnimax de grado n 1 a la funci
on
f (x) = xn

en

[1, 1]

es

pn1 (x) = 21n Tn (x)

Corolario 10.2 (Propiedad mnimax de Tn (x) ) La aproximaci

on polinomial mnimax
en [1, 1] mediante un polinomio m
onico de grado n a la funci
on cero es 21n Tn (x)
Ejemplo 10.3 . La aproximaci
on mnimax a la funci
on cero en [1, 1], mediante un
polinomio m
onico de grado 4 es,
23 T4 (x) = 23 (8x4 8x2 + 1) = x4 x2 + 0,125
Este polinomio tiene la propiedad de la alternancia y tiene valores extremos
0,125 , 0,125 , 0,125 , 0,125 y 0,125 respectivamente, en los cinco puntos dados por
yk = cos

, k = 0, 1, 2, 3, 4
4

De aqu,
1
1
yk = 1 , p , 0 , p , 1
(2)
(2)
Tambi
en por el corolario 1, la aproximaci
on mnimax por un polinomio c
ubico a la
funci
on f (x) = x4 es p3 (x) = x4 (x4 x2 + 0,125) = x2 0,125.
110

El error f (x)p3 (x) tiene la propiedad alternante en los puntos yk , obtenidos anteriormente. Se debe destacar que x2 0,125 es una aproximaci
on polinomial cuadr
atica a la
funci
on x4 en [1, 1] y tambi
en presenta cinco extremos. El teorema de la alternancia
se cumple en n + 3 puntos.
El teorema se
nala que la alternancia de los signos se cumple para a lo menos n + 2
puntos.
Si el intervalo de aproximaci
on es [0, 1] en vez de [1, 1], se debe redefinir los polinomios
de Tchebyshev a dicho intervalo. Estos polinomios trasladados se denotan por Tn (x) y
su recorrido es [0, 1]. Para obtener Tn (x) a partir de Tn (x) se debe efectuar un cambio
de intervalo, usando la transformaci
on:
s = 2x 1

o bien

x=

1
(1 + s) , x [0, 1] , s [1, 1]
2

Se obtiene:
Tn (x) = Tn (s) = Tn (2x 1)
Luego:
T0 (x) = 1 , T1 (x) = x , T2 (x) = 8x2 8x + 1 , T3 (x) = 32x3 48x2 + 18x 1
y as sucesivamente. En general el polinomio m
onico de aproximaci
on mnimax de
grado n a cero es
212n Tn (x)
Por ejemplo el polinomio de aproximaci
on mnimax de grado dos a la funci
on cero en
[0, 1] es:
23 T2 (x) = 23 (8x2 8x + 1) = x2 x + 0,125
Ejemplo 10.4 Ilustrativo.(Aproximaci
on Minimax en una base ortonormal.)
2x
Sea f (x) = e
con x [1, 3] encontrar la mejor aproximaci
on a f (x) respecto de la
base = {x + 1, x2 } con w(x) = ex en || || .
Sol.
Se debe ortonormalizar, para obtener g1 (x) y g2 (x) respecto a .
g1 (x) =

v1 (x)
x+1
=
, y
||v1 (x)||
||v1 (x)||
||v1 (x)|| = max |v1 (x)| = 4 g1 (x) =
x[1,3]

g2 (x) =

x+1
4

w2 (x)
, donde w2 (x) = v2 (x) hv2 (x), g1 (x)ig1 (x)g1 (x)
||w2 (x)||
Z
hv2 (x), g1 (x)i =
1

ex x2

x+1
17 3 1
dx =
e + e = 86,04
4
4
4

111

Luego, w2 (x) = x2 86,04

x+1
.
4

Falta determinar
||w2 (x)|| = max |w2 (x)| = max |x2 86,04
x[1,3]

Sea p2 (x) = x2 86,04

x[1,3]

x 86,04

|
4
4

x 86,04

entonces
4
4
p02 (x) = 2x

86,04
= 0 x = 10,755
4

p002 (10,755) > 0 (10,755, f (10,755) es un mnimo. El m

aximo se alcanza en x = 1.
Luego
86,04 86,04
||w2 (x)|| = |1

| = 42
4
4
y
x+1
x2 86,04
4
g2 (x) =
42
Luego, la mejor aproximaci
on es
g(x) = hf (x), g1 (x)ig1 (x) + hf (x), g2 (x)ig2 (x)
x+1
dx = 0,21
4
2
R3
x 86,04 x+1
4
dx = 0,41
hf (x), g2 (x)i = 1 e2x ex
42
hf (x), g1 (x)i =

R3
1

e2x ex

entonces la mejor aproximaci

on es:
g(x) = 0,21

x+1
x 0,41
4

112

x2 86,04 x+1
4
42

10.3.

Potencias de xn en t
erminos de {Tn (x)}

Se puede demostrar que:

k=[ n
2]
n

1n

x =2

X n
Tn2k (x).
k

(194)

k=0

donde [m] representa la parte entera de m y se

nala que el k-
esimo t
ermino de la suma
debe dividirse por dos si n es par y k = n2 .
Ejemplo 10.5 Ilustrar para el caso n = 4. En efecto, tenemos
x4

= 23

 
4
T42k (x)
k=0
k

P2

 
 
4
4
T2 (x) + 21
T0 (x)
1
2
= 81 T4 (x) + 12 T2 (x) + 38 T0 = (x)

= 23 [T4 (x) +

Esto u
ltimo se usa en la Economizaci
on de Tchebyshev cuando una funci
on se aproxima por Taylor (Maclaurin) hasta un cierto grado y se rebaja dicho grado manteniendo
un cierto orden de exactitud.
Ejemplo 10.6 f (x) = e(x) en [1, 1], con un error tolerable de 0,05, reemplazando el
t
ermino del desarrollo que contiene a x4 por los polinomios de Tchebyshev de grado
menor o igual a 4.

113

M
etodo de Diferencias Finitas
Prof: Mara Ang
elica Vega U.
M
odulo 12

12.

Resoluci
on Num
erica del Problema de Valor
de Contorno de orden 2.

Hemos visto dos metodos para determinar la solucion particular yp de la

solucion general de una ecuacion diferencial no homogenea con coeficientes
constantes y = yc + yp . Ambos metodos tienen sus ventajas y desventajas,
ademas de no ser aplicables para cualquier expresion de g(x).
El metodo que analizaremos ahora denominado M
etodo de diferencias finitas, tiene la ventaja de ser aplicable a cualquier problema de contorno, desde no lineal a ecuaciones diferenciales parciales para las que no
siempre existe un metodo analtico para determinar su solucion, pero tiene la
desventaja de dar un valor aproximado de la solucion exacta en algunos puntos del dominio de definicion de la solucion que son elegidos por el usuario.
Problema. Se trata de resolver el siguiente Problema de Valor de Con-

196

torno (P.V.C)
y 00 (x) = P (x)y 0 (x) + Q(x)y(x) + r(x) ,

0x1
(211)

y(0) = 0.
y(1) = 0.
con P (x), Q(x), r(x), funciones continuas en un intervalo abierto I.

Observemos que en el problema (211), la ecuacion diferencial es no homogenea con coeficientes variables.

12.1.

M
etodo de diferencias finitas.

Este metodo consiste en

Dividir el intervalo [0, 1] en N + 1 subintervalos de igual longitud, de

modo que:
xi = i h

i = 0, 1, , N + 1

(212)

donde,
h=

1
N +1

En cada nodo xi , se aproximan las derivadas de la ecuacion:

y 00 (xi ) = P (xi )y 0 (xi ) + Q(xi )y(xi ) + r(xi ) ,

i = 0, 1, , N + 1.
(213)

197

Existen varias formas de aproximar las derivadas, pero es recomendable usar diferencias centradas , puesto que resulatan aproximaciones
mejores,
y 00 (xi ) =

y 0 (xi )

y(xi+1 ) 2y(xi ) + y(xi1 )

h2

y(xi+1 ) y(xi1 )
=
.
2h

(214)

Reemplazamos las expresiones de (214) en (213), simplificamos la notacion escribiendo yi en lugar de y(xi ) y reordenamos los terminos de
la ecuacion en la siguiente forma :



h
h
2
r(xi )h = 1 + p(xi ) yi1 +(2+h Q(xi ))yi 1 p(xi ) yi+1 , i = 1, , N.
2
2
2

(215)

Observaci
on 12.1 Observemos que (215) es un sistema de ecuaciones lineales, cuyas inc
ognitas son los valores de la soluci
on aproximada evaluada en los siguientes puntos (extremos de intervalos),
xi

i = 1, , N.

Ejemplo 12.1 Use el metodo de diferencias finitas para resolver el

(P.V.C.):
y 00 =

2 0
(y 2) ,
t

1
y( ) = 3 ,
2
198

y(1) = 3

Soluci
on. Notemos que en esta ecuaci
on P (t) =

2
4
, Q(t) = 0 , r(t) = .
t
t

La ecuaci
on (215) en forma discreta queda,




h
8h2
h
yi1 (4 + 2h2 0)yi + 2 + 2
yi+1 =
, i = 1, , N.
22
ti
ti
ti
(216)
Si tomamos h = 0,1, tenemos que t0 = 0,5 , t1 = 0,6 , t2 = 0,7 , t3 = 0,8 y
t4 = 0,9 , t5 = 1,0. Sustituimos en (216) y obtenemos el sistema matricial :

5
73
8
0
0 y1
2( ) 3
4

60
15
12

16
8

4
4
0

y2

7
70
7
; =
= 35

1
7
8
9

4
0
y3

10
4
4
80

10
16
8
269
4
2( ) 3
.
0
0
y4
9
90
9
45
7
3

(217)

De aqu la soluci
on aproximada es,
y1 =

16
57
262
43
, y2 =
, y3 =
, y4 =
15
7
20
90

Actividades 12.1

Escriba el sistema de ecuaciones (215) en la forma matricial Ay = b

Dado el P.V.C.
y 00 + y + x = 0 , 0 < x < 1
y(0) = y(1) = 0.
Determine la soluci
on exacta del problema.
199

con

 Use el metodo de diferencias finitas, con N=3, para obtener una

soluci
on aproximada del problema.
Resuelva el P.V.C.
y 00 = 3y 0 + 2y + 2x + 3 , 0 < x < 1

con

y(0) = 2.
y(1) = 1.
use h=0.25 y compare los resultados con la soluci
on exacta.

200

Metodos Numericos para resolver Ecuaciones Diferenciales

Parciales.
Prof. Mara Angelica Vega U.
Modulo 14

Observaci
on 15.2 Este tema ha sido desarrollado, usando el texto Metodos
Numericos para Ingenieros. Steven C. Chapra - Raymond P.Canale

16.

Ecuaciones elpticas.

Se usan para caracterizar problemas en estado estacionario con valores

en la frontera.

16.1.

La ecuaci
on de Laplace

La ecuacion de Laplace, se utiliza para modelar diversos problemas que

tienen que ver con el potencial de una variable desconocida. Se usara una
placa calentada para deducir y resolver esta EDP elptica. En la figura se
muestra la placa y un elemento sobre la cara de una placa rectangular delgada de espezor 4z. La placa esta totalmente aislada excepto en sus extremos, donde la temperatura puede ajustarse a un nivel preestablecido. El
aislamiento y el espesor de la placa permiten que la transferencia de calor
222

Figura 12: Placa delgada de espesor 4z

esta limitada solamente a las dimensiones de x y y. En estado estacionario el

flujo de calor hacia el elemento en una unidad de tiempo 4t debe ser igual
al flujo de salida , es decir,
q(x)4y4z4t + q(y)4x4z4t = q(x + 4x)4y4z4t + q(y + 4y)4x4z4t
(247)
cal
donde q(x) y q(y) son los flujos de calor en x y y, respectivamente [ cm
2 s ]

Dividiendo entre 4z y 4t y reagrupando terminos, se obtiene

[q(x) q(x + 4x)]4y + [q(y) q(y + 4y)]4x = 0
Multiplicando el primer termino por

4x
4x

y el segundo por

4y
4y

se obtiene

q(x) q(x + 4x)

q(y) q(y + 4y)
4x4y +
4y4x = 0
4x
4y
223

(248)

Dividiendo entre 4x4y y tomando lmite, se llega a

q
q

=0
x y

(249)

La ecuacion (249) es una ecuacion diferencial parcial, que es una expresion

de la conservacion de la energa en la placa. Sin embargo no puede resolverse,
a menos que se especifiquen los flujos de calor en los extremos de la placa.
Debido a que las condiciones de frontera se dan para la temperatura, la
ecuacion (249) debe reformularse en terminos de la temperatura. la relacion
entre flujo y temperatura esta dada por la ley de Fourier de conducci
on
del calor, que se representa como,
qi = kC

T
t

(250)

donde,
cal
qi es el flujo de calor en la direccion de la dimension i [ cm
2 s ]

k es el coeficiente de difusividad t
ermica
es la densidad del material

cm2
s

g
cm3

C es la capacidad calorica del material [ gcal

C ]
T es la temperatura del material C que se define por
T =

224

H
CV

donde,

H es el calor (cal)
V es el volumen en cm3

Algunas veces, el termino que esta multiplicando a la derivada parcial en la

ecuacion (250) se trata como un solo termino.
k 0 = kC

(251)

cal
donde, k 0 es el coeficiente de conductividad t
ermica [ scm
C ]. En ambos

casos k y k 0 son parametros que determinan, que tan bien el material conduce
calor.
La ley de Fourier se conoce tambien como ecuaci
on constitutiva, debido a que proporciona un mecanismo que define las interacciones internas
del sistema. Inspeccionando la ecuacion (250) indica que la ley de Fourier especifica que el flujo de calor perpendicular al eje i es proporcional a gradiente
o pendiente de la temperatura en la direccion i. El signo negativo asegura
que un flujo positivo en la direccion i resulta de una pendiente negativa de
alta a baja temperatura. En la figura 13 se ilustra la representacion grafica
de un gradiente de temperatura, recpordamos que el calor se transfiere hacia
abajo, desde una temperatura alta a una baja, el flujo en a) va de izquierda

225

a derecha en la direccion i y en este caso la pendiente es negativa, es decir

un gradiente negativo se relaciona con un flujo positivo. En la figura 14 se
representa el caso inverso, donde el gradiente positivo se relaciona con un
flujo de calor negativo de derecha a izquierda. Sustituyendo, (250) en (249)

Figura 13:

T
i

<0

Figura 14:

T
i

>0

se obtiene,
2T
2T
+
=0
x2
y 2

(252)

conocida como ecuaci

on de Laplace. En el caso donde hay fuentes o perdidas de calor, dentro del dominio bidimensional, tenemos la ecuacion
2T
2T
+
= f (x, y)
x2
y 2

226

(253)

donde la funcion f (x, y) describe las fuentes o perdidas de calor. La ecuacion

(253) se conoce como ecuaci
on de Poisson.

16.2.

Soluci
on Num
erica

La solucion numerica, se basa en el metodo de diferencias finitas. En el

caso bidimensional tratamos la placa como una malla de puntos discretos.
Y
6
(x0 , yn+1 )
..
.

(xm+1, yn+1 )

yj+1

(xi , yj+1 )

yj

(xi , yj )

yj1

(xi , yj1 )

x1

xi1

xi

(x0 , y0 )

xi+1 (xm+1 , y0 ) X

Luego aproximamos las derivadas parciales en cada punto de la malla

transformando la ecuacion diferencial en una ecuacion algebraica en diferencias.

227

16.2.1.

La Ecuaci
on de Laplace.

Las diferencias centrales basadas en la malla de la figura son:

Ti+1,j 2Ti,j + Ti1,j
2T
=
2
x
4x2

(254)

Ti,j+1 2Ti,j + Ti,j1

2T
=
2
y
4y 2

(255)

las cuales tienen errores de O(4x2 ) y O(4y 2 ) respectivamente. Sustituyendo

en la ecuacion (252) se obtiene,
Ti+1,j 2Ti,j + Ti1,j
Ti,j+1 2Ti,j + Ti,j1
+
=0
4x2
4y 2

(256)

En la malla cuadrada de la figura anterior, 4x = 4y, reagrupando terminos

la ecuacion queda,
Ti+1,j + Ti1,j + Ti,j+1 + Ti,j1 4Ti,j = 0

(257)

Esta relacion, que se satisface por todos los puntos interiores de la placa, se
conoce como ecuaci
on laplaciana en diferencias. Debemos ademas especificar las condiciones de frontera en los extremos de la placa para obtener
una solucion u
nica. El caso mas simple es aquel donde la temperatura en la
frontera es un valor fijo. Este tipo de condicion se denomina condici
on de
frontera de Dirichlet. En la figura siguiente , se tiene una placa con este
tipo de condicion. En esta figura un balance en el nodo (1, 1) es de acuerdo
228

Figura 15: Condicion Dirichlet

con la ecuacion (257), se obtiene

T21 + T01 + T12 + T10 4T11 = 0

(258)

Sin embargo, T01 = 75 y T10 = 0, por lo tanto (258) se expresa como

4T11 + T12 + T21 = 75

Analogamente, se pueden desarrollar ecuaciones para los otros interiores.

229

El resultado es el sistema siguiente de nueve ecuaciones y nueve incognitas:+

4T11

T21

T12

T11 +4T21

T31

T21

+4T31

T11
T21

T22
T32
+4T12

T22

T12

+4T22

T32

T22

+4T32

T31

T13

T12

T23
T33
4T13

T22

T23

T13 +4T23
T32

T23

75

50

75

50

= 175
T33

= 100

+4T33 = 150
(259)

16.2.2.

El m
etodo de Liebmann.

En la mayora de las soluciones numericas de la ecuacion de Laplace se

tienen sistemas que son muchos mas grandes, por ejemplo para una malla de 10 por 10 nodos se tiene un sistema de 100 ecuaciones lineales. El
metodo com
unmente usado es el de Gauss-Seidel, que cuando se aplica a las
ecuaciones diferenciales parciales, tambien se conoce como el m
etodo de
Liebmann.

230

Con este metodo la ecuacion la ecuacion (258) se expresa como

Tij =

Ti+1,j + Ti1,j + Ti,j+1 + Ti,j1

4

(260)

que se resuelve de manera iterativa para j = 1...n y i = 1...m. Como la matris

asociada al sistema (257) es diagonal dominante, la sucesion de aproximaciones convergeraa a una solucion estable.
Como en el metodo coveccional de Gauus- Seidel las iteraciones se repiten
hasta alcanzar una presicion estipulada.
Ejemplo 16.1 Usando metodo de Liebmann, calcular la temperatura de la
placa calentada de la figura 15. Use sobrerrelajaci
on con = 1,5 e itere
hasta que e = 1 %.
Sol. La ecuaci
on (260) para i = 1 , j = 1 es
Tij =

0 + 75 + 0 + 0
= 18,75,
4

(261)

aplicando sobrerrelajaci
on se obtiene,
T11 = 1,5(18,75) + (1 1,5)0 = 28,125

(262)

Para i = 2 , j = 1
T21 =

0 + 28,125 + 0 + 0
= 7,03125
4
(263)

T21 = 1,5(7,03125) + (1 1,5)0 = 10,54688

231

Para i = 3 , j = 1
T31 =

50 + 10,54688 + 0 + 0
= 15,13672
4
(264)

T31 = 1,5(15,13672) + (1 1,5)0 = 22,70508

calculando sucesivamente se obtiene
T12 = 38,67188 T22 = 18,45703 T32 = 34,18579
(265)
T13 = 80,12696 T23 = 74,46900 T33 = 96,99554
Como los Tij son cero, inicialmente entonces los errores en la primera iteraci
on ser
a 100 % En la segunda iteraci
on los resultados son:
T11 = 32,51953 T21 = 22,35718 T31 = 28,60108
T12 = 57,95288 T22 = 61,63333 T32 = 71,86833

(266)

T13 = 75,21973 T23 = 87,95872 T33 = 67,68736

El error T1,1 se estima como sigue
|e11 | = |

32,51953 28,12500
|100 % = 13,5 %
32,51953

(267)

La novena iteraci
on da como resultado:
T11 = 4300061

T21 = 33,29755 T31 = 33,88506

T12 = 63,21152 T22 = 56,11238 T32 = 52,33999

(268)

T13 = 78,58718 T23 = 76,06402 T33 = 69,71050

donde el error m
aximo es 0,71 % En la figura siguiente se muestran los
resultados.
232

Figura 16: Distribucion de temperatura

16.3.

Condiciones en la frontera.

La condicion de frontera fija o de Dirichlet es una alternativa com

un en
las ecuaciones diferenciales parciales, otro tipo de condicion es la condici
on
de frontera de Newmann la que se tiene como dato la derivada en la frontera.
En el caso de la placa calentada, esto corresponde a especificar el flujo de
calor, mas que la temperatura en la frontera.
Un ejemplo es el caso donde el extremo esta aislado, en tal caso que
denominamos condici
on de frontera natural la derivada es cero. Esta conclusion se obtiene de la ecuacion (250), ya que aislar una frontera significa
233

que el flujo de calor debe ser cero.

El esquema siguiente muestra un nodo (0, j) en el extremo izquierdo de
una placa calentada, aplicando la ecuacion (257) se obtiene

T1,j + T1,j + T0,j+1 + T0,j1 4T0,j = 0

(269)

Observese que para aesta ecuacion se necesita un punto (1, j) que esta fuera
de la placa. Aunque parezca un pun to ficticio sirve para incorporar la derivada de la condicion de frontera en el problema, lo que se logra representando
la primera derivada en la dimension x en (0, j) por la diferencia dividida
finita
T1,j T1,j
T
=
x
24x2

(270)

despejando
T1,j = T1,j 24x

T
x

(271)

Ahora tenemos una relacion para T1,j que incluye la derivada, sustituyendo
en la relacion (269) se obtiene

2T1,j 24x

T
+ T0,j+1 + T1,j1 4T0,j = 0
x

De esta forma se ha incorporado la derivada a la ecuacion.

234

(272)

Ecuaciones Hiperbolicas.
Prof. Mara Angelica Vega U.
Modulo 16

18.

La ecuaci
on de ondas.

La ecuacion de ondas es un ejemplo de una ecuacion en derivadas parciales

2
2u
2 u
=
c
para 0 < x < a
t2
x2

0<t<b

(293)

o en forma mas simple

utt (x, t) = c2 uxx (x, t)

para 0 < x < a

0<t<b

(294)

con las condiciones de contorno,

u(0, t) = 0

u(a, t) = 0

0 t b,

u(x, 0) = f (x)

0 x a,

ut (0, t) = g(x)

0 < x a,

(295)

La ecuacion de ondas modela el desplazamiento u desde su posicion de equilibrio de una cuerda elastica vibrante cuyos extremos, de coordenadas x = 0
y x = a estan fijos. Aunque es posible determinar la solucion de ondas
por medio de las series de Fourier, se usara este problema como modelo de
ecuaciones hiperbolicas.
248

18.1.

Construcci
on de la ecuaci
on en diferencias.

Hacemos una particion del rectangulo R = {(x, t)/0 x a , 0 t

b} en una malla que consta de n 1 por m 1 rectangulos de lados 4x = h
y 4t = k, como la figura,
t
6
tm
..
.
tj+1
tj

tj1

xi1

xi

xi+1

xn

comenzamos por la fila de abajo, donde t = t1 = 0 y sabemos que la

solucion es u(xi , t1 ) = f (xi ). Ahora usaremos una ecuacion en diferencias
para calcular, en las filas sucesivas, las aproximaciones a la solucion exacta,
que en los puntos de la malla u(xi , tj ). Denotaremos para cada j = 2, 3, m
uij u(xi , tj )
Usaremos para aproximar utt (x, t) y uxx (x, t) las formulas:
utt (x, t) =

u(x, t + k) 2u(x, t) + u(x, t k)

+ O(k 2 )
k2
249

(296)

y
uxx (x, t) =

u(x + h, t) 2u(x, t) + u(x h, t)

+ O(h2 )
h2

(297)

El espaciado entre los puntos de la malla es uniforme en todas las filas

xi+1 = xi + h

y xi1 = xi h

y tambien uniforme en todas las columnas

tj+1 = ti + k

y tj1 = tj k

Sustituyendo (295) y (296) y usando la notacion mencionada tenemos,

ui,j+1 2ui,j + ui,j1
ui+1,j 2ui,j + ui1,j
= c2
k2
h2

(298)

que es la ecuacion en diferencias que se usara como aproximacion a la

ecuacion diferencial (293) o (294). Llamemos r =

ck
, luego tenemos:
h

ui,j+1 2ui,j + ui,j1 = r2 (ui+1,j 2ui,j + ui1,j )

(299)

Reordenando los terminos, determinaremos las aproximaciones a la solucion

en los puntos de la fila (j + 1)-esima de la malla, supuesto que conocemos
las aproximaciones a la solucion en los puntos de los puntos de las dos filas
anteriores, la j-esima y la (j 1)-esima.
ui,j+1 = (2 2r2 )ui,j + r2 (ui+1,j + ui1,j ) ui,j1

250

(300)

para i = 2, 3, , n 1. En la figura siguiente se muestra la posicion en la

malla de los cuatro valores conocidos , que aparecen en el miembro derecho
de (298) que se usan para determinar la aproximacion ui,j+1

ui,j1

r2 ui1,j

r2 ui+1,j

(2 2r2 )ui,j

ui,j1

18.2.

Condiciones Iniciales

Si se desea usar (296) para calcular las aproximaciones en los puntos de

la tercera fila de la malla, se necesitan las aproximaciones en los puntos de
la primera y segunda fila.
Los valores de la primera fila vienen dados por la funcion f (x), pero
no as los valores de la segunda fila, por lo que se usa la funcion g(x) para
determinar las aproximaciones en dichos puntos. Supongamos que x = xi
en la frontera inferior de R, aplicando teorema de Taylor de orden 1, para

251

desarrollar u(x, t) alrededor de (xi , 0) para el valoru(xi , k) , se verifica

u(xi , k) = u(xi , 0) + ut (xi , 0)k + (k 2 )

(301)

Ahora se usa uxi ,0 = f (xi ) = fi y uxi ,0 = g(xi ) = gi en (301) para obtener

las aproximaciones numericas en los puntos de la segunda fila t2 = k

ui,2 = fi + kgi

para

i = 2, 3, , n 1

(302)

En general, u(xi , t2 ) 6= ui,2 por lo que el error introducido al usar (302)

se propagara a toda la malla cuando se use el esquema dado por (299).
Para contrarrestar esta situacion es aconsejable escoger un tama
no de k
muy peque
no. Com
unmente se da el caso de que la funcion f (x) dada en el
contorno es dos veces derivable en el intervalo, con lo cual uxx (x, 0) = f 00 (x).
Esta igualdad nos permite usar la formula de Taylor de orden 2 para obtener
una aproximacion mejorada a los valores de la segunda fila de la malla. Para
ello volvemos a la ecuacion de ondas y usamos la relacion entre las derivadas
parciales segundas, obteniendo:
utt (xi , 0) = c2 uxx (xi , 0) = c2 f 00 (xi ) = c2

fi+1 2fi + fi1

+ O(h2 ) (303)
h2

Recordamos que la formula de Taylor de orden 2 es

u(x, k) = u(x, 0) + ut (x, 0)k +

252

utt (x, 0)k 2

+ O(k 2 )
2

(304)

y aplicando (304), (303)y (302) en el punto x = xi se obtiene:

c2 k 2
(fi+1 2fi + fi1 ) + O(k 2 ) + O(k 3 )
2h2

u(xi , k) = fi + kgi +

Como r =

ck
h,

(305)

la ecuacion (305) se puede simplificar obteniendose la sigu-

iente formula que proporciona aproximaciones numericas mejoradas a los

elementos de la segunda fila.
r2
(fi+1 + fi1 ),
2

ui,2 = (1 r2 )fi + kgi +

para i = 2, 3, n 1

(306)

Ejemplo 18.1 Se ilustrar

a el metodo de diferencias finitas en la resoluci
on
de la ecuaci
on de ondas de una cuerda vibrante.

utt (x, t) = 4uxx (x, t).

para 0 < x < 1sn 1

y 0 < t < 0,5

(307)

con las condiciones de contorno:

u(0, t) = 0

u(1, t) = 0

0 t 0,5,

u(x, 0) = f (x) = sin (x) + sin (2x)

0 x 1,

ut (x, 0) = g(x) = 0

0 x 1,

(308)

Por conveniencia tomaremos h = 0,1 y k = 0,05. Como c = 2 entonces

r=

2(0,05)
0,1

= 1. Como g(x) = 0 y r = 1, la f
ormula (306) para calcular los

valores de la segunda queda

ui,2 =

fi1 + fi+1
,
2
253

para i = 2, 3, 9

(309)

Sustituyendo r = 1 en la ecuaci
on (299) obtenemos la ecuaci
on en diferencias ya simplificada

ui,j+1 = ui+1,j + ui1,j ui,j1 )

(310)

Usando (309) y (310) generamos las aproximaciones a los valores de u(x, t)

que aparecen en la tabla de aproximaciones para 0 < xi < 1 y 0 tj 0,50
que coinciden en m
as de seis cifras decimales con los correspondientes a la
soluci
on exacta

u(x, t) = sin (x)cos(2t) + sin (2x)cos(4t)

Tabla de aproximaciones

254

(311)

tj

x2

x3

x4

x5

x6

x7

x8

x9

0,00

0,896802

1,538842

1,760074

1,538842

1,000000

0,363271

0,142040

0,05

0,769421

1,328438

1,538842

1,380037

0,951056

0,428980

0,000000

0,10

0,431636

0,769421

0,948401

0,951056

0,809017

0,587785

0,360616

0,1

0,15

0,000000

0,051599

0,181636

0,377381

0,587785

0,740653

0,769421

0,6

0,20

0,380037

0,587785

0,519421

0,181636

0,309017

0,769421

1,019421

0,9

0,25

0,587785

0,951056

0,951056

0,587785

0,000000

0,587785

0,951056

0,9

0,30

0,571020

0,951056

1,019421

0,769421

0,309017

0,181636

0,519421

0,5

0,35

0,363271

0,639384

0,769421

0,740653

0,587785

0,377381

0,181636

0,40

0,068364

0,181636

0,3606161

0,587785

0,809017

0,951056

0,948401

0,45

0,181636

0,210404

0,000000

0,428980

0,951056

0,380037

0,538842

0,50

0,278768

0,363271

0,142040

0,363271

1,000000

1,538842

1,760074

255

Ecuaciones Parabolicas.
Prof. Mara Angelica Vega U.
Modulo 15

17.

La ecuaci
on de conducci
on de calor.

Analogamente a la deduccion de la ecuacion de Laplace, se puede utilizar la conservacion del calor para desarrollar un balance de calor en una
barra larga aislada como en la siguiente figura Ademas de examinar el ca-

Figura 17: Barra aislada excepto en los

extremos.

235

so estacionario , este balance tambien considera la cantidad de calor que

se almacena en el elemento en el periodo 4t. El balance entonces es de la
forma:

Entradas - Salida = Acumulacion

es decir,

q(x)4y4z4t q(x + 4x)4y4z4t = 4x4y4zC4T

(273)

dividiendo por 4x4y4z y 4t se obtiene

q(x) q(x + 4x)
4T
= C
4x
4t

(274)

tomando lmite,

T
q
= C
x
t

(275)

Sustituyendo la ley de Fourier para la conduccion del calor, se obtiene

2T
T
=
x2
t

(276)

que es la ecuacion del calor.

En el caso de una ecuacion parabolica en dos variables independientes,
la malla esta abierta en los extremos en la dimension temporal.

236

6
6

yj+1

(xi , yj+1 )

yj

(xi , yj )

yj1

(xi , yj1 )

xi

(x0 , y0 )

xi1

xi+1 (xm+1 , y0 ) X

Analogamente al caso elptico, las ecuaciones parabolicas se resuelven

sustituyendo las derivadas parciales por diferencias finitas. La diferencia radica en que las ecuaciones elpticas estan acotadas en todas las dimensiones,
las ecuaciones parabolicas estan abiertas en la variable temporal

17.0.1.

M
etodos Explcitos

Usamos una diferencia centrada para aproximar la segunda derivada

espacial,
j
j
Ti+1
2Tij + Ti+1
2T
=
x2
42

(277)

que tiene un error O(4x2 ) y los superndices denotan el tiempo. La derivada

en tiempo la aproximamos por un diferencia progresiva:
T j+1 Tij
T
= i
t
4t
237

(278)

que tiene un error O(4t). Sustituyendo en la ecuacion (276) se obtiene

j
j
Ti+1
2Tij + Ti1
Tij+1 Tij
k
=
.
4x2
4t

(279)

j
j
Tij+1 = Tij + (Ti+1
2Tij + Ti1
)

(280)

de donde resulta

donde = k

4t
. Esta ecuacion permite calcular explcitamente los val(4x)2

ores de la funcion en cada nodo para un tiempo posterior j + 1 a partir

de la informacion en el nodo j y los nodos vecinos. Observemos que este
metodo es un analogo del metodo de Euler para resolver sistemas de ecuaciones diferenciales ordinarias. Un esquema de lo que ocurre en los nodos
esta representado en la siguiente figura,

Ejemplo 17.1 Usando em metodo explcito determine la distribuci

on de
temperatura en una barra larga y delgada que tiene una longitud de 10 cm
cal
y los siguientes valores k 0 = 0,49 [ segcm
C ], 4x = 2 cm y 4t = 0,1 seg.

En t = 0 la temperatura de la barra es cero y las condiciones de frontera

en T (0) = 100 C y T (10) = 50 C. Considere que la barra es de aluminio
0,49
g
con C = 0,2174 [ gcal
=
C ] y = 2,7 cm3 . Por lo tanto, k =
2,7 0,2174


cm2
0,835
y = 0,835 0,1
= 0,020875
22
seg
238

Figura 18: Metodo Explcito

Sol. Aplicando la ecuaci

on (280) se obtiene el siguiente valor en t = 0,1
seg para el nodo en x = 2 cm
T11 = 0 + 0,020875(0 2(0) + 100) = 2,0875
En los otros nodos interiores, x = 4,6 y 8 cm, los resultados son
T21 = 0 + 0,020875(0 2(0) + 0) = 0
T31 = 0 + 0,020875(0 2(0) + 0) = 0
T21 = 0 + 0,020875(50 2(0) + 0) = 1,0438

239

En t = 0,2 se obtiene para los cuatro nodos interiores

T12 = 2,0875 + 0,020875(0 2(2,0875) + 100) = 4,0878
T22 = 0 + 0,020875(0 2(0) + 2,0875) = 0,043577
T32 = 0 + 0,020875(1,0438 2(0) + 0) = 0,021788
T42 = 1,0438 + 0,020875(50 2(1,0438) + 0) = 2,0439
El aumento general de la temperatura con el tiempo indica que la difusion
del calor desde las fronteras al interior de la barra.

17.0.2.

Convergencia y estabilidad

La Convergencia significa que si 4x y 4y tienden a cero, entonces los

valores obtenidos por el metodo se aproximan a la solucion exacta.
La Estabilidad significa que los errores en cualquier etapa del calculo, no
se amplifican en etapas posteriores sino disminuyen.
Un resultado importante es que el meetodo explcito es convergente y
estable si

1
2

o
4t

17.0.3.

1 4x2
2 k

(281)

Aproximaciones temporales de orden superior

La idea general es volver expresar la EDP como un sistema de EDO se

denomina m
etodo de lneas. En efecto, una manera de mejorar el metodo
240

de Euler usado antes sera emplear un esquema de integracion mas exacto

para resolver las EDO. Por ejemplo, el metodo de Heun puede utilizarse
para obtener una exactitud temporal de segundo orden.

17.1.

Un m
etodo implcito simple

Las formulas explcitas por diferencias finitas tienen problemas con la

estabilidad. Los metodos implcitos son mas estables, pero usan algoritmos
mas complicados.
La diferencia fundamental entre ambos metodos radica en la forma de
tomar los nodos, se ilustra en la siguiente figura

Figura 19: Puntos malla

241

En un esquema explcito, aproximamos la derivada espacial para un nivel

de tiempo j, (figura a) en cambio en los esquemas implcitos la derivada
espacial se aproxima en un nivel de tiempo posterior j + 1 (figura b). Por
ejemplo la segunda derivada se aproximara por
j+1
j+1
Ti+1
2Tij+1 + Ti1
2T
=
x2
(4x)2

(282)

que tiene una exactitud de segundo orden. Al sustituir en la ecuacion original, se obtiene
j+1
j+1
Ti+1
2Tij+1 + Ti1
Tij+1 Tij
k
=
(4x)2
4t

(283)

que se puede escribir como

j+1
j+1
Ti1
+ (1 + 2)Tij+1 Ti+1
= Tij

donde =

(284)

k4t
. Esta ecuacion se aplica a todos los nodos excepto al
(4x)2

primero y al u
ltimo de los nodos interiores, los cuales se modifican al considerar las condiciones de frontera. La condicion de frontera en el extremo
izquierdo de la barra (i = 0) se expresa:
T0j+1 = f0 (tj+1 )

(285)

donde f0 (tj+1 ) es una funcion que describe como cambia con el tiempo
la temperatura de la frontera. Sustituyendo (285) en (284) se obtiene la
ecuacion en diferencias para el primer nodo interior (i = 1)
(1 + 2)T1j+1 T2j+1 = T1j + f0 (tj+1 )
242

(286)

Analogamente para (i = m)
j+1
j+1
j
+ fm+1 (tj+1 )
Tm1
+ (1 + 2)Tm
= Tm

(287)

donde fm+1 (tj+1 ) describe los cambios especficos de temperatura en el extremo derecho de la barra (i = m + 1). observamos que el sistema resultante
de m ecuaciones lineales con m incognitas es tridiagonal.

Ejemplo 17.2 Resolver el ejemplo anterior usando diferencias finitas implcitas. Para este caso = 0,020875 .
para t = 0 la ecuaci
on para el primer nodo interior es
1,04175T11 0,020875T21 = 0 + 0,020875(100)
o
1,04175T11 0,020875T21 = 2,0875
An
alogamente,

1
T1

0,020875
0
0
2,0875
1,04175

0,020875
1
1,04175
0,020875
0
T2 0

0
0,020875
1,04175
0,020875 T3 0

1
1,04375
T4
0
0
0,020875
1,04175

243

En t = 0,1 seg
T11 = 2,0047
T21 = 0,0406
T31 = 0,0209
T41 = 1,0023
En t = 0,2 seg
T12 = 3,9305
T22 = 0,1190
T32 = 0,0618
T42 = 1,9653
El metodo implcito descrito es estable y convergente, pero tiene la desventaja que la aproximacion en diferencias temporal tiene una exactitud de primer
orden y con diferencia espacial de oreden 2. Hay que mencionar que aunque
este metodo es incondicionalmente estable hay un lmite de exactitud para
el uso de pasos tiempo grande. Por tanto no es mucho mas eficiente que los
metodos explcitos para la mayora de los problemas variables en el tiempo.

17.2.

El m
etodo de Crank-Nicholson

El metodo de Crank-Nicholson es un esquema implcito alternativo que

tiene una exactitud de segundo orden, tanto para la aproximacion espacial,

244

como para la temporal. Para alcanzar esta exactitud en el tiempo se desarrollan diferencias en el punto medio del incremento tiempo, por lo tanto la
1

primera derivada temporal se aproxima en tj+ 2 por

T j+1 Tij
T
= i
t
4t

(288)

La segunda derivada en el espacio puede determinarse en el punto medio

promediando las aproximaciones por diferencias al principio tj y al final
tj+1 del incremento del tiempo.
" j
#
j
j
j+1
j+1
j+1
1 Ti+1 2Ti + Ti1 Ti+1 2Ti + Ti1
2T
=
+
x2
2
(4x)2
(4x)2

(289)

sustituyendo las ecuaciones (288) y (289) en la ecuacion wwwwwww , se

obtiene
j+1
j+1
j
j
Ti+1
+ 2(1 + )Tij+1 Ti1
= Ti+1
+ 2(1 )Tij + Ti1

donde =

(290)

k4t
. Como en el caso del metodo implcito simple, se determi(4x)2

j+1
nan las condiciones de fronteea T0j+1 = f0 (tj+1 ) y Tm+1
= fm+1 (tj+1 ) para

obtene las ecuaciones para los nodos interiores primero y u

ltim o.
Para el primer nodo interior se obtiene
2(1 + )T1j+1 )T2j+1 = f0 (tj ) + 2(1 )T1j + T2j + f0 (tj+1 )

245

(291)

y para el u
ltimo nodo interior
j+1
j
j+1
j
Tm1
+ 2(1 + )Tm
= fm+1 (tj ) + 2(1 )Tm
+ Tm1
+ fm+1 (tj+1 )

(292)
observamos que los sistemas de ecuaciones involucrados son tridiagonales,
lo que hace eficiente el proceso.

Ejemplo 17.3 Resolver el ejemplo anterior, usando el metodo de CrankNicholson.

1
T 1

0,020875
0
0
4,175
2,04175

0,020875

1
2,04175
0,020875
0

T 2 0

0
0,020875
2,04175
0,020875 T3 0

1
0
0
0,020875
2,04175
T4
2,0875
En t = 0,1 seg
T11 = 2,0450
T21 = 0,0210
T31 = 0,0107
T41 = 1,0225

246

De las ecuaciones se obtiene

T12 = 4,0073
T22 = 0,0826
T32 = 0,0422
T42 = 2,0036

247