Mtodos Numricos y Programacin

Braulio Gutirrez Pari Universidad Peruana Unin-Juliaca Facultad de Ingeniera y Arquitectura Ingenieria Civil Marzo del 2013 Chullunquiani-Per

ii

ndice general
Introduccin 1. Preliminares del Matlab 1.1. Notacin Matricial . . . . . . . . . . . . . . . 1.2. Operaciones con vectores y matrices . . . . . . 1.3. Funciones que actan sobre matrices . . . . . 1.4. Edicin de linea . . . . . . . . . . . . . . . . . 1.5. Fichero function . . . . . . . . . . . . . . . . . 1.5.1. Bucles . . . . . . . . . . . . . . . . . . 1.5.2. Sentencias condicionales . . . . . . . . 1.6. Algunas pautas para gracar funciones simples
V

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

1 1 2 3 3 8 8 8 16 19 19 20 21 22 26 29 29 31 32 36 37 41 46 50 50 53 54 55 55 56 57 57 58

2. Nociones Bsicas Sobre Errores 2.1. Factores que intervienen en la resolucin numrica . . 2.1.1. Precisin de los datos de entrada . . . . . . . 2.1.2. Representacin de los datos en el computador 2.1.3. Las operaciones numricas efectuadas . . . . . 2.2. Errores Absolutos y Errores Relativos . . . . . . . . . 3. Ceros reales de funciones reales 3.1. Aislamiento de las Raices . . . . . . . . . . . . 3.2. Renamiento: Mtodos Iterativos . . . . . . . 3.2.1. Mtodo de Biseccin . . . . . . . . . . 3.2.2. Mtodo de la Posicin Falsa . . . . . . 3.2.3. Mtodo del Punto jo (MPF) . . . . . 3.2.4. El Mtodo de Newton-Raphson . . . . 3.2.5. Mtodo de Secante . . . . . . . . . . . 3.2.6. Comparacin de los Mtodos Iterativos 3.3. Problemas . . . . . . . . . . . . . . . . . . . . 4. Resolucin de Sistemas Lineales 4.1. Aspectos Tericos . . . . . . . . . . . . . . . . 4.1.1. Conjunto Imagen de una Matriz A . . 4.1.2. Rango de una Matriz A . . . . . . . . 4.2. Mtodos Numricos para Resolver Sistemas de 4.3. Mtodo de Cramer . . . . . . . . . . . . . . . 4.4. Mtodo de Gauss . . . . . . . . . . . . . . . . 4.4.1. Resolviendo un Sistema Triangular . . iii

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ecuaciones Lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

iv 4.4.2. Conversin de un Sistema de Ecuaciones Lineales a uno perior . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5. Sensibilidad en Sistemas Lineales: Nmero de Condicin . . . 4.6. Estrategia con Pivoteamiento Parcial . . . . . . . . . . . . . . 4.7. Descomposicin LU . . . . . . . . . . . . . . . . . . . . . . . . 4.7.1. Clculo de los Factores L y U . . . . . . . . . . . . . . 4.8. Descomposicin de Cholesky . . . . . . . . . . . . . . . . . . . 4.8.1. Clculo del factor de Cholesky . . . . . . . . . . . . . . 4.9. Mtodos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . 4.10. Comparacin entre Mtodos Directos e Iterativos . . . . . . .

NDICE GENERAL Triangular Su. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59 60 62 63 64 66 66 69 71 73 74 75 77 77 78 79 82 84 89 89 91 93

5. Introduccin a Sistemas No Lineales 5.1. El Mtodo de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1. Convergencia y Rapidez del Mtodo de Newton . . . . . . . . . . . . . 6. Interpolacin y Aproximacin 6.1. Interpolacin . . . . . . . . . . . . . . . . 6.1.1. Interpolacin Polinomial . . . . . . 6.1.2. Modos de Obtener el Polinomio pn 6.2. Aproximacin: Caso Discreto . . . . . . . . 6.2.1. Mtodo de Mnimos cuadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7. Integracin numrica 7.1. Mtodo de los Trapecios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2. Mtodo de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliografa

Introduccin
La matemtica est comprendida de varias partes, cada una de ellas tiene importancia propia, pero algunas tambin son fundamentales en diferentes disciplinas. Muchos problemas de la vida real pueden ser representados por formulaciones matemticas, las cuales son denominadas modelos matemticos. Usando argumentos matemticos tericos, algunas veces es posible garantizar la existencia de soluciones para esos modelos matemticos, pero encontrar manualmente esas soluciones puede resultar extremadamente difcil y a veces imposible. Estudiar mtodos numricos desde un punto de vista general nos permitir analizar mecanismos de clculo capaces de otorgar soluciones, o aproximaciones a las soluciones, all donde las herramientas tericas fracasaban. Estos mecanismos numricos de clculo deben caminar de la mano con el computador, pues en su mayora requieren de muchos pasos y frecuentemente estn orientados a la resolucin de problemas de grandes dimensiones. La importancia del estudio de mtodos numricos es indiscutible, pues gran parte de las investigaciones en ciencias, ingenieras, economa, etc., recurren a tcnicas numricas para la obtencin de resultados.

vi

INTRODUCCIN

Captulo 1 Preliminares del Matlab

1.1. Notacin Matricial

Sea Rmn el espacio vectorial de las matrices reales de dimensin m n, de la forma A= a1,1 a2,1 . . . am,1 a1,2 a1,n a2,2 a2,n . . ... . . . . am,2 am,n

donde la matriz anterior se denota como A = (ai,j ), i = 1, ..., m y j = 1, ..., n Los elementos de una matriz se introducen entre corchetes. Las las se separan mediante un punto y coma (;) y los elementos separados por los espacios en blanco o comas. >> A=[1 2 3; 3 1 2; 1 1 0] Una vez denida una matriz o un vector, se pede acceder a sus elementos o submatrices con las rdenes

Orden V (i) V (end) A(i,j ) A(:,j ) A(i,:) A(:) A(v ,w)

Salida Coordedenada i del vector V ltima coordenada del vector V Elemento de la matriz A que ocupa la posicin i,j Columna j de la matriz A Fila i de la matriz A Todos los elementos de A, en una sola columna Submatriz de A que contiene las las indicadas en las coordenadas de v y las columnas indicadas en w A(i,:) = [ ] Elimina la la i de la matriz A A(:,j ) = [ ] Elimina la columna j de la matriz A Puede denirse tambin ciertas matrices con las siguientes rdenes 1

CAPTULO 1. PRELIMINARES DEL MATLAB

Orden ones(n) ones(m,n) zeros(n) zeros(m,n) eye(n) eye(m,n)

Salida Matriz cuadrada n n de unos Matriz m n de unos Matriz cuadrada n n de ceros Matriz m n de ceros Matriz identidad n n Matriz m n con unos en la diagonal principal y el resto de ceros

1.2.

Operaciones con vectores y matrices

Si A y B son matrices con las dimensiones adecuadas y es un escalar, las operaciones habituales se efectan con las siguientes rdenes

Operacin A+B AB AB A/B A\B A n

Resultado Suma A y B Resta A y B Multiplica A por B Divisin derecha matricial Divisn izquierda matricial potencia matricial

Ademas de las operaciones mencionadas, en MATLAB se denen otras operaciones a las que llamaremos operaciones elemento a elemento

Operacin A. B A./B A.\B A.^B A0 A.0

Resultado Multiplicacin de arreglos Divisin derecha de arreglos Divisin izquierda de arreglos Potenciacin de arreglos transpuesta compleja transpuesta matricial

1.3. FUNCIONES QUE ACTAN SOBRE MATRICES

1.3.

Funciones que actan sobre matrices

Las siguientes funciones permiten obtener informacin sobre las matrices o vectores que tienen como argumentos de entrada Funciones size(A) size(A, 1) size(A, 2) length(V ) det(A) trace(A) inv(A) sum(A) sum(V ) diag(A) Salida Vector con las dimensiones de la matriz A Nmero de las de la matriz A Nmero de columnas de la matriz A Nmero de coordenadas del vector V Determinante de la matriz A Traza de la matriz A Devuelve la inversa de A Devuelve un vector la en el que el elemento i contiene la suma de todos los elementos de la columna i de A. Devuelve la suma de coordenadas del vector V Devuelve la diagonal de la matriz A.

1.4.
, ; Esc % clc

Edicin de linea
Contina la entrada de la expresin en la lnea siguiente Permite separar instrucciones en la misma ejecucin Ejecuta la instruccin pero no muestra el resultado Borra la lnea Todo lo que aparece en la lnea despus del smbolo % se considera un comentario Borra el contenido de la ventana de comandos y lo coloca el prompt (>>) en la pri mera lnea de la ventana

Obviamente, las funciones matemticas habituales tambin estn en MATLAB con la nica particularidad de que actan sobre vectores a matrices elemento a elemento Funcin ex ln(x) log10 (x) log 2 (x) x |x| MATLAB exp(x) log(x) log10(x) log2(x) sqrt(x) abs(x) Funcin sen(x) cos(x) tan(x) cot(x) sec(x) csc(x) MATLAB sen(x) cos(x) tan(x) cot(x) sec(x) csc(x)

Comandos que se usan en MATLAB == > < >= <= = igual a mayor que menor que mayor o igual que menor o igual que distinto de

4 y los siguientes conectivos lgicos

CAPTULO 1. PRELIMINARES DEL MATLAB

& y | o no xor o excluyente Ejemplo 1.1 1 4 7 3 A= 2 5 1 5 3 6 5 1

>> A=[1 4 7 3;2 5 1 5;3 6 5 1] A= 1 4 7 3 2 5 1 5 3 6 5 1

>> A(3,2)= O % reemplaza la tercera fila y la segunda columna con cero A= 1 4 7 3 2 5 1 5 3 0 5 1 >> A(4,3)= 8 % A(4,3) no existe; se modifica las dimensiones de la matriz para aadir el elemento A(4,3) en una cuarta fila. El resto se llena con ceros. A= 1 2 3 0 4 5 0 0 7 1 5 8 3 5 1 0

>> A(1,5)= 9 % A(1,5) no existe; se modifica las dimensiones de la matriz para aadir el elemento A(1,5) en una quinta columna. El resto se llena con ceros A= 1 2 3 0 4 5 0 0 7 1 5 8 3 5 1 0 9 0 0 0

>> A(:,2) =[] % elimina la segunda columna de A A= 1 2 3 0 7 1 5 8 3 5 1 0 9 0 0 0

1.4. EDICIN DE LINEA A=[A,[1:4]] % aade una columna a la matriz A A= 1 2 3 0 7 1 5 8 3 5 1 0 9 0 0 0 1 2 3 4 1 2 3 4 A= 5 6 7 8 9 0 5 2

Ejemplo 1.2

>> A=[1 2 3 4;5 6 7 8;9 0 5 2] A= 1 2 3 4 5 6 7 8 9 0 5 2

>> [m n] =size(A) % tamao de la matriz m= 3 n= 4 Cmo seleccionar a los elemento de una matriz >> A(2,3) % elemento de la segunda fila y tercera columna de A ans= 7 >> A(3,:) % elemento de la tercera fila de A ans= 9 0 5 2 >> A(:,2) % elemento de la segunda columna de A ans= 2 6 0 Cmo acceder a una submatriz de una matriz dada >>A([2,3],:) % submatriz de A formada por las filas de 2 y 3 ans= 5 6 7 8 9 0 5 2

CAPTULO 1. PRELIMINARES DEL MATLAB >> A(:,[2,3,4]) % submatriz de A formada por las columnas 2,3 y 4 ans= 2 3 4 6 7 8 0 5 2 >> A([2,3],[1,2]) % submatriz formado por filas 2 y 3, columnas 1 y 2 ans= 5 6 9 0 >> A(2:3,2:4) % submatriz de A formada por filas (2 a 3) y columnas (2 a 4) ans=

6 7 8 0 5 2 Las matrices formadas por una sola la o columna se denominavectores la y columna respectivamente. Las las y columnas de una matriz se pueden completar utilizando vectores la y columna. >> b=[2 4 6] % vector columna b= 2 4 6 >> b(3) % tercer elemento del vector ans= 6 >> [A b] % agregamos una columna a la matriz A ans= 1 2 3 4 2 5 6 7 8 4 9 0 5 2 6 >> v=[1 3 5 7] % vector fila V= 1 3 5 7 >> v(4) % cuarto elemento del vector ans= 7 >> [A v] % agregamos una fila a la matriz A

1.4. EDICIN DE LINEA ans= 1 5 9 1 2 6 0 3 3 7 5 5 4 8 2 7 2 4 B= 6 8

Ejercicio 1.1 Sean los arreglos 1 2 A = 3 1 3 C = 4 7 0 5 8 2 6 9 1 5 7 D= Solucin: >> D=[A. B;C] D= -1 2 3 4 -5 6 -7 8 9 0 2 -1 1 3 5 b) Construya el arreglo E= Solucin: >> E=[[A; C] B] E= -1 4 -7 0 2 2 -5 8 2 4 3 6 9 -1 6 1 3 5 7 8 Ejercicio 1.2 Construya una matriz mgica de orden 7 (A=magic(7)) y efecte las siguientes operaciones 1. Obtenga en un arreglo P los elementos de A comprendidos entre las las 2 y 5 y las columnas 1 y 4. 2. Obtenga en un arreglo Q las tres ltimas columnas de A. 3. Obtenga en un arreglo R las tres primeras las de A 4. Crear un arreglo B que contenga las las de A con las las 1 y 4 intercambiadas. 5. Incrementar la la 4 del arreglo B en 5 veces la la 7 2 4 6 8 7 A B C At C B

a) Construya el arreglo

CAPTULO 1. PRELIMINARES DEL MATLAB 6. Asignar a las columnas 3 y 6 de A, las las 2 y 4 del arreglo B respectivamente 7. Eliminar la la 3 y la columna 5 del arreglo B. 8. Intercambiar las columnas 1 y 7 del arreglo A. 9. Listar los elementos del arreglo A como un nico vector columna.

1.5.

Fichero function

Dentro de la organizacin de un programa es muy comn la realizacin de tareas que pueden servir para diferentes programas o simplemente la separacin en etapas del programa global que se pueden abordar independientemente. Una de las formas de realizar esta divisin en MATLAB es la travs de las function. Las function se construyen en cheros m. function [Argumentos de salida]= nombre_funcin (Arg. Entrada) A continuacin daremos los comandos ms usados en programas

1.5.1.

Bucles

son iteraciones que se utilizan para controlar el ujo de un programa for-end . La sintaxis de este comando es for i=vi:inc:vf instrucciones end While-end . La sintaxis de este comando es while relacin instrucciones end

1.5.2.

Sentencias condicionales
if relacin instrucciones end

if-end . La sintaxis de este comando es

if-else-end if relacin instruccin 1 else instruccin 2 end;

1.5. FICHERO FUNCTION

Ejemplo 1.3 Construya un programa donde al introducir un nmero n,determine si pertenece o no al intervalo [a, b] Solucin: function Intervalo_1(a,b,x) if x>=a & x<=b disp(esta dentro del intervalo) else disp(esta fuera del intervalo) end Ejemplo 1.4 Construir una funcin raices.m que calcule las raices de un polinomio de segundo grado ax2 + bx + c = 0 y calcule las raices de x2 3x + 2 = 0 Solucin: function [x1,x2]= raices(a,b,c) D = b^2-4*a*c x1=(-b + sqrt(D))/(2*a); x2=(-b - sqrt(D))/(2*a); Ejemplo 1.5 Construya un programa donde puedas calcular el rea y el volumen de un cilindro y calcule para r=5, h=12 Solucin: function [area,volumen]= cilindro(r,h) area=2*pi*r*h+2*pi*r^2; volumen=pi*r^2*h; La sentencia for-end Los bucles que utiliza for en MATLAB tiene el siguiente aspecto for contador = Valor inicial : incremento : valor nal sentencia del cuerpo del bucle end Ejercicio 1.3 A n de aclarar las ideas ingrese los siguientes comandos y comente sus resultados for i=1:8 i end Ejercicio 1.4 Ingrese los siguientes comandos y comente sus resultados for i=1:3:19 i end Ejemplo 1.6 hacer un programa donde puedas calcular la suma de los n primeros nmeros naturales

10 Solucin:

CAPTULO 1. PRELIMINARES DEL MATLAB

function s=Sumita(n) suma=0; for x=1:n, suma = suma + x; end; s=suma; Ejemplo 1.7 Hga un programa que reciba un vector y entregue el promedio del valor de sus elementos Solucin: function p=promedio(x) n=length(x); suma=0; for i=1:n suma=suma+x(i); end p=suma/n; Ejemplo 1.8 Construya un programa que evalue la funcin f (x) = Solucin: function y= f1(x) n=length(x); for i=1:n if x(i)<=0 y(i)=2*(sen(2*x(i)))^2; else y(i)=1-exp(-x(i)); end end x 1, 1 x2 , Ejemplo 1.9 Construya un programa que evalue la funcin: f (x) = 1 , x+1 Solucin: function y= f2(x) n=length(x); for i=1:n if x(i)<=-2 y(i)=x(i)-1; elseif (x(i)>-2 & x(i)<0) y(i)=1-x(i).^2; else y(i)=1/(x(i)+1); end end x 2. 2 < x < 0, x0 2 sen2 (2x), 1 ex , x0 x>0

1.5. FICHERO FUNCTION Ejemplo 1.10 Escriba una funcin para elegir el mayor entre dos nmeros Solucin: function M=mayor(a,b) if a>b M=a; else M=b; end

11

Ejemplo 1.11 Escriba una funcin que dada dos vectores A, B y encuentre los elementos de la interseccin de ambos vectores Solucin: function C= interseccion(A,B) m=length(A); n=length(B); k=1; for i=1:m for j=1:n if A(i)==B(j) C(k)=A(i); k=k+1; end end end Ejemplo 1.12 Haga un programa que me calcule el valor absoluto Solucin function va= ValorAbsoluto(n) if n >= 0 va = n; else va = -n; end

Ejemplo 1.13 Construir una funcin raices2.m que, que muestre el polinomio de segundo grado ax2 + bx + c = 0 si tiene raices complejas o reales Solucin

12

CAPTULO 1. PRELIMINARES DEL MATLAB

function [x1,x2]= raices1(a,b,c) D = b^2-4*a*c; for D>=0 disp (Raices Reales) x1=(-b + sqrt(D))/(2*a); x2=(-b - sqrt(D))/(2*a); else disp (Raices Complejas) x1=(-b + sqrt(D))/(2*a); x2=(-b - sqrt(D))/(2*a); end Ejemplo 1.14 Crear un programa tal que dada la matriz A cuadrada de dimensin n devuelva una nueva matriz B cuyos elementos son elementos de A salvo los de la diagonal principal que son ceros. Solucin: function B= Inters(A) [m n]=size(A); if (m == n) error (la matriz no es cuadrada) end for i=1:m for i=1:n if i=j B(i,j)=A(i,j); else B(i,j)=0; end end end Ejemplo 1.15 Haga un programa que calcule la traza de una matriz cuadrada.A de dimensin n Solucin: function suma= Traza(A) [m n]=size(A); if (m == n) error (la matriz no es cuadrada) end suma=0; for i=1:m for j=1:n if i==j suma=suma+A(i,j); end end end

1.5. FICHERO FUNCTION Ejercicio 1.5 construir un programa que me calcule la norma euclidiana en Rn Solucin: function N= norma(v) n=length(v); suma=0; k=1; while k <=n suma=suma+v(k)^2; k=k+1; end N=sqrt(suma); Ejemplo 1.16 Haga un programa donde dado dos matrices cuadradas se calcule la suma Solucin: function C= SumMatrices(A,B) [m1,n1]=size(A); [m2,n2]=size(B); if (m1=m2) | (n1=n2) error (error en las dimensiones de las matrices) end for i=1:m1 for j=1:n1 C(i,j)=A(i,j)+B(i,j); end end

13

Ejemplo 1.17 Haga un programa que multiplique una matriz con un vector controlando las dimensiones Solucin: Sea la matrices A R
mn n X y x R , el producto es yi = ai,j xj n j =1

function y= PMatrizVector(A,x) [m n]=size(A); p=size(x,1); if (p=n) error (error en las dimensiones) end y=zeros(m,1); for i=1:m for j=1:n y(i)=y(i)+ A(i,j)*x(j); end end Ejemplo 1.18 Haga un programa que multiplique una matriz cuadrada triangular superior son un vector columna controlando las dimensiones

14 Solucin:

CAPTULO 1. PRELIMINARES DEL MATLAB

function k= PTriangSupVect(A,v) n=size(A,1); p=size(v,1); k=zeros(n,1); if (p=n) error (error en las dimensiones) end for i=1:n k(i)=A(i,i:n)*v(i:n); end Ejemplo 1.19 Haga un programa que multiplique dos matrices cuadradas triangulares y superiores, controlando las dimensiones Solucin: function r= ProdTriangSUp(A,B) n=size(A,1); p=size(B,1); if (p=n) error (error en las dimensiones) end r=zeros(n); for i=1:n for j=1:n r(i,j)=A(i,i:j)*B(i:j,j); end end Ejemplo 1.20 Crear un programa que efectue el producto punto o interior de dos vectores Solucin: function v= Producto_Punto(x,y) n=length(x); m=length(y); if (n=m) error (error en las dimensiones de los vectores) end v=0; for i=1:n v=v+x(i)*y(i); end Ejemplo 1.21 Haga un programa donde dado dos matrices cuadradas se calcule el producto Solucin: Dada las matrices A Rmp y B Rpn , el producto es otra matriz P = AB

1.5. FICHERO FUNCTION Rmn

15

function P= producto(A,B) [m,p]=size(A); [q,n]=size(B); if (p=q) error (producto matricial incompatible) end for i=1:m for j=1:n S=0; for k=1:p S=S+A(i,k)*B(k,j); end P(i.j)=S end end

otra forma function P= producto(A,B) [m,p]=size(A); [q,n]=size(B); if (p=q) error (producto matricial incompatible) end C=zeros(m,n): for i=1:m for j=1:n C(i,j)=A(i,:)*B(:,j); end end Ejemplo 1.22 describir la implementacin en un archivo M funcin, que permita evaluar la funcin 2y x2 f (x, y ) = x + p x2 + y 2 + 103 y evale en los puntos Ejercicio 1.6 1. (x, y ) = (3, 4) 2. (x, y ) = (, 2) function z=fun(x,y) z=x+(2*y-x^2)/sqrt(x^2+y^2+1e-3); >> z=fun(3,4) z= 2.8000

Solucin:

16 >> z=fun(pi,sqrt(2)) z= 1.0979

CAPTULO 1. PRELIMINARES DEL MATLAB

Ejemplo 1.23 Realizar un programa que permita pasar de grados farentheit a Celsius Solucin function C= Celsius(o F) C=5/9*(o F-32);

1.6.

Algunas pautas para gracar funciones simples

sen x x

Gracar una funcin f : R 7 R sobre un intervalo [a, b] usando el comando fplot Graquemos la funcin f : R 7 R sobre un intervalo [8, 10] , donde f (x) = procedemos a invocar en la ventana de comandos >> fplot(sin(x)/x,[-8,10]) Obtenemos el siguiente grco

La orden de dibujo ms simple es plot Ejemplo 1.24 Supongamos que conocemos la hora en que se conoce la temperatura en una estacin meteorolgica 9 12 15 18 21 24 h T 15.7 17.2. 19.5 16.7 15.5 12.6 En la ventana de comandos ejecutamos.

1.6. ALGUNAS PAUTAS PARA GRAFICAR FUNCIONES SIMPLES >> h=[4 8 12 16 20 24]; >> T=[15.6 16.2 18.5 17.1 16.5 15.3]; >> plot(h,T)

17

18

CAPTULO 1. PRELIMINARES DEL MATLAB

Captulo 2 Nociones Bsicas Sobre Errores

En este primer captulo an no estudiaremos los mtodos numricos, sino que veremos algunos conceptos bsicos y sealaremos algunos factores que intervendrn en la resolucin de problemas mediante el computador. Existen diversas fases cuando intentamos resolver un problema mediante mtodos numricos. Puede suceder que los resultados nales obtenidos no sean justamente los esperados, aunque todas las fases hayan sido ejecutadas correctamente, los motivos pueden ser varios y los estudiaremos a continuacin.

2.1.

Factores que intervienen en la resolucin numrica

La mayora de los mtodos numricos que veremos aqu tienen un carcter iterativo, esto signica que iniciarn con una estimativa de la solucin, digamos x(0) , para luego (0) (1) inicial (2) construir una sucesin de valores x , x , x , ... de modo que l mk x(k) = x , donde x es la solucin. Por lo general, x(k+1) ser calculada a partir de x(k) mediante un procedimiento debidamente fundamentado. Frecuentemente, este procedimiento requerir un gran nmero de clculos, por lo que ser necesario el auxilio de un computador. Lamentablemente, los computadores presentan algunos inconvenientes que, si no se controlan, pueden ocasionar respuestas catastrcas. Ejemplo 2.1 Considere el trecho de un programa en Matlab, tal como se muestra en la gura adjunta. Observe que en teora, debera cobrarse S/,1000000. Programa principal a=1; b=0.0000000000000001; c=1; if a+b>c disp(cobrar S/.1000000); else disp(pagar S/. 2000000); end Un error numrico grave cometido por el computador. Sin embargo, debido al error cometido por el computador, se terminara pagando S/.2 000 000. 19

20

CAPTULO 2. NOCIONES BSICAS SOBRE ERRORES

Laboratorio 2.1 Haga un programa en algn lenguaje de programacin que usted conozca, de modo que en la prctica corrobore el Ejemplo anterior. Al resolver un problema por mtodos numricos, los resultados obtenidos pueden depender de: 1. La precisin de los datos de entrada 2. La forma cmo stos son representados en el computador 3. Las operaciones numricas efectuadas Cada uno de estos temas sern explicados a continuacin.

2.1.1.

Precisin de los datos de entrada

Cuando ingresamos los datos de un problema, ya sea para calcular en un papel o trabajar en el computador, ellos contienen una imprecisin inherente, quiere decir que no hay cmo evitarlos. El siguiente ejemplo aclara esta armacin. Ejemplo 2.2 Sabemos que para calcular el rea de un crculo, tenemos que ingresar numricamente el radio r y el valor de . El valor de r quiz pueda ser conocido exactamente (r = 2), pero apenas podemos conocer una aproximacin de con un nmero nito de dgitos. As, aproximando el valor de por 3,14, el rea del crculo ser: 3, 14 (2)2 = 12, 56 m2 Si consideramos 3, 1416, entonces el rea del crculo estar dada por: 3, 1416 (2)2 = 12, 5664 m2 Pero si consideramos 3, 141592654, entonces el rea ser 3, 141592654 (2)2 = 12, 566370616 m2 Claramente las imprecisiones de los datos de entrada ocasionan imprecisiones en los resultados. En el ejemplo anterior vimos que el mejor resultado se obtuvo en el ltimo caso. Pero, cuando usamos un computador, cuntos dgitos decimales reconoce ste? El siguiente ejemplo intentar aclarar esta situacin. Ejemplo 2.3 Usando Matlab 6.0 en un PC de 32-bit con un sistema operativo Windows XP, hicimos la siguiente operacin: 0,00000000000001 + 1 y el resultado obtenido fue 1,00000000000001, lo cual es satisfactorio. Luego, hicimos: 0,000000000000001 + 1 y obtuvimos como respuesta 1. Qu sucedi?

2.1. FACTORES QUE INTERVIENEN EN LA RESOLUCIN NUMRICA

21

El ejemplo anterior nos hace ver que en el primer caso, cuando se usan los 14 dgitos a la derecha del punto, el sistema de cmputo no comete error. Mientras que, si se usan 15 dgitos, el sistema nos otorga una respuesta errada. La razn se debe a que todo computador trabaja con un nmero nito y bien reducido de dgitos, en nuestro caso 14 a la derecha del punto, si el nmero de dgitos sobrepasa lo esperado, el sistema lo trunca o redondea, dependiendo del sistema utilizado. Laboratorio 2.2 Utilice un paquete o lenguaje de programacin en un computador, para corroborar los resultados obtenidos en el Ejemplo 1.3.

2.1.2.

Representacin de los datos en el computador

Un nmero con representacin decimal nita puede tener una representacin innita en el sistema binario. Como un computador trabaja con el sistema binario y con una cantidad ja de dgitos, necesariamente trabajar con una aproximacin, por lo tanto no se obtendrn resultados exactos. El sistema con el que trabajamos comunmente es el decimal, un nmero x en este sistema lo representaremos algunas veces, cuando se preste a confusin, por (x)10 . Por otro lado, el sistema con el que trabaja un computador hoy en da es el binario, un nmero y en este sistema ser representado por (y )2 . As por ejemplo. (5)10 y (101)2

representa el nmero 5 en el sistema decimal y binario, respectivamente. El siguiente ejemplo muestra cmo esta representacin aparentemente inofensiva, puede generar terribles errores cuando se trabaja en el computador. Ejemplo 2.4 Considere la siguiente sumatoria: S= Para ai = 0, 5, el resultado exacto debera ser S = 0, 5
30000 X i=1 30000 X i=1

ai

1 = 0, 5 30000 = 15000

Despus de implementar un pequeo programa en computador, el resultado fue tambin 15000. Claramente, no hay por qu preocuparse en este caso, los resultados son los mismos. Pero, para ai = (0, 11)10 , el resultado exacto debera ser S = (0, 11)10 30000 = 3300 Sin embargo, el resultado obtenido por el computador fue S = 3300, 00000000063 cmo explicar la diferencia de resultados en este caso?

22

CAPTULO 2. NOCIONES BSICAS SOBRE ERRORES

Esto se debe a que el nmero (0, 11)10 , cuya representacin en el sistema decimal es nito, tiene una representacin binaria innita: (0, 0001110000101000111101)2 Si el computador trabajara con 14 dgitos despus del punto, el nmero debera ser cortado o redondeado, lo cual representa ya un error. Todos los clculos subsiguientes sern, afectadas por este hecho. Laboratorio 2.3 Utilizando algn lenguaje de programacin, haga un programa para ejecutar lo tratado en el ejemplo 1.4. Por lo general, el error ocurrido depende de la representacin de los nmeros en la mquina utilizada. La representacin de un nmero depende de la base elegida o disponible en la mquina en uso, y, del nmero mximo de dgitos usados en su representacin. Cualquier clculo que envuelva nmeros que no pueden ser representados a travs de un nmero nito de dgitos, no otorgar como resultado un valor exacto. Cuanto mayor sea el nmero de dgitos utilizados, mayor ser la precisin obtenida. Como vimos en el ejemplo 1.4, un nmero puede tener representacin nita con respecto a una base, pero una representacin innita en otra base. La base decimal es la que emplearemos generalmente, pero antiguamente fueron empleadas otras bases, como la base 12 y la base 60. Un computador opera normalmente en el sistema binario. Observe los que pasa cuando un usuario interacta con el computador: Los datos de entrada son enviados al computador por el usuario en el sistema decimal, esa informacin es convertida al sistema binario por el computador, y, todas las operaciones son realizadas en ese sistema. Los resultados nales sern convertidas para el sistema decimal y, nalmente, sern transmitidos hacia el usuario. Todo este proceso es una fuente de errores que afectan el resultado nal de los clculos.

2.1.3.

Las operaciones numricas efectuadas

Pero errores no slo ocurren en la impresin de los datos de entrada y su representacin binaria. Errores ocurren tambin en las operaciones numricas efectuadas por un sistema de cmputo (binario). Conversin de Nmeros en los Sistemas Decimal y Binario Cualquier nmero entero en la base , de la forma (aj aj 1 , ..., a2 a1 a0 ) , donde 0 ak 1 y k = 0, ..., j, puede ser escrito en la forma polinomial aj j + aj 1 j 1 + ... + a2 2 + a1 1 + a0 0 Mediante esa representacin, podemos convertir fcilmente un nmero entero representado en el sistema binario para el decimal, e inversamente. Por ejemplo: (10111)2 puede ser representado por 1 24 + 0 23 + 1 22 + 1 21 + 1 20 Reordenando y resaltando la base 10 tenemos (10111)2 = 1 24 + 0 23 + 1 22 + 1 21 + 1 20 = 2(23 + 2) + 3 = 2 101 + 3 100 = (23)10

2.1. FACTORES QUE INTERVIENEN EN LA RESOLUCIN NUMRICA

23

Laboratorio 2.4 En algn lenguaje de programacin, haga un programa tal que, dado un nmero entero binario, retorne su equivalente decimal. E inversamente, dado un entero decimal, otorgue su equivalente binario. Como Convertir un Nmero Fraccionario de Representacin Decimal a Binario? Consideremos ahora la conversin de un nmero fraccionario de base 10 para la base 2. Por ejemplo, r = 0, 125, s = 0, 777..., t = 0, 414213562..., etc. Notemos que r tiene una representacin nita, pero s y t tienen representaciones innitas. En trminos generales, dado un nmero entre 0 y 1 en el sistema decimal, Cmo obtener su representacin binaria? Consideremos el nmero decimal fraccionario r = 0, 125, existen dgitos binarios d1 , d2 , ..., dj , ... tal que (0, d1 d2 ...dj ...)2 ser su representacin binaria en la base 2. As, (0, 125)10 = d1 21 + d2 22 + ... + dj 2j + ... Multiplicando cada trmino de la expresin (2.1) por 2, tendremos 2 0, 125 = d1 + d2 21 + d3 22 + ... + dj 2j +1 + ... Por lo tanto, d1 representa la parte entera de 2 0, 125, que es igual a 0, mientras que d2 21 + d3 22 + ... + dj 2j +1 + ... representa la parte fraccionaria de 2 0, 125, que es 0, 250. Aplicando ahora el mismo procedimiento para 0, 250, tendremos 0, 250 = d2 21 + d3 22 + ... + dj 2j +1 + ... 2 0, 125 = 0, 5 = d2 + d3 21 + d4 22 + ... + dj 2j +2 + ... de donde d2 = 0. Repitiendo el procedimiento para 0, 5 tenemos 0, 5 = d3 21 + d4 22 + ... + dj 2j +2 + ... 2 0, 5 = 1 = d3 + d4 21 + ... + dj 2j +3 + ... de donde d3 = 1. Como la parte fraccionaria de 2 0, 5 es cero, el proceso de conversin termina. En resumen tenemos: d1 = 0, d2 = 0 y d3 = 1. Por lo tanto, el nmero (0, 125)10 tiene representacin nita en la base 2: (0, 125)10 = (0, 001)2 Laboratorio 2.5 Usando los procedimientos anteriores para convertir nmeros fraccionarios decimales, a binarios, haga un programa usando algn lenguaje de programacin y verique que: 1. El nmero (0, 5)10 tiene una representacin binaria nita (0, 1)2 (2.1)

24

CAPTULO 2. NOCIONES BSICAS SOBRE ERRORES 2. El nmero (0, 11)10 tiene una representacin binaria innita (0, 0001110000101000111101)2 3. Verique cuntos dgitos el computador est considerando.

En hecho de que un nmero no tenga representacin nita en el sistema binario, puede ocasionar la ocurrencia de errores aparentemente inexplicables en los clculos efectuados en sistemas computacionales binarios. Un computador que opera en el sistema binaio, necesariamente tendr que almacenar una aproximacin para (0, 11)10 , debido a que slo posee una cantidad ja de posiciones para guardas los dgitos de la mantisa de un nmero. Al ser esta aproximacin usada para realizar los clculos, no puede esperarse un resultado exacto. Laboratorio 2.6 (Precisin de una mquina) La precisin de la mquina se dene como el menor nmero > 0 en aritmtica de punto otante, tal que (1 + ) > 1. Este nmero depende totalmente del sistema de representacin de la mquina: base numrica, total de dgitos en la mantisa, de la forma cmo son realizadas las operaciones y del compilador utilizado. Es importante conocer la precisin de la mquina porque en varios algoritmos se requiere ingresar como dato de entrada un valor positivo, prximo de cero, para que sea usado como criterio de comparacin para la detencin del algoritmo. El siguiente algoritmo determina dicha precisin: Paso 1: A 1, S 2

Paso 2: Mientras S > 1, hacer A A 2 S 1+A

Paso 3: Hacer prec = 2A e imprimir prec . 1. Haga un programa en algn lenguaje de programacin que ejecute el algoritmo anterior. 2. Discuta su signicado prctico. Aritmtica de Punto Flotante Cualquier computador o calculadora representa un nmero en un sistema denominado aritmtica de punto otante. En este sistema, el nmero r ser representado en la forma: (.d1 d2 ...dt ) e donde : La base en que la mquina opera t : El nmero de dgitos en la mantisa, 0 dj 1, j = 1, ..., t, d1 6= 0 e : El exponente en el intervalo [u; u] , el valor de u depende de la mquina con que se est trabajando. En una computadora, slo una pequea cantidad de nmeros son representados exactamente, por lo general, la representacin ser realizada por medio de truncamiento o redondeo.

2.1. FACTORES QUE INTERVIENEN EN LA RESOLUCIN NUMRICA

25

Ejemplo 2.5 Considere una mquina que opera en el sistema = 10, t = 3, e [5; +5] . Los nmeros no nulos representados en este sistema sern de la forma (.d1 d2 d3 ) 10e , 0 dj 9, d1 6= 0, e [5; +5]

El menor nmero, no nulo y en valor absoluto, expresado en esta mquina ser: m = 0, 100 105 Mientras que el mayor nmero es: M = 0, 999 10+5 Ahora, en esta misma mquina. Consideremos el subconjunto de nmeros reales caracterizados por: G = {x R : m |x| M } Puede ocurrir varios casos: x G Por ejemplo, si x = 235, 89 = 0, 23589 10+3 . Este nmero posee 5 dgitos en la mantisa. Debido a que t = 3, este nmero no ser considerado de forma exacta en esta mquina. Si la mquina usa truncamiento, entonces el nmero ser representado como 0, 235 10+3 . Pero si la mquina usa redondeo, entonces el nmero ser representado por 0, 236 10+3 . |x| < m Por ejemplo, si x = 0, 345 107 . Este nmero no puede ser representado en esta mquina porque el exponente e es menor que 5. La mquina en estas condiciones retorna una advertencia de underow . |x| > M Por ejemplo, x = 0, 875 109 . En este caso, el exponente es mayor que 5 y la mquina no lo puede representar, advierte la ocurrencia de overow . Algunos lenguajes de computador permiten que las variables sean declaradas en doble precisin. En este caso, tal variable ser representada en el sistema de aritmtica de punto otante de la mquina, pero con aproximadamente el doble de dgitos disponibles en la mantisa. Debemos resaltar que en estas condiciones, el tiempo de ejecucin y requerimientos de memoria aumentan considerablemente. La adicin en aritmtica de punto otante requiere el alineamiento de los puntos decimales de los dos nmeros. Para eso, la mantisa de menor exponente debe ser desplazada para la derecha. Ejemplo 2.6 Sumar Alineando los puntos decimales, tenemos x = 0, 937 104 y y = 0, 1272 102 y y = 0, 001272 104

x = 0, 937 104 Entonces,

El resultado almacenado despus del truncamiento ser 0, 9382 104 . Mientras que despus del redondeo ser 0, 9383 104 .

x + y = (0, 937 + 0, 001272) 104 = 0, 938272 104

26

CAPTULO 2. NOCIONES BSICAS SOBRE ERRORES

El cero puede representarse con una mantisa nula y cualquier exponente. Por lo general, se utiliza el menor exponente posible de la mquina. Caso contrario, si se usa cualquier exponente para denotar el cero, se puede perder dgitos signicativos, tal como muestra el siguiente ejemplo: Ejemplo 2.7 Supongamos que tenemos una mquina que opera con base 10 y 4 dgitos en la mantisa. Si denotramos al cero por 0, 0000 104 , al sumarlo al nmero y = 0, 3134 102 : 0, 0000 104 + 0, 3134 102 = 0, 0000 104 + 0, 003134 104 = (0, 0000 + 0, 003134) 104 = 0, 003134 104 El resultado despus del truncamiento sera 0, 0031 104 = 0, 3100 102 . Esto signica que fueron perdidos 2 dgitos del valor exacto. Ejemplo 2.8 Represente los siguientes nmeros en un sistema de aritmtica de punto otante (con redondeo y con truncamiento) de 3 dgitos, cuando = 10, m = 104 y M = 10+4 : 3, 14 10, 052 238, 15 2, 71828 0, 000007 718235, 82

Para el primer caso, con truncamiento nos resulta 0, 314 10, mientras que con redondeo 0, 314 10. Para el segundo caso, con truncamiento obtenemos 0, 100 102 y con redondeo 0, 100 102 . Y as sucesivamente: 0, 238 103 0, 271 10 underow overow 0, 238 103 0, 271 10 underow overow

Todo esto nos da una idea de los posibles errores que pueden suceder, ya sea por desconocimiento nuestro o por la limitacin del computador, en el proceso de la resolucin numrica de problemas. Debemos advertir que an es posible realizar un anlisis ms completo del manejo de errores, pero eso lo veremos en otra ocasin.

2.2.

Errores Absolutos y Errores Relativos

Denimos el error absoluto como la diferencia entre el valor exacto de un nmero x y su valor aproximado x : EAx = x x (2.2) Frecuentemente, slo nos interesa la magnitud de este error. As por ejemplo, si x = 12, 60 y x = 12, 81, el error absoluto es EAx = 12, 60 12, 81 = 0, 21. Mientras que la magnitud de este error es |0, 21| = 0, 21. Esta idea puede ser extendida para comparar la proximidad de vectores. Por ejemplo, consideremos los vectores en R3 dados por 6, 1 6, 2 x = 5, 8 y x = 5, 9 11, 3 10, 9

2.2. ERRORES ABSOLUTOS Y ERRORES RELATIVOS Entonces, usando la norma euclidiana, la magnitud del error estara dada ahora por p kx x k = (6, 1 6, 2)2 + (5, 8 5, 9)2 + (11, 3 (10, 9))2 = 0, 4243

27

No obstante, el error absoluto denido en (2.2) quiz no tenga inters prctico en este caso. En general, apenas el valor de x es conocido, lo que hace imposible obtener el error absoluto exacto. Lo que se puede hacer en este caso es obtener una cota superior o una estimativa para el mdulo del error absoluto, tal como muestra el siguiente ejemplo. Ejemplo 2.9 Conocindo que h3, 14; 3, 15i , podemos tomar para x un valor dentro de este intervalo y tendremos que: |EA | = | x | < 0, 01 En este caso diremos que el error absoluto de x con respecto a , en mdulo, es menor que 0, 01. Ms an, diremos que el nmero x est representado con precisin menor que 0, 01. Ejemplo 2.10 En la prctica, a veces no es recomendable controlar alguno procesos basados en el error absoluto. Por ejemplo, si usted gana un premio de S/. 10 000 000, y cuando va a recogerlo le dicen que slo tiene S/. 9 999 990, entonces a usted puede que no le importe la diferencia, pues apenas hay un error absoluto de S/. 10. Pero qu pasa si usted gan S/. 20 de premio, y cuando usted va a recogerlo le dicen que apenas tiene S/. 10, observe que el error absoluto sigue siendo S/. 10, probablemente no le agrade nada esta ltima situacin, pues se trata de la mitad del premio. Para evitar situaciones como la anterior, en la prctica es mejor utilizar otro criterio para medir el error, ste es conocido como error relativo. El error relativo es denido como el error absoluto dividido por el valor exacto, es decir: ERx = EAx xx = x x

Frecuentemente, tambin se suele trabajar con el mdulo de este valor. Observe que el error relativo respecto al primer premio de S/. 10 000 000 es ER10 000 000 = 10 0, 000001 10 000 000

Mientras que el error relativo respecto al segundo premio es ER20 = 10 0, 5 20

Con esto, digamos que en este caso se midi el error con ms justicia. Ejercicio 2.1 Convierta los siguientes nmeros decimales para su forma binaria: 26, 1278 y 0, 1217 Ejercicio 2.2 Convierta los siguientes nmeros binarios para su forma decimal: (101101)2 , (0, 111111101)2 y (0, 1101)2 .

28

CAPTULO 2. NOCIONES BSICAS SOBRE ERRORES

Laboratorio 2.7 El siguiente algoritmo calcula de una forma aproximada la raz n-sima de un nmero no negativo a. Ingresa: a, n y > 0, donde es la precisin deseada ( = 109 ) x=a Mientras |xn a| > (controlando la magnitud del error absoluto) x = x (xn a)/(nxn1 ), x 6= 0 Retorna x (una aproximacin para n a) En algn lenguaje de programacin, haga un programa para ejecutar este algoritmo. Modique el programa para que retorne tambin el nmero de pasos (iteraciones). Cmo utilizara el error relativo para controlar el algoritmo? Laboratorio 2.8 (Clculo de ex ) En algn lenguaje de programacin, haga un programa para calcular ex mediante la serie de Taylor con n trminos. El valor de x y el nmero de trminos de la serie, n, deben ser dados en la entrada de su programa. Para valores negativos de x, el programa debe calcular ex de dos formas: En una de ellas el valor de x es usado directamente en la serie de Taylor y, en la otra forma, el valor usado en la serie ser y = x, y en seguida, se calcula el valor de ex por medio de e1 x. 1. Experimente su programa con varios valores de x (x prximo de cero y distante de cero) y, para cada valor de x, experimente el clculo de la serie con varios valores de n. Analice los resultados. 2. (Dicultades con el clculo del factorial) El clculo de k! necesario en la serie de Taylor puede ser hecho de modo a evitar la ocurrencia de overow. Para esto es necesario k analizar cuidadosamente el k-simo trmino, x . Intente combinar el clculo del numerak! dor con el del denominador y realizar divisiones intermedias. Estudie una manera de realizar esta operacin de modo que k! no se sobrecargue. 3. Con la modicacin del segundo item, la serie de Taylor puede ser calculada con los trminos que se desee. Cul sera el criterio para detener su programa e interrumpir el clculo de la serie?

Captulo 3 Ceros reales de funciones reales

En esta seccin vamos a resolver la ecuacin representada por f (x) = 0 (3.1)

donde f : [a; b] R 7 R. Resolver tal ecuacin signica encontrar [a; b] de modo que f ( ) = 0.

Algunas de las tcnicas para resolver esta ecuacin, es decir, encontrar una raz de la ecuacin (3.1) o simplemente un cero de f , requiere de un procedimiento que comprende esencialmente de dos fases. Fases para encontrar una raiz de la ecuacin f (x) = 0 1. Aislamiento de las raices : (Intervalo que contenga una raz) 2. Renamiento : (Aproximacin a la raz deseada) Cada uno de estas fases sern explicados a continuacin.

3.1.

Aislamiento de las Raices

Si es una raz de f , el procedimiento de aislamiento de una raz consiste en obtener un intervalo [a; b] que contenga .

Una primera alternativa sera mediante una observacin grca, probablemente con ayuda de un computador o una calculadora.

Ejemplo 3.1 Para el aislamiento de una raz de f (x) = x 5ex , procedemos a gracar utilizando Matlab. Observamos que en el intervalo [1; 2] se encuentra una raz. La gura 2.1 29

30 muestra esta situacin.

3

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

0 eje x -1 raz

-2

-3

-4 f(x) = sqrt(x) - 5*exp(-x) -5

Figura 2.1: Aislando grcamente una raz otra alternativa para el aislamiento de una raz consiste en analizar el cambio de signo de los valores de la funcin, tal como explica el siguiente ejemplo. Ejemplo 3.2 Considere nuevamente la funcin f (x) = donde se encuentre por lo menos una raz. x 5ex . Determine un intervalo

Solucin: Tenemos que Df = R+ 0 (los reales positivos y el cero), construimos una tabla de valores con el signo de f (x) para determinados valores de x : x 0 1 2 3 f (x) + + Analizando la tabla, vemos que f admite por lo menos una raz en el intervalo [1; 2] . Para ver si esa raz es nica en ese intervalo, podemos analizar el signo de la derivada de f : 1 f 0 (x) = + 5ex > 0, x > 0 2 x Vemos que f es estrictamente creciente en R+ . Por lo tanto, podemos concluir que f admite una nica raz en Df y sta se encuentra en el intervalo [1; 2] . Ejercicio 3.1 Sea f : R 7 R una funcin continua en R. Si sabemos que una raz de f est en [a; b] . Qu podra usted decir sobre si: 1. f (a)f (b) < 0 2. f (a)f (b) > 0 3. f (a)f (b) 0 4. f (a)f (b) 0

3.2. REFINAMIENTO: MTODOS ITERATIVOS

31

Laboratorio 3.1 Usando un computador y un programa adecuado, graque las siguientes funciones y determine aquellos intervalos que incluyan alguna raz: 1. f (x) = 1000 sen(x3 + 1)/ log(5 + x2 ), x [2; 1] xR xR

2. g (x) = ln(x2 + 1) 200(x + 10)3 + 9x2 + 5,

3. h(x) = 0, 00037x11 (x )2 + x2 + 5x 100,

3.2.

Renamiento: Mtodos Iterativos

Una vez que tenemos un intervalo que contenga la raz, el siguiente paso es construir un mecanismo que nos otorgue aproximaciones razonables a la solucin exacta.En esta seccin veremos algunos mtodos numricos clsicos los cuales nos otorgarn aproximaciones a una raz de f . En su mayora, estos mtodos son iterativos, es decir, que inician con una estimativa de la solucin inicial y utilizan sta para encontrar la siguiente, y as por delante, hasta obtener una aproximacin a la solucin. Cuando se use un mtodo iterativo, debemos considerar un criterio para detener el algoritmo respectivo. En los mtodos que buscan una raz, stos se repetirn hasta que xk sea prxima a la raiz exacta con presicin > 0,esto ocurrir si: 1. xk < , o 2. f (xk ) < [a; b] y ba<

Pero, Cmo efectuar el primer item si no se conoce ?. Una forma es reducir el intervalo que contiene a la raz en cada iteracin. Al conseguirse un intervalo [a; b] tal que

entonces x [a; b] , |x | < . En estas condiciones, cualquier x [a; b] podra ser la aproximacin requerida.

Figura 2.2: f ( ) = 0

32

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

El orden de grandeza de los nmeros con que trabajamos puede darnos poca informacin, tal como mostraba el ejemplo 1.10, es aconsejable en estos casos utilizar el error relativo. Por ejemplo, podemos considerar xk prximo de una raz, si f (xk ) < L

Otro aspecto que debemos tener en cuenta es el mximo nmero de iteraciones permitidas por el algoritmo. Esto ayuda a evitar que el programa en computador trabaje indenidamente, sobre todo en el caso cuando el algoritmo no converge. Antes de todo, debemos hacer una aclaracin con respecto a mtodo y algoritmo. Entenderemos por mtodo a un procedimiento con las justicaciones matemticas necesarias para resolver un determinado problema. Mientras que por algoritmo, entenderemos como un resumen del mtodo, una especie de receta. Existen varios mtodos numricos para obtener un cero real de una funcin real, algunos simplemente requieren que la funcin sea continua, mientras que otros requieren que la funcin sea diferenciable. En lo que resta de este captulo analizaremos cada uno de los mtodos ms populares que existen hoy en da. En nuestro caso, analizar comprender la construccin del mtodo, estudiar las propiedades de convergencia y la rapidez del mismo. Para estudiar la convergencia debemos dar las hiptesis para que el mtodo garantice una solucin. Por otro lado, para analizar la rapidez del mtodo, es necesario tener en consideracin dos criterios: El nmero de iteraciones: Dada la precisin deseada , determinar el nmero de iteraciones, k, para que el algoritmo respectivo se detenga. La rapidez: Una vez garantizada la convergencia.Determinar cul es la tasa o rapidez (velocidad) de convergencia con que trabaja el algoritmo. Lo ms deseable es obtener el nmero de iteraciones que un algoritmo requiere para resolver el problema, al menos una cantidad aproximada de ste, pero no siempre es posible tal hazaa. A veces es posible obtener slo la tasa de convergencia del algoritmo, esto tambin dar informacin sobre el desempeo del mismo, lo cual permitir realizar comparaciones para decidir por el algoritmo ms eciente para un determinado problema. Los detalles relacionados a estos conceptos sern explicados a medida que vayamos avanzando.

donde L = f (xk ) , para algn x escogido en una vecindad de la raz .

3.2.1.

Mtodo de Biseccin

Sea f una funcin contnua en el intervalo [a; b] tal que f (a)f (b) < 0. El objetivo de este mtodo es reducir la amplitud de este intervalo que contiene la raz hasta alcanzarse una

3.2. REFINAMIENTO: MTODOS ITERATIVOS precisin requerida, (b a) < , usando para eso una sucesiva divisin de [a; b] a la mitad.

33

Figura 2.3: Mtodo de Biseccin Algoritmo 3.1 (Biseccin) Dados a y b tales que f (a)f (b) < 0. Sea > 0 la precisin deseada 1. Si (b a) < . Eligir x [a; b] y terminar el algoritmo. Caso contrario, ir al paso 2. 2. Hacer k = 1 e ir al paso 3. 3. Hacer c =
a+b 2

e ir al paso 4.

4. Si f (a)f (c) > 0, hacer a = c e ir al paso 5. Caso contrario, hacer b = c e ir al paso 5. 5. Si b a < , elegir x [a; b] y nalizar el algoritmo. Caso contrario, hacer k = k + 1 y volver al paso 3. Estudio de la Convergencia del Mtodo de Biseccin Bajo las hiptesis establecidas, es claro que el mtodo de la biseccin construir una sucesin {xk } que converge a una raz. Para probar esto analticamente procedemos del siguiente modo. Supongamos que [a0 ; b0 ] sea el intervalo inicial y la nica raz de f en el interior de ese intervalo. El mtodo de la biseccin genera tres sucesiones: La sucesin {ak }kN : No decreciente y acotada superiormente por b0 . Luego, existe r R tal que limk ak = r. La sucesin {bk }kN : No creciente y limitada inferiormente por a0 . Luego, existe s R tal que limk bk = s. La sucesin {ck }kN : Genera por ck =
ak +bk , 2

donde ak < ck < bk , para todo k = 1, 2, ... .

Observe que el tamao de cada intervalo es la mitad del intervalo anterior. As, para k = 1, 2, ... a0 b1 a1 = b0 2 a1 a0 b2 a2 = b1 = b02 2 2 b2 a2 b0 a0 b3 a3 = 2 = 22 . . .

34

CAPTULO 3. CEROS REALES DE FUNCIONES REALES bk ak = Entonces,

k

b0 a0 2k (b0 a0 ) =0 2k

(3.2)

l m (bk ak ) = l m

Como {ak }kN y {bk }kN son convergentes:

k

m ak = 0 s r = 0 l m bk l
k

Por lo tanto, s = r, Sea = s = r el lmite de las dos sucesiones. Dado que para todo k = 1, 2, ... el punto ck hak ; bk i , entonces
k

l m ck =

Resta probar apenas que es un cero de la funcin f, o sea, f () = 0. Observe que en cada iteracin k = 1, 2, ..tenemos que f (ak )f (bk ) < 0. Entonces, dado que hemos asumido f continua en [a; b] 0
k

m f (ak ) l m f (bk ) l m f (ak )f (bk ) = l

k k k 2 k

m bk ) = f (r)f (s) = f ()f () = f ( l m ak )f ( l = (f ())

de donde concluimos que f () = 0. Por tanto, limk ck = es un cero de f, como habamos asumido que en el interior haba una nica raz, tenemos que = . Estimacin del Nmero de Iteraciones Dada una precisin > 0 y un intervalo inicial [a; b] , es posible saber, a priori, cuntas iteraciones sern efectuadas por el mtodo de la biseccin (algoritmo 2.1) hasta que se cumpla b a < . Vimos en (3.2) que para k = 1, 2, ... bk ak = bk1 ak1 b0 a0 = 2 2k (3.3)

Observe que el algoritmo (2.1) se detendr cuando bk ak < , segn la ecuacin (3.3) esto equivale a encontrar un valor de k de modo que b0 a0 < 2k Esto a su vez equivale a decir que 2k > Lo cual implica que k > log2 (b0 a0 ) log2 () (3.4) b0 a0

3.2. REFINAMIENTO: MTODOS ITERATIVOS Por lo tanto, en el algoritmo 2.1, si el nmero de iteraciones k es por lo menos [log2 (b0 a0 ) log2 ()] + 1

35

el intervalo [ak ; bk ] conteniendo a la raz verica bk ak < , en estas condiciones cualquier x [ak ; bk ] debera ser la aproximacin a la raz buscada , pues |x | bk ak < . Podemos concluir sobre este algoritmo que el nmero de iteraciones depender de la longitud del intervalo [a; b] y de la precisin deseada . Como muestra la ecuacin (3.4), ese nmero de iteraciones no debera ser grande, debido a la presencia del logaritmo. Ejemplo 3.3 Usando el algoritmo 2,1, se desea encontrar una aproximacin a un cero de la funcin denida por f (x) = x log10 x 1 la cual est en el intervalo [2; 3] y con precisin = 102 . Cuntas iteraciones debera efectuar el algoritmo? Solucin: Segn (3.3) vemos que k log2 (3 2) log2 (102 ) + 1 = log2 1 log2 (102 ) + 1 = log2 (102 ) + 1 = 7 Luego, el algoritmo debera detenerse con k = 7 iteraciones.

Conclusin 2.1 (sobre el mtodo de biseccin) Si f : R 7 R es continua en el intervalo [a; b] y f (a)f (b) < 0 : El mtodo de biseccin genera una sucesin que converge a la raz de f (x) = 0, esto se consigue mediante las reducciones sucesivas del intervalo de bsqueda hasta una precisin deseada. Cada iteracin no requiera clculo complicados. No se requieren derivadas. Las hiptesis no son rigurosas. Converge razonablemente rpido (comparado a otros mtodos se le considera lento). Laboratorio 3.2 En algn lenguaje de programacin de su preferencia, implemente el algoritmo de la biseccin, y resuelva las siguientes ecuaciones: 1. ln(x2 + 1) = 200(x + 10)3 9x2 5, en R 2. 1000 sen(x3 + 1)/ log(5 + x2 ) = 0, en [2; 1] 3. 0, 00037x11 (x )2 + x2 = 5x + 100, en R

36

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

3.2.2.

Mtodo de la Posicin Falsa

Sea f : R 7 R continua en [a; b] tal que f (a)f (b) < 0. Suponga que el intervalo [a; b] contiene una nica raz de la ecuacin f (x) = 0. Podemos esperar conseguir una raz aproximada usando las informaciones sobre los valores de f disponibles a cada iteracin.

Por ejemplo, en la gura 2.4 se aprecia que al ser |f (a)| pequeo en comparacin a |f (b)|, podemos sospechar que la raz se encuentra ms cercana al punto a que al punto b. Luego, en cada iteracin, en vez de tomar ck como el punto medio, como lo haca el mtodo de biseccin, podemos tomarlo de la siguiente manera: a |f (b)| + b |f (a)| |f (b)| + |f (a)|

ck =

que en realidad es una media aritmtica ponderada entre a y b, con pesos |f (b)| y |f (a)| . Despus de unos clculos, tenemos: af (b) bf (a) f (b) f (a)

ck =

Lo que resta del mtodo de la posicin falsa es anlogo al mtodo de biseccin, la parte donde no se encuentra la raz debera ser desechada y el intervalo debera ser reducido hasta una precisin deseada.

Figura 2.4: xito en el mtodo de la posicin falsa: |f (a)| pequeo y raz cercana de a Lamentablemente, las cosas no son como parecen, pues el valor de |f (a)| puede ser pequeo y sin embargo, la raz puede estar muy lejos de a, tal como se aprecia en la gura 2.5. Esto muestra que la sospecha puede estar totalmente errada, lo que retrasara la convergencia del mtodo, de ah el nombre.

3.2. REFINAMIENTO: MTODOS ITERATIVOS

37

Figura 2.5: Fracaso en el mtodo de la posicin falsa: |f (a)| pequeo y raz alejada de a. Ejercicio 3.2 Haga un algoritmo que resuma el mtodo de la posicin falsa. Laboratorio 3.3 En algn lenguaje de programacin de su preferencia, implemente el algoritmo de la posicin falsa, pruebe con varios ejemplares y compare el nmero de iteraciones con el mtodo de biseccin. Ejercicio 3.3 Investigue sobre la convergencia y el nmero de iteraciones requeridas por el mtodo de la posicin falsa.

3.2.3.

Mtodo del Punto jo (MPF)

Sea f : R 7 R una funcin. Se dice que es un punto jo de f, si f () = . Este concepto es el mismo para una funcin vectorial de variable vectorial. El mtodo del Punto Fijo consiste en lo siguiente: 1. Transforma la ecuacin f (x) = 0 en una ecuacin equivalente: x = (x) 2. Dado un punto inicial x0 R, genera una sucesin {xk } de aproximaciones hacia , la raz buscada, mediante la relacin xk+1 = (xk ) Observe que la funcin debe ser una funcin que cumpla: f ( ) = 0 si, y slo si, ( ) = . De este modo, resolver el problema de encontrar una raz de una ecuacin se convierte en un otro problema de hallar un punto jo de una funcin. Aunque a simple vista no parezca, ms adelante veremos que esta idea trae ciertas ventajas. La funcin con esa caracterstica se denomina funcin iteracin asociada a la ecuacin f (x) = 0. Naturalmente, pueden existir muchas funciones iteracin asociadas a una sola ecuacin.

38

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

Ejemplo 3.4 Dada la ecuacin x2 + 2x 10 = 0. Algunas candidatas a funcin iteracin son las siguientes: 1. (x) = 5 2. (x) = 3. (x) =
10 , x+2 10 x x2 2

x 6= 2 x 6= 0

2,

Denicin 3.1 (Forma general de una funcin iteracin) Una funcin iteracin asociada a la ecuacin f (x) = 0, est dada de una forma general por (x) = x + A(x)f (x) con la condicin que en , punto jo de , se tenga A( ) 6= 0. Teorema 3.4 Si es una funcin iteracin de la ecuacin f (x) = 0, entonces f ( ) = 0 si, y slo si, ( ) = . Prueba. (.) Sea tal que f ( ) = 0. Como ( ) = . + A( )f ( ), claramente tenemos ( ) = . () Si ( ) = , entonces . + A( )f ( ) = implica que A( )f ( ) = 0, esto a su vez implica que f ( ) = 0, pues A( ) 6= 0 por denicin. Estudio de la Convergencia del Mtodo del punto Fijo Dependiendo de la eleccin de la funcin iteracin , el mtodo del punto jo puede o no convergir a la solucin de la ecuacin f (x) = 0. El siguiente teorema establece las condiciones sucientes para que esta convergencia suceda. Teorema 3.5 Sea una raz de la ecuacin f (x) = 0, aislada en un intervalo abierto I centrado en . Sea una funcin iteracin asociada a esta ecuacin. Si 1. y 0 son funciones contnuas en I. 2. |0 (x)| M < 1, para todo x I. 3. x0 I. Entonces, la sucesin {xk } generada por la regla xk+1 = (xk ), hacia . Prueba. La prueba consta de dos partes: 1. Si x0 I, entonces xk I, para todo k = 0, 1, 2, ... 2. l mk xk = k = 0, 1, 2, ...converge

3.2. REFINAMIENTO: MTODOS ITERATIVOS

39

Primero, como es una raz exacta de la ecuacin f (x) = 0, se tiene que f ( ) = 0 = ( ). Adems, para cualquier k = 0, 1, 2, ...se tiene xk+1 = (xk ), desde aqu xk+1 = (xk ) ( ) Como es continua y diferenciable en I, por el teorema del valor medio, existe ck entre xk y , tal que xk+1 = (xk ) ( ) = 0 (ck )(xk ) k = 0, 1, 2, ... |xk+1 | = |0 (ck )| |(xk )| M |xk | < |xk | k = 0, 1, 2, .. (3.5)

Luego,

Es decir, la distancia entre xk+1 y es menor que la distancia entre xk y , como I est centrado en , vemos que si xk I. Luego, si x0 I, claramente xk I para todo k = 1, 2, .. Para probar la segunda parte, l mk xk = , desde (3.5) vemos que: |x1 | M |x0 | |x2 | M |x1 | M 2 |x0 | . . . |xk | M k |x0 | de donde l mk |xk | = 0, pues 0 M < 1. Por lo tanto, l mk xk = . Ejemplo 3.5 Sea la funcin f (x) = x2 +2x 10 cuya raz es 2, 3166. Dadas las funciones iteracin x2 1 (x) = 5 2 y 10 2 (x) = , x 6= 2 x+2 Observe que |01 (x)| = |x| = |x| < 1 x h1, 1i

Luego, no existe un intervalo I centrado en tal que |01 (x)| < 1 para todo x I. El teorema 2.2 no arma nada con respecto de la convergencia de la sucesin {xk } generada por xk+1 = 1 (xk ), pues 1 no cumple la hiptesis, el mtodo del punto jo puede convergir o no cuando utilice 1 como funcin iteracin. Por otro lado, si usamos la funcin iteracin 2 , la situacin es diferente. Observe que D E D E 10 0 < 1 x ; 10 2 10 2; + |2 (x)| = (x + 2)2 h; 5, 1622i h1, 1622; +i

Luego, existe un intervalo I centrada en , tal que |02 (x)| < 1 para todo x I . En este caso, el teorema asegura la convergencia del MPF tomando x0 I. Ejercicio 3.4 Analice el caso para la funcin f (x) = x2 + 2x 10, cuya raz es 2, 3166, cuando se usa como funcin iteracin: 3 (x) = 10 2, x x 6= 0

Adems, analice el otro caso, cuando consideramos como raz a 4, 1622.

40

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

Algoritmo 3.2 (Punto Fijo) Considere la ecuacin f (x) = 0 y la ecuacin equivalente x = (x). Supongamos que ya es conocida explcitamente y las hiptesis sucientes del teorema 2.2 son satisfechas. Dada una precisin deseada > 0, el algoritmo se detendr cuando |f (xk )| < . La aproximacin al punto jo (la raz buscada) ser x . 1. Datos iniciales: a) x0 , la aproximacin inicial b) 1 y 2 las precisiones deseadas 2. Si |f (x0 )| < 1 , hacer x = x0 , nalizar el algoritmo 3. k = 1 4. x1 = (x0 ) 5. Si |f (x1 )| < 1 o si |x1 x0 | < 2 , hacer x = x1 , nalizar el algoritmo 6. x0 = x1 7. k = k + 1. Volver a paso 4. Laboratorio 3.6 En algn lenguaje de programacin de su preferencia, implemente el MPF y ejecute el programa sobre el problema del ejemplo 2.5. Observe que el mtodo diverge cuando usamos 1 , pruebe con varios puntos iniciales e interprete los resultados Tasa de Convergencia del Mtodo de Punto Fijo Cuando vimos el mtodo de biseccin, notamos que era posible obtener una cota inferior para el nmero de iteraciones a ser realizadas por el algoritmo. Pero eso no siempre es posible, as, necesitamos de algunos parmetros que nos indiquen con qu rapidez la sucesin generada por un algoritmo converge a la solucin deseada. Eso nos permitir calicar a un algoritmo como lento o rpido. Denicin 3.2 (Tasa o rapidez de convergencia) Sea la sucesin r(k) generada por al. Asumamos que r(k) 6= r para todo k = 0, 1, 2, .., gn algoritmo, de modo que l mk r(k) r la tasa de convergencia del algoritmo es el supremo P de los nmeros no negativos p satisfaciendo: (k+1) r r l m p = < k kr (k) r k

Si P = 1 y el radio de convergencia < 1, la sucesin se dice que tiene una tasa de convergencia lineal (por lo menos lineal). Si P 1 y = 0, la sucesin tiene tasa de convergencia sper lineal. Si P = 2 y < , entonces diremos que la sucesin tiene una tasa de convergencia cuadrtica Una manera natural de ver esta situacin es la siguiente: supongamos que r(k) es una sucesin generada por un algoritmo la cual converge a la solucin r , donde el algoritmo presenta una tasa de convergencia lineal, entonces: (k+1) r r r(k) r , k k0 | | {z } {z }
ek+1 ek

3.2. REFINAMIENTO: MTODOS ITERATIVOS

41

nos dice que el error ek+1 cometido en la iteracin k + 1 es menor (linealmente) que el error ek en la iteracin k, cuando k0 es grande. En el caso de una tasa de convergencia cuadrtica, es claro que el error cometido en la iteracin k + 1 es aproximadamente el cuadrado del error cometido en la iteracin anterior. Esto indica que para valores grandes de k0 y en las proximidades de r , el error ek+1 disminuye considerablemente con respecto a ek : 2 (k+1) r r r(k) r , | {z } | {z }
ek+1 e2 k

k k0

As, un algoritmo con tasa de convergencia cuadrtica convergir con mayor rapidez hacia un punto de acumulacin, que uno que posee tasa de convergencia sper lineal. Por otro lado, un algoritmo con tasa de convergencia sper lineal ser ms rpido que uno con tasa de convergencia lineal. Ms adelante veremos que el Mtodo Secante tiene una tasa de convergencia P = 1, 618... Proposicin 3.1 (Tasa de convergencia de MPF) Asumamos que la sucesin generada por el Mtodo del Punto Fijo converge a , la raz de f , entonces la tasa de convergencia es por lo menos lineal Prueba. Desde (3.5), tenemos |xk+1 | = |0 (ck )| |xk | , Luego |xk+1 | = |0 (ck )| , |xk | ck entre xk y , k = 0, 1, 2, ... ck entre xk y , k = 0, 1, 2, ...

Tomando lmites, por la continuidad de y 0 , tenemos |xk+1 | m ck = |0 ( )| = M < 1 = l m |0 (ck )| = 0 l k |xk | k k l m As, vemos que el mtodo del punto jo posee una tasa de convergencia lineal, por este motivo mucha veces se le considera lento. Sin embargo, el mtodo de biseccin posee esa misma velocidad, tal como lo arma el siguiente ejercicio. Ejercicio 3.5 Muestre que la tasa de convergencia del mtodo de la biseccin es lineal.

3.2.4.

El Mtodo de Newton-Raphson

Cuando estudiamos el mtodo del punto jo, notamos que: 1. Una de las condiciones para la convergencia es que |0 (x)| M < 1, para todo x I, donde I es un intervalo centrado en la raz. 2. La convergencia del mtodo ser ms rpida cuando menor sea |0 ( )| .

42

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

La segunda armacin se debe al siguiente hecho: cuando analizamos la tasa de convergencia del MPF, vimos que |xk+1 | l m = |0 ( )| < 1 k |xk |

Entonces, acelerar la convergencia del MPF se conseguira escogiendo una funcin iteracin de modo que 0 ( ) = 0, pues en este caso la tasa de convergencia sera por lo menos sper-lineal.

Hacia este objetivo, dada la ecuacin f (x) = 0 cuya raz es . Consideramos la forma general de la funcin iteracin (x) = x + A(x)f (x), donde A( ) 6= 0, con la nueva condicin 0 ( ) = 0. As, (x) = x + A(x)f (x) = 0 (x) = x + A0 (x)f (x) + A(x)f 0 (x) Luego, 0 ( ) = x + A0 ( )f ( ) + A( )f 0 ( ) = 0 ( ) = 1 + A( )f 0 ( ) de donde 0 ( ) = 0 1 + A( )f 0 ( ) = 0 A( ) = Esto nos motiva a denir A( ) = 1 f 0 ( ) 1 f 0 ( ) , f 0 ( ) 6= 0

Entonces, dada la ecuacin f (x) = 0, la funcin iteracin deseada ser de la forma: (x) = x Observe que 0 (x) = 1 (f 0 (x))2 f (x)f 00 (x) f (x)f 00 (x) 0 = 2 2 . f (x) 6= 0 0 0 (f (x)) (f (x)) f (x) . f 0 (x) 6= 0 f 0 (x) (3.6)

y como f ( ) = 0, se tiene que 0 ( ) = 0. Por lo tanto, usando la funcin iteracin denida en (3.6) obtenemos un caso particular del mtodo del punto jo, a ste se le denomina mtodo de Newton. En resumen, dado x0 R, el Mtodo de Newton consiste en construir una sucesin denida por la regla f (xk ) , f 0 (xk ) 6= 0, k = 0, 1, 2, .. (3.7) xk+1 = xk 0 f (xk ) Otro Enfoque del Mtodo de Newton-Raphson Desde el punto de vista geomtrico, el mtodo de Newton puede ser visto como la solucin de un problema difcil, mediante la sucesiva resolucin de problemas fciles.Es decir, dada una aproximacin inicial xk R a la raz buscada, el problema dicil ser hallar una raz de la ecuacin no lineal f (x) = 0, mientras que el problema fcil asociado ser resolver la ecuacin Lk (x) = 0, donde L es una funcin lineal afn que es parecida, al menos localmente, a la funcin no lineal f en torno al punto xk.

3.2. REFINAMIENTO: MTODOS ITERATIVOS

43

As, sea el problema (difcil) que consiste en hallar una raz de f (x) = 0 y x0 R una aproximacin inicial. Por el teorema de Taylor, existe > 0 tal que f (x) L0 (x) = f (x0 ) + f 0 (x0 )(x x0 ) para todo x hx0 , x0 + i. Luego, denotando por x1 la solucin de la ecuacin lineal L0 (x) = 0 y asumiendo que f 0 (xk ) 6= 0, entonces L0 (x) = 0 si, y slo si, f (x0 ) + f 0 (x0 )(x x0 ) = 0 de donde x1 = x0 f (x0 ) f 0 (x0 ) (3.8)

Esperando que x1 sea una mejor aproximacin que x0 a la solucin de f (x) = 0. Este procedimiento puede ser repetido iterativamente, crendose una sucesin {xk } k=0 , donde xk+1 = xk f (xk ) , f 0 (xk ) f 0 (xk ) 6= 0, k = 0, 1, 2, ..

Bajo algunas hiptesis se consigue que l mk xk = .

Figura 2.6: Iteraciones Newton Convergencia del Mtodo de Newton-Raphson A continuacin damos las condiciones bajo las cuales se asegura la convergencia del mtodo de Newton. Teorema 3.7 Sean f ,f 0 y f 00 contnuas en un intervalo abierto I que contiene en su interior I conteniendo la raz de f (x) = 0, donde f 0 ( ) 6= 0, entonces existe un intervalo abierto I la sucesin {xk } generada por la raz , tal que si x0 I, xk+1 = xk convergir para . f (xk ) , f 0 (xk ) f 0 (xk ) 6= 0, k = 0, 1, 2, ..

44

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

Prueba. Observe que el mtodo de newton es en realidad un MPF con la funcin iteracin f (x) (x) = x f 0 (x) .As, para probar la convergencia debemos probar que existe un I I, centrado en , tal que: 1. y 0 son contnuas en I 2. |0 (x)| M < 1, para todo x I
f (x) x)f (x) 0 (x) = f((f y por hiptesis f 0 ( ) 6= 0. Como f 0 es Vemos que (x) = x f 0 (x) , 0 (x))2 contnua en I,es posible obtener un intervalo abierto I1 I, I1 , tal que f 0 (x) 6= 0 para todo x I1 .
00

As, en el intervalo I1 se tiene que f ,f 0 y f 00 son continuas y f 0 (x) 6= 0, x I1 . Por lo tanto, y 0 son contnuas en I1 .
x)f (x) 0 0 Como 0 (x) = f((f 0 (x))2 , es contnua en I1 y ( ) = 0, entonces es posible obtener otro intervalo abierto I2 I1 , centrado en , tal que |0 (x)| < 1 para todo x I2 .
00

I2 centrado en donde y 0 Por lo tanto, hemos encontrado un intervalo abierto I 0 As, si x0 I, la sucesin {xk } generada por la regla de son contnuas y | (x)| < 1, x I. 0 correspondencia xk+1 = xk f (xk )/f (xk ) converge hacia la raz de f (x) = 0. En pocas palabras, lo que dice el teorema anterior es que el mtodo de Newton converge slo si el punto inicial es tomado lo sucientemente prximo de la solucin , esta propiedad se conoce como convergencia local, la cual es una desventaja. Ms adelante veremos que, en cierta forma, ese defecto se compensa con la rapidez con que converge el mtodo. Algoritmo 3.3 (Newton-Raphson) Dada la ecuacin f (x) = 0. Suponga que las hiptesis de suciencia dadas en el teorema 2.3 son satisfechas. El algoritmo otorgar una aproximacin x a la raz Paso inicial: Sea x0 una aproximacin inicial de y 1 , 2 > 0 las precisiones deseadas = x0 y nalizar. Caso contrario hacer: Paso principal: Si |f (x0 )| < 1 , hacer x 1. k = 1 2. x1 = x0 f (x0 )/f 0 (x0 ) = x1 y nalizar. Caso contrario, hacer 3. Si |f (x1 )| < 1 o si |x1 x0 | < 2 , hacer x a) x0 = x1 b) k = k + 1. Volver al paso 2. Laboratorio 3.8 En algn lenguaje de programacin de su preferencia, implemente el algoritmo de Newton y experimntelo con diversos ejemplares. Compare sus resultados con los mtodos anteriormente estudiados. Tasa de Convergencia del Mtodo de Newton-Raphson

3.2. REFINAMIENTO: MTODOS ITERATIVOS

45

Dada una raz exacta de la ecuacin f (x) = 0. Sea {xk } la sucesin generada por el mtodo de newton, tal que l mk xk = . Debido a que el mtodo de Newton es un caso particular de MPF, entonces debe tener por lo menos una tasa de convergencia lineal. Nosotros mostraremos que es mucho ms que eso, debido a la condicin 0 ( ) = 0, el mtodo alcanzar una tasa cuadrtica, es decir, en la denicin 2.2, P = 2. Supongamos adems que se satisfacen todas las hiptesis del teorema 2.3. Observe que xk+1 = xk f (xk ) f (xk ) f (xk ) = xk+1 = xk 0 = ek+1 = ek 0 0 f (xk ) f (xk ) f (xk )

donde ek = xk . Usando la serie de Taylor para f en torno al punto xk : f (x) = f (xk ) + f 0 (xk )(x xk ) + 1 f 00 (ck )(x xk )2 , ck est entre x y ck . Para x = , tenemos 2 1 0 = f (xk ) + f 0 (xk )( xk ) + f 00 (ck )( xk )2 2 de donde 1 f (xk ) = f 0 (xk )(xk ) f 00 (ck )(xk )2 2 0 Dividiendo entre f (xk ), obtenemos f (xk ) f 00 (ck )(xk )2 = ( x ) k f 0 (xk ) 2f 0 (xk ) f 00 (ck )e2 k = ek 2f 0 (xk ) Luego, f (xk ) f 00 (ck )e2 k = ek 0 = ek+1 0 2f (xk ) f (xk ) de donde ek+1 f 00 (ck ) = e2 2f 0 (xk ) k (f 0 (x))2 (f 0 (x)f 00 (x) + f (x)f 000 (x)) (f (x)f 00 (x))((f 0 (x))2 )0 (x) = (f 0 (x))4
00

(3.9)

Despus de unos clculos, vemos que

de donde 00 es continua en y 00 ( ) = Llevando (3.9) al lmite, tenemos |ek+1 | = k |ek |2 l m

f 00 ( ) f 0 ( )

|f 00 (ck )| k |2f 0 (ck )| mk ck )| 1 |f 00 ( )| 1 1 |f 00 (l = = |00 ( )| = = 0 0 2 |f (l mk ck )| 2 |f ( )| 2 l m

Por lo tanto, el mtodo de Newton una tasa de convergencia cuadrtica.

46

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

3.2.5.

Mtodo de Secante

Uno de los incovenientes en el mtodo de Newton es el clculo de la derivada de f en cada iteracin, es decir, evaluar f 0 (xk ). Una manera de enfrentar esto es considerar una aproximacin de la derivada: f (xk ) f (xk1 ) f 0 (xk ) , xk xk1 xk xk1 donde xk y xk1 son dos estimativas de la raz exacta . Para este caso, la funcin iteracin Newton queda establecida como (xk ) = xk f (xk )
f (xk )f (xk1 ) xk xk1

= xk

f (xk )(xk xk1 ) f (xk ) f (xk1 )

O si lo preere, dadas dos estimativas iniciales para la raz de la ecuacin f (x) = 0, x0 y x1 , el mtodo secante consiste en construir una sucesin {xk } denida por la regla xk+1 = xk f (xk )(xk xk1 ) , f (xk ) f (xk1 ) k = 1, 2, 3, ... (3.10)

Las condiciones para la convergencia del mtodo secante son prcticamente las mismas que las del mtodo de Newton. Si bien la dicultad del clculo explcito de la derivada fue evitada, lamentablemente el precio que se paga por esto es la disminucin en la tasa de convergencia, que es sper-lineal. Una interpretacin grca puede ser vista en la gura 2.7.

Figura 2.7: MtodoSecante Convergencia del Mtodo Secante Al igual que el Mtodo de Newton, la convergencia de Mtodo Secante est asegurada cuando los puntos iniciales x0 y x1 fueron tomados lo sucientemente prximos de la raz buscada . El teorema a seguir formaliza esta armacin. Teorema 3.9 Sea una raz de f (x) = 0 y suponga que f C 2 en una vecindad en torno de donde f 0 ( ) 6= 0. Si los puntos iniciales x0 y x1 estn sucientemente prximo de , entonces l mk xk = .

3.2. REFINAMIENTO: MTODOS ITERATIVOS Prueba. Desde la ecuacin (3.10), tenemos xk+1 = xk + f (xk )(xk xk1 ) f (xk ) f (xk1 ) f ()f (x
xk

47

= ( xk1 )( xk )

xk1 f (xk )f (xk1 ) xk xk1

k ) f (xk )f (xk1 ) xk xk1

(3.11)

f 00 ( k ) = ( xk1 )( xk ) 2f 0 ( k )

donde k est entre , xk y xk1 , mientras que k est entre xk y xk1 . Como f C 2 , existe una vecindad I = [ , + ], con > 0, tal que f 0 y f 00 son continuas y f 0 (x) 6= 0 para todo x I . Por tanto, existe M > 0 tal que 00 f ( k ) M = m ax xI 2f 0 ( k ) Elijamos los puntos iniciales x0 , x1 I, de modo que M | x0 | < 1 y M | x1 | < 1 Denamos ahora t = m ax {M | x0 | , M | x1 |} , claramente t < 1. As, por (3.11) tenemos f 00 ( k ) M | x2 | = M ( x1 ) ( x0 ) 2f 0 ( k ) = M 2 | x1 | | x0 | t2 t < 1 | x2 | <

t = m ax {| x0 | , | x1 |} < M As, x2 I. Usando este argumento, tenemos por induccin que xk I, para k 2. En general se tiene M | x2 | = M 2 | x1 | | x0 | t2 M | x3 | = M 2 | x2 | | x1 | t3 M | x4 | = M 2 | x3 | | x2 | t5 . . . M | xk+1 | = M 2 | xk | | xk1 | tFk de donde Fk es el k simo trmino de la sucesin de Fibonacci. Por lo tanto l mk M | xk | = 0.

de donde

Tasa de Convergencia del Mtodo Secante A continuacin analizaremos la tasa de convergencia que el Mtodo Secante posee. Debemos recalcar que algunos autores denominan a esta propiedad como el orden de convergencia.

48

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

Teorema 3.10 Sea una raz de f (x) = 0. Asumamos que la sucesin generada por el Mtodo Secante converge hacia . Si f C 2 y f 0 ( ) 6= 0, entonces el Mtodo Secante tiene una tasa de convergencia p = 1+2 5 = 0, 618... y p1 00 xk+1 f ( k ) l m p = k ( xk ) f 0 ( k ) Prueba. Desde (3.11) tenemos xk+1 f 00 ( k ) = ( xk+1 ) ( xk ) f 0 ( k ) (3.12)

donde k est entre , xk y xk1 , mientras que k est entre xk y xk1 . Resolviendo (3.12) obtenemos 1 n n ( xk ) = Ap B r (3.13) K 1 1 r p f 00 ( k ) (( x0 )K )r rp ( x1 )K 1+ 5 1 5 donde K = 2 , p = , r = , , A = y B = , la p f 0 ( k ) 2 2 ( x1 )K (( x0 )K ) solucin de (3.12) dada por (3.13) se puede comprobar por la sustitucin directa. Por lo tanto, xk+1 p1 rn B (r p) p = K ( xk ) y p1 00 xk+1 f ( k ) p1 rn m K l m B (r p) = l m p = l k ( xk ) k k 2f 0 ( k )

xk+1 Como se puede observar, l mk ( existe para p = 1+2 5 1, 61803398, el cual en cierta xk )p forma nos recuerda al nmero ureo, 0, 61803398..., concluimos entonces que la tasa o rapidez (orden) de convergencia del Mtodo Secante es aproximadamente 1, 618.

Como se puede observar, el Mtodo Secante posee una tasa de convergencia un tanto inferior al Mtodo de Newton, que tiene una tasa cuadrtica. No obstante, la tasa de convergencia de 1, 618 es superior a una lineal, y en la prctica esa disminucin en la velocidad de convergencia y su evaluacin en cada iteracin. Laboratorio 3.11 En algn lenguaje de programacin de su preferencia, implemente el mtodo secante y experimntelo con diversos ejemplares. Compare en la prctica sus resultados con los mtodos anteriores estudiados. Ejercicio 3.6 Experimente y compare los mtodos estudiados en este captulo, hallando una raz de la ecuacin f (x) = 0, donde f (x) = x x ln(x) Ejercicio 3.7 Experimente y compare los mtodos estudiados en este captulo, hallando un raz de la ecuacin f (x) = 0, donde f (x) = ex 4x2

3.2. REFINAMIENTO: MTODOS ITERATIVOS Ejercicio 3.8 Sea

49

x2 + x (ln(x) 1) 2 Halle sus puntos crticos (puntos donde f 0 (x) = 0) usando un mtodo iterativo estudiado en este captulo. f (x) = Ejercicio 3.9 Halle un punto donde la funcin alcanza un mnimo sobre el intervalo [10, 20] . f (x) = x2 + x3 (ln(x) 3) + 850

Ejercicio 3.10 (ciclaje en el mtodo de newton) Las iteraciones del Mtodo de Newton entrarn en un ciclo ilimitado si Esto a su vez sucede si f satisface xn+1 a = (xn a)

f (x) = (x a) f 0 (x) La anterior expresin es una ecuacin diferencial ordinaria separable de la forma f (x) 1 = 0 f (x) 2(x a) cuya solucin es p f (x) = sign(x a) |x a| xa

donde la funcin sign es la funcin signo, el. cero de f es claramente x = a. Grcamente f para el caso a = 2. Seguidamente, use un punto inicial x0 , donde x0 6= a, y ejecute al algoritmo de Newton. Ejercicio 3.11 Se considera la funcin F (x) = x5 + 2x. Mediante el Mtodo de Newton, hallar el menor nmero positivo x (con tres decimales) para el cual F (x) = 4. Ejercicio 3.12 Probar, mediante el mtodo de Newton, que la ecuacin se puede utilizar para aproximar 1/a si x0 es una estimacin inicial del recproco de a. Ntese que este mtodo de aproximar recprocos utiliza slo operaciones de suma y multiplicacin. 1 (considerar f (x) = x a). Pruebe para los casos: 1.
1 3

xn+1 = xn (2 axn )

2.

1 11

Ejercicio 3.13 Pruebe que, al ser aplicado el mtodo de newton en la resolucin ax + b = 0, donde a 6= 0, ste requiere slo una iteracin, sin importar qu punto inicial fue tomado. Ejercicio 3.14 En los casos siguientes, aplicar el mtodo de Newton con la estimacin inicial propuesta, explicar por qu falla el mtodo. 1. y = 2x3 6x2 + 6x 1, 2. y = 4x3 12x2 + 12x 3, donde x0 = 1 donde x0 =
3 2

3. y = x3 + 3x2 x + 1, donde x0 = 1 4. y = 3 x 1, donde x0 = 2

50

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

3.2.6.

Comparacin de los Mtodos Iterativos

El esfuerzo computacional para la ejecucin de cada uno de los mtodos depende de varios factores, los ms importantes son: 1. La complejida de los clculos, sobre todo para la derivada. 2. El nmero total de iteraciones 3. Condiciones para la convergencia El mtodo de la biseccin y el mtodo de la posicin falsa exigen pocas condiciones para garantizar la convergencia, el incoveniente est en que el nmero de iteraciones puede ser grande. Observe que su tasa de convergencia es lineal. Los mtodos de punto jo frecuentemente son ms rpidos, pero a cambio exigen mucha hiptesis para la convergencia. El ms rpido es el mtodo de newton, pero requiere el clculo de la derivada y demanda, al igual que los mtodos de punto jo, hiptesis rigurosas para su convergencia. El mtodo secante puede ser prctico cuando el clculo de la derivada es complicado, pero no es tan rpido como el mtodo de Newton. Se puede concluir que la eleccin del mtodo ms eciente depende de la ecuacin que se intente resolver. Cada mtodo tiene sus ventajas y desventajas. Como un comentario adicional, despus de llevar al computador cada uno de estos mtodos y experimentarlos con diversos ejemplares, probablemente el estudiante halle que las diferencias de tiempo de ejecuacin, entre un programa y otro, sea insignicante cuando se aplica a la resolucin de una ecuacin, y ese afn por buscar el mtodo ms rpido parecera no tener sentido. Esa percepcin es equivocada, pues estos mtodos deben verse como subrutinas de otros mtodos iterativos ms sosticados, para otro tipo de problemas, donde la prdida de una fraccin de segundos restara el desempeo del mtodo en su conjunto.

3.3.

Problemas

Problema 3.1 Considera las dos vigas de 30m y 20m cruzndo a una altura de 8m del suelo como se muestra en la gura 2.8. Determine el ancho del pasadizo, H.

Figura 2.8: Problema de las Vigas

3.3. PROBLEMAS

51

Problema 3.2 La concentracin c de una bacteria contaminada en un lago decrece segn la expresin c(t) = 80e2t + 20e0,5t , siendo t el tiempo en horas. Determine el tiempo que se necesita para que el nmero de bacterias se reduzca a 7. Problema 3.3 Una determinada sustancia se desintegra segn la ecuacin A(t) = P e0,0248t , donde P es la cantidad inicial en el tiempo t = 0 y A la cantidad resultante despus de t aos. Si inicialmente se depositan 500 miligramos de dicha sustancia, Cunto tiempo debe transcurrir para que quede el 1 por ciento de sta? Utilizar el Mtodo de Newton. Problema 3.4 Una medicina administrada a un paciente produce concentracin en la sangre dada por c(t) = Atet/3 mg / ml, t horas despus de que se hagan administrado A unidades. La mxima concentracin sin peligro es de 1 mg / ml, y a esta cantidad se le denomina concentracin de seguridad. 1. Qu cantidad debe ser inyectada para alcanzar como mximo esta concentracin de seguridad? Cundo se alcanza este mximo? 2. Una cantidad adicional se debe administrar al paciente cuando la concentracin baja a 0, 25 mg / ml . Determnese con un error menor de 1 minuto cundo debe ponerse esta segunda inyeccin. Problema 3.5 El crecimiento de poblaciones grandes puede modelarse en periodos cortos suponiendo que el crecimiento de la poblacin es una funcin contnua en t mediante una ecuacin diferencial cuya solucin es: N (t) = N0 et + t e 1

donde N (t) es el nmero de individuos en el tiempo t (medido en aos), es la razn de natalidad, N0 es la poblacin inicial y es una razn constante de inmigracin, que se mide en nmeros de inmigrantes al ao. Supngase que una poblacin dada tiene un millon de individuos inicialmente y una inmigracin de 400, 000 individuos al ao. Se observa que al nal del primer ao la poblacin es de 1506000 individuos. Se pide: 1. Determinar la tasa de natalidad. 2. Hacer una previsin de la poblacin al cabo de tres aos.

52

CAPTULO 3. CEROS REALES DE FUNCIONES REALES

Captulo 4 Resolucin de Sistemas Lineales

Sistemas de ecuaciones lineales aparecen en la resolucin de diversos problemas de la vida real, frecuentemente, algunos de esos problemas involucran un gran nmero de variables y ecuaciones, Claramente, el mtodo basado en la sustitucin de variables, empleado para resolver pequeos sistemas de dos o tres variables, resulta absoleto en estos casos, por lo que en necesario conocer mtodos numricos especializados para hacer frente a estas situaciones. Estos mtodos numricos no podran trabajar sin la intervencin de un computador, debido a la cantidad de clculos envueltos en ese procedimiento. Todo esto sin lugar a dudas torna importante el estudio de sistema de ecuaciones lineales desde el punto de vista numrico. Pero la resolucin de sistemas de ecuaciones lineales no slo tiene importancia prpia, sino que constituye tambin una herramienta indispensable en la resolucin de sistemas de ecuaciones no lineales. Ms an, es la base para la resolucin de muchos otros problemas que surgen en diferentes reas. Considere el siguiente sistemas de ecuaciones lineales compuesto de dos ecuaciones y dos variables: a1,1 x1 + a1,2 x2 = b1 a2,1 x1 + a2,2 x2 = b2 En sistemas simples como ste, podemos notar lo siguiente 1. Solucin nica

3x1 + 2x2 = 10 x1 x2 = 5

x 1 x 2

4 1

2. Soluciones innitas 3x1 3x2 = 15 x1 x2 = 5 3. Ninguna solucin

x 1 x 2

5+t 1

tR

2x1 2x2 = 8 x1 x2 = 5

Grcamente, cada ecuacin representa a una recta, un punto que satisface ambas ecua53

54

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES ciones al mismo tiempo deber estar en la interseccin de ambas rectas (gura 3.1)

Figura 3.1: Representacin grca de sistemas de ecuaciones lineales Por otro lado, cuando tenemos 3 variables, cada ecuacin representara un plano, un punto que satisface tres ecuaciones de un sistema simultneamente estar en la interseccin de tres planos. Para el caso n-dimensional, la situacin es anloga, pero el conjunto denido por cada ecuacin se denomina n-hiperplano .

4.1.

Aspectos Tericos
ecuaciones lineales y n variables es la formulacin + + . ... . . + + + . . . a1,n xn a2,n xn . . . = = . . . b1 b2 . . .

Para un caso general, un sistema de m matemtica siguiente: a1,1 x1 + a1,2 x2 a2,1 x1 + a2,2 x2 . . . . . . . . . a x + a x m,1 1 m,2 2 donde A Rmn , b Rm y x Rn es a1,1 a1,2 a2,1 a2,2 A= . . ... . . . . am,1 am,2

(4.1)

+ am,n xn

= bm (4.2)

Usando la notacin matricial, el mismo sistema puede ser visto por Ax = b el vector de las variables o incgnitas, denidas por: b1 x1 a1,n b2 x2 a2,n , b = y x = . . . . . . . . . am,n bm xn

Una solucin para el sistema (4.1) es un vector x = [ x1 ...x n ]t Rn el cual satisface simultneamente cada una de las m ecuaciones que conforman el sistema: En un sistema lineal puede suceder lo siguiente: 1. El sistema lineal tiene nica solucin. 2. El sistema lineal posee innitas soluciones. 3. El sistema lineal no posee soluciones Cuando n = m y A es inversible en (4.2), entonces la solucin es nica, en consecuencia, dicha solucin puede ser obtenida haciendo x = A1 b. No obstante, se debe advertir que este procedimiento tiene slo un valor terico, pues desde el punto de vista computacional es considerado costoso por el excesivo nmero de operaciones involucradas para calcular la inversa de A. En contraste, existen mtodos ms apropiados que no requieren el clculo explcito de A1 , como los que veremos ms adelante.

4.1. ASPECTOS TERICOS

55

4.1.1.

Conjunto Imagen de una Matriz A

Dada una matriz A Rmn , denimos la imagen de A, denotado por Im(A), como el conjunto: Im(A) = {y Rm : x Rn , y = Ax}

Algunas veces a este conjunto se le denomina espacio imagen de A, debido a lo siguiente: Ejercicio 4.1 Muestre que Im(A) es un subespacio vectorial de Rm . Desde el punto de vista de las columnas de A Rmn , note que
n X Ax = aj xj j =1

donde aj representa la j -columna de A y xj es la j -componente del vector x Rn . Luego, resolver el sistema Ax = b signica obtener escalares x1 , ..., xn que permitan escribir b Rm como una combinacin lineal de las columnas de A, es decir: b = a1 xl + ... + an xn Por esta razn, al conjunto Im(A) tambin se le conoce como espacio columna de A. En consecuencia, la dimensin de la imagen de A, denotada por dim(Im(A)), est determinada por el mximo nmero de columnas linealmente independientes de A.

4.1.2.

Rango de una Matriz A

El rango de la matriz A Rmn , denotado por rango(A), es denido por rango(A) = dim(Im(A)) Una manera prctica de encontrar el rango(A) consiste en determinar el mayor nmero de columnas linealmente independientes de A. En lgebra Lineal se conoce que rango(A) =rango(At ), por consiguiente, para encontrar el rango(A) basta determinar tambin el mayor nmero de las linealmente independientes de A, esto ltimo suele ser de gran ayuda en muchos casos. Con respecto a la resolucin del sistema de ecuaciones lineales Ax = b, podemos destacar lo siguiente: 1. Si m = n, puede pasar lo siguiente: a) Si rango (A) = n, el sistema Ax = b tiene solucin nica, pues A es inversible. En estas condiciones, se dice que el sistema es compatible y determinado. b) Si rango (A) < n, puede pasa lo siguiente 1) Si b Im(A), el sistema Ax = b admitir innitas soluciones. En estas condiciones, se dice que el sistema es compatible e indeterminado. 2) Si b / Im(A), el sistema no admite solucin. En estas condiciones, se dice que el sistema es incompatible.

56

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES 2. Si m > n, aunque rango(A) = n, el sistema puede no tener solucin, debido a que es muy frecuente que b / Im(A). 3. Si m < n, es muy probable que el sistema se compatible e indeterminado si b Im(A). En el caso que b / Im(A), el sistema ser incompatible.
rango(A)=n b Im(A) so l. n ica b / Im(A) so l. n ica Rango b Im(A) innitas soluc. innitas soluc innitas soluc deciente b / Im(A) incompatible incompatible incompatible Cuadro 3.1: Soluciones de un sistema de ecuaciones lineales

Matriz A Rmn Rango completo

m<n
rango(A)=m intas soluc.

m=n
rango(A)=n sol. nica

m>n

El cuadro 3.1 esquematiza las soluciones de un sistema de ecuaciones lineales.

4.2.

Mtodos Numricos para Resolver Sistemas de Ecuaciones Lineales

Los mtodos numricos que estudiaremos en este captulo requieren que el sistema (4.1) est constituido de n las y n columnas, es decir n = m. Al nal de este captulo, en el problema 3.4, ser discutida una estrategia para enfrentar sistemas lineales donde m < n. Por otro lado sistemas donde n > m se denominan sobredeterminados, ellos sern discutidos en el siguiente captulo. Los mtodos para resolver sistemas de ecuaciones lineales de n n pueden ser de dos tipos, directos e iterativos. Directos: Si la solucin existe, otorgan la solucin exacta del sistema lineal despus de un nmero nito de operaciones, excepto errores de redondeo. Iterativo: Dada una aproximacin inicial x0 , generan una sucesin de vectores {xk } k=0 . Si la solucin existe, bajo ciertas condiciones, esta sucesin converge a la solucin.

A simple vista, parece ser que la eleccin de un mtodo directo es la ms adecuada. No obstante, como veremos ms adelante, los mtodos directos tienen un alto costo computacional, requiriendo alrededor de n3 operaciones elementales para resolver un sistema lineal de n n. En contraste, los mtodos iterativos parecen menos atractivos, sin embargo, slo requieren alrededor de n2 operaciones elementales por cada iteracin1 , eso los torna ms viables en la resolucin de sistemas de ecuaciones de grandes dimensiones2 . En las siguientes secciones analizaremos algunos de los mtodos ms importantes para la resolucin numrica de sistemas de ecuaciones lineales.
1 2

Pero pueden requerir muchas iteraciones. Frecuentemente, cuando n > 1000.

4.3. MTODO DE CRAMER

57

4.3.

Mtodo de Cramer

La Regla de Cramer"puede ser enfocada como un mtodo directo. Proviene de un teorema en lgebra lineal, mediante el cual se puede obtener la solucin de un sistema lineal de ecuaciones en trminos de determinantes. Recibe este nombre en honor a Gabriel Cramer (1704-1752). Si Ax = b es un sistema de ecuaciones lineales, donde A Rnn es inversible y b Rn es un vector columna, entonces la solucin del sistema se calcula as: xj = det(Aj ) det(A) j = 1, ..., n (4.3)

donde Aj es la matriz que resulta de reemplazar la j -colimna de A por b. Ejercicio 4.2 El determinante de una matriz real cuadrada de orden n est denido por
n X j +1 det(A) = (1) a1,j det(A [1, j ]) j =1

(4.4)

donde a1,j es el elemento en la 1-la y la j -columna de A, mientras que A [1, j ] es la submatriz obtenida al eliminar la 1-la y la j -columna de A. Muestre que el nmero de multiplicaciones necesarias para hallar det(A) es aproximadamente n! Es posible calcular el determinante de una matriz sin utilizar directamente la denicin dada en (4.4), pero eso tomar por lo menos alrededor de n3 operaciones elementales. Por lo tanto, para nes prcticos, mtodos como el de Cramer deben estar fuera de nuestro interes, debido a que resulta muy costoso en la resolucin numrica de sistemas de ecuaciones lineales3 .

4.4.

Mtodo de Gauss

Este mtodo numrico es el ms conocido y encuadra dentro de los mtodos directos. Dado el sistema lineal de n n a1,1 x1 + a1,2 x2 a2,1 x1 + a2,2 x2 . . . . . . . . . a x + a x n,1 1 n,2 2

+ + . ... . . +

+ + . . .

a1,n xn a2,n xn . . .

= = . . .

b1 b2 . . .

(4.5)

+ an,n xn

= bn

El mtodo de eliminacin de Gauss consiste de transformar el sistema de (4.5), de un modo equivalente, a un sistema triangular superior: a0 x + a01,2 x2 + + a01,n xn = b01 1,1 1 a02,2 x2 + + a02,n xn = b02 (4.6) . . . ... . . . . . . . . . a0n,n xn = b0n
3

La regla de Cramer tiene un valor terico, pues se utiliza en la demostracin de muchas propiedades.

58

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES

El benecio de esto es que se puede resolver el sistema triangular (4.6) de modo eciente, as que de la ltima ecuacin de (4.6) tenemos b0n xn = 0 an,n Luego, xn1 puede ser obtenido mediante xn1 = b0n1 a0n1,n xn a0n1,n1

Y as sucesivamente, obtenemos xn2 , xn3 , ..., x2 , nalmente b01 a01,2 x2 a01,3 x3 ... a01,n xn x1 = a01,1

4.4.1.

Resolviendo un Sistema Triangular

El siguiente algoritmo resuelve un sistema de ecuaciones lineales de orden n, el cual ya est en la forma tringular superior. Algoritmo 4.1 (Resolucin de un sistema triangular superior) Dado un sistema tringular superior Ax = b de orden n, con elementos de A sobre la diagonal no nulos. Los valores de las variables xn , xn1 , ..., x2 , x1 son obtenidos mediante: xn = bn /an,n para k = (n 1), ..,1 s=0 Para j = (k + 1), ..., n s = s + ak,j xj xk = (bk s)/ak,k Ejercicio 4.3 Anlogamente al algoritmo 3.1, disear un algoritmo que resuelva un sistema triangular inferior de orden n. Ejercicio 4.4 Cuntas operaciones elementales (sumas, restas, multiplicaciones, divisiones y comparaciones) son necesarias para la ejecuacin del algoritmo 3.1 y el algoritmo planteado en el ejercicio 3.3? Laboratorio 4.1 En algn lenguaje de programacin, implemente el algoritmo planteado en el ejercicio 3.3.. La conversin de un sistema de orden n a un sistema equivalente, y triangular superior, es posible en virtud al siguiente teorema del lgebra lineal: Teorema 4.2 Sea Ax = b un sistema de ecuaciones lineales. Aplicando sobre las ecuaciones de este sistema una sucesin de operaciones elementales: 1. Cambiar dos ecuaciones 2. Multiplicar una ecuacin por una constante no nula 3. Adicionar un mltiplo de una ecuacin a otra ecuacin = obtenemos un nuevo sistema Ax b, el cual es equivalente4 al sistema original Ax = b.
4

Equivale, en el sentido que tiene las mismas soluciones

4.4. MTODO DE GAUSS

59

4.4.2.

Conversin de un Sistema de Ecuaciones Lineales a uno Triangular Superior

El siguiente algoritmo usa el Teorema 3.1 y convierte un sistema de ecuaciones lineales Ax = b en un sistema triangular superior equivalente. Algoritmo 4.2 Dado el sistema lineal Ax = b de n n. Suponga que el elemento ak,k 6= 0 al inicio de la etapa k : Para k = 1, ..., n 1 Para i = k + 1, ..., n mi,k = ai,k /ak,k ai,k = 0 Para j = k + 1, ..., n ai,j = ai,j mi,k ak,j bi = bi mi,k bk Denicin 4.1 El nmero mi,j = ai,k , ak,k i = 2, ..., n, k = 1, ..., n 1

introduciendo en el algoritmo 3,2, lo denominaremos (i, k )-multiplicador de la matriz A. Adems, ak,k se denomina el k-simo pivote. Ejemplo 4.1 Considere el siguiente sistema 1 2 3 2 2 4 A = A(0) = 1 1 2 1 0 0 Aplicando el Algoritmo 3.2, para y m4,1 = 1, de donde: 1 0 A(1) = 0 0 de ecuaciones lineales Ax = b, donde 0 2 10 (0) 4 y b = b = 6 1 2 0 2 0 = 4 2 2 0 = 4 2 2 0 = 4 3

k = 1, los multiplicadores respectivos son m2,1 = 2, m3,1 = 1 2 3 0 2 2 10 1 5 1 2 3 2

y b(1)

Finalmente, para k = 3 tenemos m4,3 = 0, 25 y 1 2 3 0 0 2 2 10 A(3) = 0 0 4 4 0 0 0 7

para k = 2, m3,2 = 0, 5, m4,2 = 1 y 1 2 3 0 0 2 2 10 A(2) = 0 0 4 4 0 0 1 8

y b(2)

y b(3)

60

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES

Resta slo aplicar el Algoritmo 3.1 para resolver el sistema triangular superior A(3) x = b(3) , de donde obtuvimos 6/7 25/7 x= 10/7 3/7

Un incoveniente que puede suceder en la aplicacin de los algoritmos 3.2 y 3.1, es el clculo del multiplicador mi,k , pues se necesita que ak,k 6= 0 en cada iteracin. Pero el simple hecho que ak,k sea pequeo puede ocasionar que el multiplicador mi,k tome valores inmensamente grandes, los que puede ocasionar a su vez mal condicionamiento5 del sistema.

Ejemplo 4.2 (Sistema mal condicionado) Utilice algn lenguaje de programacin para resolver, usando los algoritmos 3.2 y 3.1, el siguiente sistema de ecuaciones: 11x1 + 2x2 = 5 16 10 x1 + 0, 5x2 = 9 Despus de aplicar consecutivamente los algoritmos, obtuvimos como solucin; 7, 26691434300103 1016 x= 2, 5 Al vericar si realmente es solucin, vemos que 5 5 6= b = Ax = 9 8, 51691434300102 Este fenmeno se debe a que el pivote A1,1 = 11, cuando comparando con A2,1 = 1016 , es muy pequeo. Observe que en estas condiciones, el multiplicador m2,1 = a2,1 1016 = a1,1 11

es muy grande, lo que ocacionar imprecisiones en los clculos siguientes, pues el computador trabaja con precisin nita. Pero qu signica mal condicionamiento ? Ese concepto es de gran ayuda en la resolucin de sistemas lineales. A continuacin formalizaremos ese trmino.

4.5.

Sensibilidad en Sistemas Lineales: Nmero de Condicin

Considere el sistema Ax = b, donde A es inversible y b es un vector no nulo. Digamos que x sea la nica solucin. Consideremos ahora una perturbacin b de b, el sistema lineal perturbado ser Ay = b + b (4.7)
La matriz A tiende a ser no inversible en la prctica, aunque en la teora lo sea. Esto puede producir imprecisiones en los clculos realizados por el computador.
5

4.5. SENSIBILIDAD EN SISTEMAS LINEALES: NMERO DE CONDICIN

61

Observe que el sistema (4.7) tambin tiene nica solucin, digamos que sta sea y . Denotamos x =y x , en consecuencia, y =x + x . Es natural esperar que cuando b sea pequeo, entonces x tambien lo sea. Para cuanticar el tamao de vectores, usaremos la norma vectorial euclidiana6 kk . As, la medida de b relativa a b es kbk / kbk , mientras que la medida de kx k / kx k . Por lo tanto, en trminos ms precisos, esperamos que cuando kbk / kbk sea pequeo, entonces kx k / kx k tambin lo sea. As, como y es solucin del sistema (4.7), entonces Ay = b + b = A( x + x ) = b + b = A(x ) = b = x = A1 (b) Considerando la norma matricial inducida por la vectorial, kAk = m axx6=0 x = A1 (b) = kAx k = A1 (b) A1 kbk b = x = kbk = kAx k kAk kx k = Por (4.8) y (4.9) tenemos kAk 1 kx k kbk
kAxk , kxk

tenemos (4.8)

Por otro lado,

(4.9)

kx k kbk A1 kAk kx k kbk 1 Observe que, si kbk / kbk es pequeo y kA k kAk es un nmero razonablemente pequeo, entonces kx k / kx k debera ser tambin pequeo. Sin embargo, si kA1 k kAk es extremadamente grande, a pesar que kbk / kbk sea pequeo, no hay garanta que kx k / kx k tambin lo sea. En otras palabras, si kA1 k kAk es grande, podemos advertir que el sistema Ax = b puede ser susceptible a grandes alteraciones en la solucin si b es ligeramente perturbado. Esto es crucial, pues el mtodo de Gauss altera el vector b en cada iteracin. Como se puede apreciar, el papel del nmero kA1 k kAk juega un rol importante en la determinacin de la estabilidad del sistema de ecuaciones lineales. A continuacin lo formalizamos. Denicin 4.2 (Nmero de condicin) Dada la matriz A Rnn inversible. El nmero de condicin de A, denotado por cond(A), es cond(A) = A1 kAk

Si cond(A) es grande, diremos que la matriz es mal condicionada, caso contrario, ella ser bien condicionada. Ejercicio 4.5 Mostrar que, si kk es una norma matricial inducida por la norma vectorial eucldea en Rn , cond(A) 1. Solucin. Observe que la norma matricial usada es consistente, as: I = A1 A = 1 = kI k = A1 A A1 kAk = cond(A)

Observacin 4.1 Algunas veces se utiliza tambin la siguiente expresin, rcond(A), la cual est denida por rcond(A) = 1/ cond(A) Claramente, si rcond(A) es prximo a cero, la matriz A ser mal condicionada.
6

Si x R, kxk =

p 2 x2 1 + ... + xn .

62

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES

4.6.

Estrategia con Pivoteamiento Parcial

Como habamos comentado enteriormente, uno de los inconvenientes en la aplicacin del mtodo de Gauss consista en el clculo de los multiplicadores, si el multiplicador mi,k = ai,k /ak,k es demasiado grande, ste podra ocasionar clculos imprecisos en el sistema de cmputo. Algunas de las estrategias para evitar tal inconveniente son: 1. Estrategia con pivoteamiento parcial 2. Estrategia con pivoteamiento total Aqu veremos en qu consiste la estrategia de pivoteamiento parcial: En el inicio de ta etapa k de la eliminacin gaussiana, escoger como pivote el elemento de k1 k1 mayor mdulo entre los coecientes: ai,k , para i = k, k + 1, ..., n, donde ai,k representa la (i, k )-componente en la iteracin k 1. Permutar la la k e i, si fuera necesario. Ejemplo 4.3 Considere un sistema de orden n, donde n = 4 y la iteracin k = 2. Observe que 3 2 1 1 5 (1) (1) 0 1 0 3 6 A b = 0 3 5 7 7 0 2 4 0 15 Note que A(1) b(1) representa la situacin del sistema en la primera iteracin. Ahora, para el inicio de la segunda iteracin 1 1. Escoger el pivote: m axj =2,3,4 a1 j,2 = a3,2 = 3, entonces el pivote es 3

2. Cambiamos las lneas 2 y 3, de donde obtenemos 3 2 1 1 5 (1) (1) 0 3 5 7 7 b = A 0 1 0 3 6 0 2 4 0 15 1 1 = 3 3 2 2 = = 3 3

Observe que en este caso los multiplicadores respectivos sern: m3,2 = m4,2

As, al escoger el mayor elemento en mdulo entre los candidatos a pivote, se consigue que los multiplicadores m, en mdulo, estn entre cero y uno, lo que evita la propagacin de errores de redondeo.

4.7. DESCOMPOSICIN LU

63

La estrategia con pivoteamiento parcial no elimina del todo la acumulacin de errores de redondeo, existe otra estrategia denominada estrategia con pivoteamiento completo. En contraste al pivoteamiento parcial, que busca el mejor pivote en una porcin de cada columna en cada iteracin, la estrategia de pivoteamiento completo analiza toda matriz. A pesar que en teora esto elimina dinitivamente las imprecisiones numricas que puedan ocurrir, su uso no es comn en la prctica pues requiere mucho esfuerzo computacional, es decir, requiere muchas operaciones elementales (comparaciones) para su ejecucin. Problema 4.1 Investigue la estrategia de pivoteamiento completo.

4.7.

Descomposicin LU

Esta estrategia consiste en descomponer la matriz A como el producto de dos matrices factores: A = LU, donde L es una matriz triangular inferior con unos sobre su diagonal, es decir 1 0 0 2,1 1 0 L= . . . . . . . . . . . . 1 n,1 n,2 y U es una matriz triangular superior U = u1,1 0 . . . 0 u1,2 u2,2 . . . 0 ... u1,n u2,n . . . un,n

De ese modo, resolver el sistema Ax = b, o sino LU x = b, es equivalente a resolver consecutivamente: 1. Ly = b 2. Ux = y Una de la primeras ventajas que tiene resolver sistemas mediante este mtodo, en contraste con el mtodo de Gauss, es la siguiente: supongamos que se quiere resolver los sistemas Ax = b y Ax = b0 , una vez obtenidos tales factores L y U, bastara resolver consecutivamente los sistemas triangulares: Ly = b Ux = y (4.10) y Ly = b0 Ux = y (4.11)

Observe que los mismos factores L y U fueron usados para resolver los sistemas (4.10) y (4.11), los cuales son fciles de resolver usando el algoritmo 3.1 para sistemas triangulares superiores.

64

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES

4.7.1.

Clculo de los Factores L y U

El clculo de las matrices L y U de modo que A = LU, esta basado en los multiplicadores mi,j introducidos en las eliminacin gaussiana (denicin 3.1). Supongamos que tenemos la matriz 0 a0 a1,1 a0 1,2 1,3 =A a0 a0 A0 = a0 2,1 2,2 2,3 0 0 0 a3,1 a3,2 a3,3 Primera iteracin: Los respectivos multiplicadores estn dados por
0 m2,1 = a0 2,1 /a1,1 0 m3,1 = a0 3,1 /a1,1

para deshacernos de la variable x1 en las i-las, para i = 2, 3, hacemos

0 a1 1,j = a1,j

j = 1, 2, 3 i = 2, 3, a1 1,2 a1 2,2 a1 3,2 a1 1,3 a1 2,3 1 a3,3 0 0 1 0 0 1 j = 1, 2, 3

0 0 a1 1,j = a1,j mi,1 a1,j

y obtenemos

Observe que A1 = M 0 A0 , donde

a1 1,1 1 A = 0 0

1 0 M = m2,1 m3,1

Segunda iteracin: Para eliminar x2 de las restantes lineal, calculamos el multiplicador 1 m3,2 = a1 3,2 /a2,2 y hacemos
1 a2 1,j = a1,j

j = 1, 2, 3 j = 2, 3 m3,2 a1 2,j a2 1,2 a2 2,2 0 j = 1, 2, 3 a2 1,3 a2 2,3 2 a3,3 0 0 1

a2 2,j a2 3,j y obtenemos

= =

a1 2,j a1 3,j

Observe que donde

a2 1,1 2 A = 0 0

A2 = M 1 A1 1 0 1 1 M = 0 0 m3,2

Adems, A2 es triangular superior y podemos denir U = A2 . Luego: U = M 1 A1 = M 1 (M 0 A0 ) = (M 1 M 0 )A0 = (M 1 M 0 )A

4.7. DESCOMPOSICIN LU es decir, A = (M 1 M 0 )1 U = (M 0 )1 (M 1 )1 U y como (M 0 )1 (M 1 )1 1 1 1 0 0 1 0 0 1 0 = m2,1 1 0 0 m3,1 0 1 0 m3,2 1 1 0 0 1 0 0 1 1 0 = m2,1 = m2,1 1 0 0 m3,1 0 1 m3,1 0 m3,2 1

65

0 1 m3,2

de modo que A = LU :

concluimos que la matriz triangular inferior debera ser 1 0 0 1 0 L = m2,1 m3,1 m3,2 1

0 0 1

Teorema 4.3 Dada una matriz cuadrada A de orden n, sea la matriz Ak constituida de las primeras k lneas y k columnas de A. Suponga que det(Ak ) 6= 0 para k = 1, ..., n 1. Entonces, existe una nica matriz triangular inferior L = [mi,j ] , con mi,i = 1, para i = 1, ..., n, y una nica matriz triangular superior U = [ui,j ] tales que A = LU. Adems, det(A) = u1,1 u2,2 ...un,n . Ejercicio 4.6 Investigue la demostracin del teorema 3.2. Ejemplo 4.4 Resolver el siguiente sistema de ecuaciones lineales utilizando la descomposicin LU. x1 + 2x2 + 3x3 = 10 2x1 + 5x2 x3 = 20 x1 + 2x2 + x3 = 6 Observe que 1 2 3 A = 2 5 1 1 3 1 m2,1 = 2, y m3,2 = 4, Por lo tanto m3,1 = 1, y 10 b = 20 6

Luego,

1 2 A= 1

2 3 5 1 3 1 1 2 = 1

1 2 3 A2 = 0 1 7 0 0 32

1 2 3 A1 = 0 1 7 0 4 4

1 = m2,1 m 3,1 0 0 1 1 0 0 4 1 0

0 0 1 0 A2 = m3,2 1 2 3 1 7 = LU 0 32

66

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES

Puesto que Ax = b lo podemos expresar como (LU )x = b, usando el algoritmo 3.1 y el ejercicio 3.3 para resolver sistemas triangulares, calculamos consecutivamente: 10 Ly = b = y = 0 16 y 3/2 U x = y = x = 7/2 1/2

4.8.

Descomposicin de Cholesky

Una estrategia que tiene suma importancia, sobre todo en optimizacin, es la descomposicin de Cholesky. Denicin 4.3 (Matriz denida positiva) Se dice que una matriz A Rnn simtrica y de orden n es denida positiva, si xt Ax > 0, para todo x Rn y x 6= 0. Teorema 4.4 (Descomposicin de Cholesky) Si A es una matriz de orden n, simtrica y denida positiva, entonces existe una nica matriz triangular inferior n y con diagonal positiva, tal que A = GGt Ejercicio 4.7 Investigue la demostracin del teorema 3.3.

4.8.1.

Clculo del factor de Cholesky

g1,1 0 . . . 0 g1,2 g2,2 . . . 0 ... g1,n g2,n . . . gn,n

La columna 1:

El factor G ser obtenido de la ecuacin matricial: g1,1 a1,1 a1,2 a1,n 0 a2,1 a2,2 a2,n g2,1 g2,2 . = . . . . ... ... . . . . . . . . . . gn,1 gn,1 an,1 an,1 an,n a1,1 a2,1 . . . an,1 g1,1 g2,1 . . . gn,1 0 g2,2 . . . gn,1 ... 0 0 . . . gn,n

Dada un matriz simtrica y denida positiva: a1,1 a1,2 a2,1 a2,2 A= . . . . . . an,1 an,1

...

a1,n a2,n . . . an,n

0 0 . . . gn,n

g1,1 0 . . . 0

g1,1 g1,1 g2,1 g1,1 . . . gn,1 g1,1

4.8. DESCOMPOSICIN DE CHOLESKY Luego g1,1 = La columna 2: a1,2 a2,2 a3,2 . . . an,2 de donde a1,1 aj,1 g1,1

67

gj,1 =

para j = 2, ..., n

g1,1 g2,1 = . . . gn,1

0 g2,2 . . . gn,1

...

0 0 . . . gn,n

g1,2 g2,2 0 . . . 0

g1,1 g1,1 2 2 g2,1 + g1 ,1 g3,1 g2,1 + g3,2 g2,2 . . . gn,1 g1,1 + gn,2 g2,2

Como

a1,2 g1,2 g1,1 = a1,2 = g1,2 = g2,1 = g1,1 q 2 2 2 g2 g2,2 = a2,2 g2 ,1 ,1 + g1,1 = a2,2 = gj,1 g2,1 + gj,2 g2,2 = aj,2 j = 3, ..., n

y los gj,1 = 1, ..., n, ya fueron calculados, entonces gj,2 = (aj,2 gj,1 g2,1 ) g2,2 j = 3, ..., n

columna k : Los elementos de la columna k de G son: [0 0 gk,k gk+1,k gn,k ]t haciendo ak,k ak+1,k . . . an,k ak,1 ak,2 . . . g1,1 g2,1 = . . . gn,1 k = 2, ..., n gk,1 gk,2 . . .

0 g2,2 . . . gn,1

...

0 0 . . . gn,n

tenemos de donde

gk,k 0 . . . 0

2 2 2 ak,k = gk, 1 + gk,2 + + gk,k

gk,k = Y como

k1 X 2 ak,k gk,j i=1

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

y como los elementos gi,k , i = 1, ...., k 1, ya fueron calculados, tenemos Pk1 aj,k i g g j,i k,i =1 gj,k = j = k + 1, ..., n gk,k

aj,k = gj,1 gk,1 + gj,2 gk,2 + + gj,k gk,k

68

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES

Algoritmo 4.3 Sea A una matriz simtrica y denida positiva. para k = 1, ..., n suma = 0; Para j = 1, ...(k 1) 2 suma = suma + gk,j r = ak,k suma............( ) gk,k = r para i = k + 1, ..., n suma = 0 para j = 1, ..., k 1 suma = suma + gi,j gk,j gi,k = (ai,k suma)/gk,k Por lo general, ver si una matriz simtrica es denida positiva usando la denicin es una tarea prcticamente imposible. Sin embargo, podemos usar el algoritmo de Cholesky para vericar si A es denida positiva. Si en () se tiene que r 0, entonces la descomposicin no es posible y A no es denida positiva. Caso contrario, el algoritmo otorga la matriz triangular inferior G tal que A = GGt .

La descomposicin de Cholesky requiere alrededor de n3 /6 operaciones de multiplicacin para la descomposicin. Este nmero es aproximadamente la mitad del nmero de operaciones necesarias para la realizacin de la eliminacin en la descompisicin LU, pues requera n3 /3. Al igual que la descomposicin LU , una vez conocido A = GGt , es posible resolver el sistema lineal asociado: Ax = b, del siguiente modo: como Ax = b GGt x = 0 entonces 1. Resolver Gy = b 2. Resolver Gt x = y Laboratorio 4.5 En algn lenguaje de programacin de su preferencia, implemente el algoritmo de Cholesky, de modo que: cuando la matriz A sea denida positiva, devuelva la matriz triangular inferior G. Caso contrario, debera emitir un mensaje anunciando que A no es denida positiva. Ejercicio 4.8 Usando el algoritmo de Cholesky. Verique que la matriz 27 12 10 24 12 9 16 29 10 16 42 70 24 29 70 137 es realmente denida positiva.

4.9. MTODOS ITERATIVOS

69

es o no denida positiva (use el algoritmo)

Ejercicio 4.9 Una duda comn es, Podra ser denida positiva una matriz la cual tiene algunas componentes negativas? Verique si la matriz 27 12 10 24 12 9 16 11 10 16 42 2 24 11 2 137

Ejercicio 4.10 Una matriz simtrica A Rnn se llama semidenida positiva, si xt Ax 0, para todo x Rn . Probar que, para cualquier matriz B Rmn , la matriz C = B t B es semidenida positiva. Solucin: Observe que C Rmm es simtrica, pues C t = (B t B )t = B t (B t )t = B t B = C y donde kk es la norma euclidea. xt Cx = xt (B t B )x = (Bx)t (Bx) = kBxk2 0

Ejercicio 4.11 Probar que, si A Rnn es denida positiva, entonces A es inversible. Lo recproco no siempre es vlido, de un contraejemplo. Solucin: Recordemos el teorema de lgebra lineal que nos deca: na matriz A Rnn es inversible si, y slo si, el sistema Ax = 0 tiene como nica solucin la trivial". Ahora, con respecto al ejercicio, realizaremos una prueba por contradiccin. Asumamos que A es denida positiva y no inversible. Entonces, el sistema Ax = 0 tiene una solucin no trivial, digamos que existe x 6= 0 tal que Ax = 0. Luego. 0<x t Ax =x t 0 = 0 Lo cual es una contradiccin. Esto completa la prueba. Para vericar la segunda parte del ejercicio, basta hacer notar que la matriz 1 0 0 1 es inversible pero no es denida positiva (en realidad es denida negativa).

4.9.

Mtodos Iterativos

La idea de los mtodos iterativos para resolver sistemas de ecuaciones est inspirada en el mtodo de punto jo. Dado el sistema lineal de ecuaciones Ax = b (4.12)

donde A Rnn y b Rn . Los mtodos iterativos consisten en expresar el sistema (4.12) en una ecuacin matricial equivalente: x = Dx + d (4.13)

70

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES

donde D Rnn y d Rn . Note que en este caso (x) = Dx + d es una funcin iteracin n-dimensional. As, los mtodos de este tipo construyen una sucesin {xk Rn } k=0 denida por la regla xk+1 = (xk ) = Dxk + d Observe que, si l mk xk = x , entonces x es un punto jo de . Bajo ciertas condiciones, x debera ser tambin la solucin del sistema original Ax = b. Mtodo de Gauss-Jacobi Este mtodo se caracteriza por transformar el sistema Ax = b a la forma x = Dx + d, del siguiente modo: Sea el sistema original a1,1 x1 a2,1 x1 . . . an,1 x1 + + . . . a1,2 x2 a2,2 x2 . . . + + . ... . . + + + . . . a1,n xn a2,n xn . . . = = . . . b1 b2 . . .

+ an,2 x2

+ an,n xn

= bn

Suponiendo que ai,i 6= 0 para i = 1, ..., n, despejamos las variables x1 , ..., xn de las n ecuaciones, respectivamente: x1 = x2 1 (b1 a1,2 x2 a1,3 x3 ... a1,n xn ) a1,1 1 = (b2 a2,1 x2 a2,3 x3 ... a2,n xn ) a2,2 . . . 1 = (bn an,1 x1 an,2 x2 ... an,n1 xn1 ) an,n

x1 De esta forma tenemos x1 x2 . = . . xn

b1 a1,1 b2 a2,2

. . .

bn an,n

a2,1 a2,2 + . . . an,1 an,n

1,2 a a1,1 0 . . .

1,3 a a1,1 2,3 a a2,2 . . .

...

1,n a a1,1 2,n a a2,2 . . .

x1 x2 . . . xn

an,2 a n,n

an,3 a n,n

Es decir, la funcin iteracin quedara establecido por x = (x) = Dx + d donde 0

1,2 a a1,1 0 . . . 1,3 a a1,1 2,3 a a2,2 . . .

a2,1 a2,2 D= . . . an,1 an,n

...

1,n a a1,1 2,n a a2,2 . . .

an,2 a n,n

an,3 a n,n

y d=

b1 a1,1 b2 a2,2

. . .

(4.14)

bn an,n

4.10. COMPARACIN ENTRE MTODOS DIRECTOS E ITERATIVOS

71

El mtodo de Gauss-Jacobi consiste en lo siguiente: dado x0 , una aproximacin inicial a la solucin del sistema lineal Ax = b, obtener x1 , ..., xk , ... por medio de la relacin recursiva xk+1 = Dxk + d donde D y d estn denidos en este caso por (4.14). O sino:
+1 = xk 1 k+1 x2 = 1 (b a1,1 1 1 (b a2,2 2 k k a1,2 xk 2 a1,3 x3 ... a1,n xn ) k k a2,1 xk 2 a2,3 x3 ... a2,n xn ) . . .

(4.15)

+1 = xk 1

1 (b an,n n

k k an,1 xk 1 an,2 x2 ... an,n1 xn1 )

Una caracterstica de los mtodos iterativos, a diferencia de los mtodos directos, es que slo convergen si algunas hiptesis son satisfechas. Criterio de convergencia para el mtodo de Gauss-Jacobi Dado el sistema lineal Ax = b, sea Pn i =

j =1 j 6=i

Si < 1, entonces el mtodo de Gauss-Jacobi genera una sucesin de vectores xk k=0 , la cual converge a la solucin del sistema lineal Ax = b, independientemente de la eleccin del vector inicial x0 .

|ai,i |

|ai,j |

y = m ax {i }
1in

4.10.

Comparacin entre Mtodos Directos e Iterativos

Convergencia: Los mtodos directos son procesos nitos y, tericamente, encuentran la solucin exacta de cualquier sistema de ecuaciones, desde que sta exista. En la prctica no es bien cierto esta ltima armacin, pues debemos tener en cuenta los errores de clculo del sistema de cmputo y el mal condicionamiento del sistema. Por otro lado, los mtodos iterativos apenas convergen cuando son satisfechos algunos requerimiento. Esto ltimo parece tornarlos poco atractivos en la prctica, pero debemos tener en cuenta el espacio requerido en memoria-computador, observe por ejemplo la ecuacin (4.15), notar que para operar el mtodo de Gauss-Jacobi, slo necesitamos almacenar en memoria la matriz A y el vector b, el resto de clculo no requieren la construccin de matrices auxiliares como en el mtodo directo basado en la descomposicin LU Errores de Redondeo: Mtodos directos modican el sistema original, por ejemplo, en la triangularizacin hecha por el mtodo de Gauss, adicionalmente a los errores cometidos por un sistema de cmputo. Los datos del problema original son alterados y la solucin no es precisa desde un punto de vista computacional. Por otro lado, los mtodos iterativos no modican el sistema original y, una vez asegurada la convergencia, los errores numricos son menores a los mtodos directos. Problema 4.2 Investigue la prueba del criterio de convergencia del mtodo de Gauss-Jacobi. Problema 4.3 Investigue sobre el mtodo iterativo de Gauss-Seidel.

72

CAPTULO 4. RESOLUCIN DE SISTEMAS LINEALES

Problema 4.4 Elaborar una estrategia para resolver un sistema de ecuaciones lineales indeterminado, es decir, donde el nmero de ecuaciones es menor que el nmero de variables (m < n). Una de la opciones sera la siguiente: Asumamos que el sistema est denido por Ax = b, donde A Rmn , b Rn y m < n. Inicialmente, asumamos que el rango de A es completo, es decir, m. Por ese motivo, existe una submatriz cuadrada e inversible B , tal que, posiblemente despus de un intercambio compatible de posiciones de columnas deA y de variables de x, tenemos A = [B N ] xB x = xN As. Ax = b es equivalente a xB = b BxB + NxN = b xB = B 1 b B 1 NxN [B N ] xN Luego, las soluciones x Rn seran generadas dando valores arbitrarios a x N Rnm , mediante: 1 B b B 1 NxN x B = x = x N x N As por ejemplo, haciendo x N = 0 Rnm , tendramos: 1 B b x = 0

el cual claramente, es una solucin para el sistema Ax = b. Se debe advertir que, a diferencia de la resolucin de sistemas cuadrados, la estrategia expuesta no es eciente, pues encontrar las columnas de A que conformarn B puede ser extremadamente costoso. Esto se debe a que resolver sistemas indeterminados, y tambin sobredeterminados, es realmente difcil. El estudiante debera investigar otras opciones adems de la expuesta aqu.

Captulo 5 Introduccin a Sistemas No Lineales

En este captulo revisaremos brevemente el clsico mtodo de Newton para resolver sistemas de ecuaciones no lineales. Este mtodo es iterativo y, en su versin bsica, converge a un a solucin del sistema slo cuando el punto inicial x0 es tomado lo sucientemente prximo de . Cuando esto sucede, se dice que el mtodo tiene la propiedad de convergencia local, tal requerimiento es muy discriminatorio y frecuentemente el mtodo converge inclusive cuando el punto inicial est distante de la solucin. Un ejemplo de un sistema de 3 ecuaciones no lineales con 3 incgnitas es una expresin matemtica de la siguiente forma 7x1 x2 + 5x2 x2 3 sen x1 12 = 0 4 2 x1 + cos x2 + 2x3 38 = 0 6x1 + 2x2 x3 + 34 = 0 (5.1)

Resolver problemas de este tipo no es un asunto fcil, sobre todo cuando el sistema en cuestin es de gran dimensin. Observe que el sistema (5.1) puede ser representado por: f1 (x) = 0 f2 (x) = 0 f3 (x) = 0 donde fi : R3 7 R, i = 1, 2, 3. Si denimos f1 (x) F (x) = f2 (x) f3 (x) F (x) = 0 (5.2)

donde Fi : R3 7 R3 , entonces el sistema de ecuaciones (5.1) puede ser visto a su vez como (5.3) La ecuacin (5.3) es la forma en que nosotros identicaremos un sistema de ecuaciones no lineal. El Jacobiano de F en el punto x, denotado por F 0 (x) o J (x), est denido por f1 f1 f1 (x) x ( x ) ( x ) x1 x n 2 f2 (x) f2 (x) f2 (x) x x x n 2 J (x) = F 0 (x) = 1. . . . . . . . . . . .
fn (x) x1 fn (x) x2

fn (x) xn

73

74

CAPTULO 5. INTRODUCCIN A SISTEMAS NO LINEALES

fi donde x : Rn 7 R, i, j = 1, ..., n, son derivadas parciales y las asumimos como funciones j continuas en Rn .

5.1.

El Mtodo de Newton

El Mtodo de Newton consiste en resolver un problema difcil mediante sucesivas resoluciones de un problema fcil, se espera que cada solucin del problema fcil sea una mejor aproximacin para el problema difcil. As, sea x0 Rn una aproximacin inicial a una solucin del sistema (problema difcil) F (x) = 0 denido en (5.3), si tomamos una funcin L0 tal que L0 (x) F (x), para todo x en una vecindad de x0 , entonces esperamos que la solucin del sistema L0 (x) = 0 (problema fcil) sea una mejor aproximacin que x0 para una solucin de F (x) = 0. El mtodo de Newton es de carcter iterativo y est basado en la aproximacin lineal (segn Taylor) de la funcin F en torno al punto actual x0 : L0 (x) = F (x0 ) + J (x0 )(x x0 ) El prximo punto, x1 , es denido como la solucin de L0 (x) = 0 Si J (x0 ) es inversible, entonces (5.5) tiene solucin nica: x1 = x0 J 1 (x0 )F (x0 ) Visto de otro modo, una iteracin Newton consiste en calcular x1 previamente conocido x0 , mediante la resolucin consecutiva del sistema lineal: J (x0 )d0 = F (x0 ) x1 = x0 + d0 (5.6) (5.7) (5.5) (5.4)

Lo que dice (5.7) es que una vez conocido el paso de newton d0 , es posible calcular x1 , Con esto, podemos utilizar x1 como nuevo punto y calcular x2 y repetir este procedimiento inicial para k = 0, 1, 2, 3, ...hasta que se cumpla F (xk ) < , donde > 0 es una precisin eseada y kk es alguna norma en Rn . Observacin 5.1 De un modo directo, el mtodo de Newton puede ser visto del siguiente modo: calcular xk+1 una vez conocido xk , mediante la siguiente regla: hasta que F (xk ) < . xk+1 = xk J 1 (xk )F (xk ), k = 0, 1, 2, 3, ...

(5.8)

En realidad, obtener xk+1 directamente desde (5.8) debera ser evitado para nes computacionales, ya que calcular J 1 (xk ) es considerado costoso. Por otro lado, en (5.6) podran usarse tcnicas ms baratas para resolver el sistema sin tener que calcular explicitamente J 1 (xk ), una alternativa sera usar el mtodo de Gauss, la descomposicin LU, o la Cholesky si J (xk ) es denida positiva. En conclusin, el mtodo de Newton es aplicado originalmente para resolver un sistema de ecuaciones F (x) = 0, donde F : Rn 7 Rn tiene funciones componentes fi , i = 1, ..., n, las cuales admiten derivadas parciales continuas en Rn . El siguiente algoritmo resume el mtodo.

5.1. EL MTODO DE NEWTON

75

Algoritmo 5.1 (Algoritmo Bsico de Newton) Sea x0 Rn un punto inicial lo sucientemente cerca de la solucin y > 0 el parmetro de precisin deseado: Paso 1: Si F (xk ) < , detenerse, xk es la aproximacin buscada. Caso contrario, ir al paso 2. Paso 2: Resolver el sistemas J (xk )dk = F (xk ) y obtener el paso de newton dk . Ir al paso 3. Paso 3: Calcular xk+1 = xk + dk , hacer k k + 1. Volver al paso 1. Naturalmente, el sistema F (x) = 0 puede no tener solucin, entonces es posible modicar el algoritmo 4.1 para que nalice despus de un nmero determinado de iteraciones indicando la posibilidad de infactibilidad. Ejemplo 5.1 Consideremos el siguiente sistema de ecuaciones no lineales de 3 incgnitas y 3 ecuaciones: 7x1 x2 + 5x2 x2 3 sen x1 12 = 0 4 2 x1 + cos x2 + 2x3 38 = 0 6x1 + 2x2 x3 + 34 = 0 Observe que 7x1 x2 + 5x2 x2 3 sen x1 12 2 3 F (x) = x4 1 + cos x2 + 2x3 8 6x1 + 2x2 x3 + 34

(5.9)

7x2 x2 7x1 + 5 2x3 sen x1 3 cos x1 4x3 2 cos x2 sen x2 6x2 J (x) = 1 3 6 2 1 10 0 20 y como parmetro de precisin = 109 , una Si usamos como punto inicial x = 50 implementacin computacional bsica nos otorga: 4, 23134959407946 x13 = 1, 56752981158965 5, 47684281234392 Es fcil vericar que kF (x13 )k 5, 293308530720633 108 , donde kk es la norma eucldea. Esto indica que x13 es una buena aproximacin para una raz de F (x) = 0.

5.1.1.

Convergencia y Rapidez del Mtodo de Newton

El mtodo de Newton para sistemas n-dimensionales posee propiedades similares al mtodo estudiado con ms detalles en el captulo 2, el cual resolva sistemas de una variable. Resultados tericos muestran que, si el punto incial x0 Rn es elegido lo sucientemente cerca de la solucin del sistema F (x) = 0, entonces el mtodo de Newton converge y lo hace con una tasa de convergencia cuadrtica. Frecuentemente, para caracterizar este hecho uno utiliza la siguiente expresin: el mtodo de Newton es localmente cuadrtico. Los detalles sobre este tema no sern tratados en este curso.

76

CAPTULO 5. INTRODUCCIN A SISTEMAS NO LINEALES

Laboratorio 5.1 En algn lenguaje de programacin de su preferencia, realice una funcin que ejecute el Algoritmo de Newton. Laboratorio 5.2 Conocer J (xk ) presupone el clculo y evaluacin de las primeras derivadas parciales de las funciones fi (xk ), para i = 1, ..., n, un procedimiento analtico que debemos hacer siempre que sea posible. Esa tarea puede a veces ser considerado un inconveniente debido al esfuerzo de clculo y al probable error humano, una alternativa sera usar una aproximacin mediante diferencias nitas. Ejercicio 5.1 Hallar dos raices ms del sistema no lineal dado en (5.9) Ejercicio 5.2 Al aplicar el mtodo de Newton para resolver F (x) = 0, donde F (x) = Ax b, A Rnn , b Rn

Cuntas iteraciones debera hacer el mtodo de Newton? Por qu?. Ejercicio 5.3 Encuentre un punto crtico de f, es decir, resuelva Of (x) = 0, donde 1 f (x) = a + bt x + xt Qx 2 y donde Q Rnn es simtrica y denida positiva, b Rn y a R. En ese punto f alcanza un mnimo?Porqu? Ejercicio 5.4 Sean 1 2 5 1 3 6 2 5 8 25 7 8 2 3 5 9 7 1 10 200 7 6 12 15

Resolver

A=

y b=

Minimizar 1 kAx bk2 2 x R3

1 2

donde kk es la norma eucldea. Sugerencia: Basta hallar un punto crtico de f : R3 7 R, donde f (x) = Observe que f (x) = At Ax At b

kAx bk2 . (5.10)

Adems de usar el mtodo de Newton para resolver (5.10), use tambin un mtodo directo para sistema de ecuaciones lineales y compare sus resultados.

Captulo 6 Interpolacin y Aproximacin

Es frecuente encontrar libros de clculos elemental funciones representando el comportamiento de algn hecho real, lo primero que viene a la mente es la siguiente pregunta de donde sali esa funcin? Algunas de estas funciones se obtienen de la solucin de una ecuacion diferencial que representa algn fenmeno fsico, pero algunos tambin resultan de simples observaciones. Suponga que disponemos de los siguientes datos, resultados de observaciones: (0, 3,42), (1, 3,45), (2, 3,47), (3, 3,5), (4, 3,48), (5, 3,46), (6, 3,44), (7, 3,39), (8, 3,49), (9, 3,56) donde el par (ti , ci ) representa el precio del dlar ci (soles) en el tiempo ti (horas). Estos datos pueden ser organizados del siguiente modo ti ci 0 1 2 3 4 5 6 7 8 9 3,42 3,45 3,47 3,5 3,48 3,46 3,44 3,39 3,49 3,56

Supngase que se quiere calcular: 1. El precio del dlar en el tiempo t = 3,5h 2. La hora en que el precio del dlar est en 3,465 soles El estudio de la interpolacin y aproximacin adems de responder a estos requerimientos, nos permite encontrar funciones que representen estos datos, estas funciones frecuentemente son continuas y hasta diferenciables, esto a su vez nos permitir manejar matemticamente un problema que consista de apenas algunos datos.

6.1.

Interpolacin

Interpolar una funcin consiste en aproximar una funcin f por otra funcin g, escogida dentro de una clase denida de antemano y satisfaciendo algunas propiedades. La funcin g es entonces usada en reemplazo de f . Considere n + 1 puntos distintos: x0 , x1 , ..., xn y los valores respectivos de f en esos puntos: f (x0 ), f (x1 ), ..., f (xn ). La forma de interpolacin de f que veremos a seguir consiste en 77

78 determinar una funcin g tal que

CAPTULO 6. INTERPOLACIN Y APROXIMACIN

g (x0 ) = f (x0 ) g (x1 ) = f (x1 ) . . . g (xn ) = f (xn ) Una manera prctica y muy utilizada es considerar la funcin interpolante g como polinomio. No obstante, g puede ser una funcin racional, trigonomtrica, etc. En la gura 5.1 se muestra un ejemplo para el caso n = 5.

Figura 5.1: La funcin g interpolando f en 6 puntos

6.1.1.

Interpolacin Polinomial

Dados n + 1 pares (x0 , f (x0 )), (x1 , f (x1 )), ..., (xn , f (xn )), queremos aproximar f por un polinomio pn de grado menor o igual a n, tal que f (xk ) = p(xk ), k = 0, 1, ..., n (6.1)

Naturalmente, podemos preguntarnos ahora: existe siempre un polinomio que satisface tales condiciones?, y si existe, el polinomio es nico?. Para responder a esto, representaremos el polinomio p por p(x) = an xn + an1 xn1 + ... + a1 x + a0 As, al encontrar los valores an , an1 , ..., a1 , a0 :
n1 + ... + a1 x0 + a0 = f (x0 ) p(x0 ) = an xn 0 + an1 x0 n n1 p(x1 ) = an x1 + an1 x1 + ... + a1 x1 + a0 = f (x1 ) . . . n1 + ... + a1 xn + a0 = f (xn ) p(xn ) = an xn n + an1 xn

Observe que el sistema anterior es de la forma A = b (6.2)

6.1. INTERPOLACIN donde A= xn 0 xn 1 . . . xn n

1 xn 0 1 xn 1 . . . 1 xn n

79

...

x0 x1 . . . xn

1 1 . . . 1

Note que es el vector de incgnitas. La matriz A se le conoce como la matriz de Vandermonde, y por lo tanto, si x0 , x1 , ..., xn son todos distintos, entonces Y (xi xj ) 6= 0 det(A) =
0i<j n

n n1 . . . 0

y b=

f (x0 ) f (x1 ) . . . f (xn )

Esto signica que A es invertible y, en consecuencia, el sistema lineal (6.2) tiene solucin y sta es nica. Esto constituye la prueba del siguiente teorema: Teorema 6.1 Existe un nico polinomio pn de grado menor o igual a n, tal que f (xk ) = pn (xk ), desde que xk 6= xj , k 6= j, k, j = 0, 1, ..., n. k = 0, 1, ..., n

6.1.2.

Modos de Obtener el Polinomio pn

Una vez que sabemos las condiciones que garantizan la existencia y unicidad de pn , procedemos a ver las maneras cmo encontrarlo. Existen varias formas de realizar esta tarea, la primera es mediante la solucin directa del sistema A = b. No obstante, existen otras formas, como la de Lagrange y de Newton. Desde un punto de vista terico, todas ellas llevan al mismo polinomio, pero la eleccin de una u otra depende de algunos factores, tales como la dicultad de los clculos, estabilidad del sistema lineal y tiempo de ejecucin. Obtencin de pn Mediante la Resolucin A = b Consiste en resolver por algn mtodo, directo o iterativo, el sistema de ecuaciones lineales en (6.2). Ejemplo 6.1 Vamos a interpolar f en los puntos x0 = 2, x1 = 0, x2 = 1 y x3 = 4, usando un polinomio p3 (x) = a3 x3 + a2 x2 + a1 x + a0 y los siguientes datos: xi 2 0 1 4 f (xi ) 6 4 5 8 podemos calcular la matriz de Vandermonde, pero vamos a repetir el procedimiento: p3 (2) p3 (0) p3 (1) p3 (4) = = = = a3 (2)3 + a2 (2)2 + a1 (2) + a0 = 6 = f (2) a3 (0)3 + a2 (0)2 + a1 (0) + a0 = 4 = f (0) a3 (1)3 + a2 (1)2 + a1 (1) + a0 = 5 = f (1) a3 (4)3 + a2 (4)2 + a1 (4) + a0 = 8 = f (4)

80 Luego

CAPTULO 6. INTERPOLACIN Y APROXIMACIN 1 1 , 1 1 3 2 = 1 0 6 4 y b= 5 8

La solucin puede ser obtenida resolvindo el sistema lineal A = b, por cualquier tcnica ya estudiada. Obtenemos 3 1/9 2 5/9 1 = 5/9 4 0 1 5 5 p3 (x) = x3 + x2 + x + 4 9 9 9 La gura 5.2 muestra la situacin de este ejemplo. de donde

(2)3 (0)3 A= (1)3 (4)3

(2)2 (0)2 (1)2 (4)2

(2) (0) (1) (4)

Figura 5.2: Polinomio p3 interpolando 4 puntos Obtencin de pn Mediante la Forma de Lagrange Dados n + 1 puntos diferentes x0 , x1 , ..., xn , donde y0 = f (x0 ), y1 = f (x1 ), ..., yn = f (xn ). Sea pn un polinomio de grado menor o igual que n, tal que interpola a f en los puntos x0 , x1 , ..., xn , es posible representar pn en la forma pn (x) = y0 L0 (x) + y1 L1 (x) + ... + yn Ln (x) donde los polinomios Lk son de grado n. Queremos que para cada i se cumpla la condicin pn (xi ) = yi , es decir: pn (xi ) = y0 L0 (xi ) + y1 L1 (xi ) + ... + yn Ln (xi ) = yi i = 0, 1, ..., n

Una forma simple de elegir Lk consiste en exigir la condicin Lk (xj ) = 0 si k 6= j y Lk (xk ) = 1, donde j, k = 0, 1, ..., n, esto se consigue deniendo Lk (x) = (x x0 )(x x1 )...(x xk1 )(x xk+1 )...(x xn ) (xk x0 )(xk x1 )...(xk xk1 )(xk xk+1 )...(xk xn )

6.1. INTERPOLACIN

81

Como el numerador de Lk es el producto de n factores (x xi ), i = 0, 1, ..., n, i 6= k, entonces el polinomio Lk es de grado n. Con esto, el polinomio pn tiene grado menor o igual que n, como requerimos. Ms an, para xi , i = 0, 1, ..., n, tenemos pn (xi ) = yi Li (xi ) = yi y vemos que pn realmente interpola f en tales puntos. Por lo tanto, la forma de Lagrange para el polinomio interpolador es:
n X yk Lk (x) pn (x) = k=0

donde Lk (x) =

i=0 (x xi ) i6=k Yn i=0 (xk xi ) i6=k

Yn

Ejemplo 6.2 La siguiente tabla resume los datos de alguna observacin xi 2 0 3 5 yi = f (xi ) 3 1 4 9 Queremos interpolar f usando un polinomio interpolador segn la forma de Lagrange. As, el polinomio buscado es p3 (x) = y0 L0 (x) + y1 L1 (x) + y2 L2 (x) + y3 L3 (x) Como y0 , ..., y3 son conocidos, resta slo conocer los polinomios L0 , ..., L3 : L0 = L1 = L2 = L3 = (x x1 )(x x2 )(x x3 ) (x0 x1 )(x0 x2 )(x0 x3 ) (x x0 )(x x2 )(x x3 ) (x1 x0 )(x1 x2 )(x1 x3 ) (x x0 )(x x1 )(x x3 ) (x2 x0 )(x2 x1 )(x2 x3 ) (x x0 )(x x1 )(x x2 ) (x3 x0 )(x3 x1 )(x3 x2 ) = = = = (x 0)(x 3)(x 5) x3 8x2 + 15x = (2 0)(2 3)(2 5) 70 3 2 (x (2))(x 3)(x 5) x 6x x + 30 = (0 (2))(0 3)(0 5) 30 3 (x (2))(x 0)(x 5) x 3x2 10x = (3 (2))(3 0)(3 5) 30 3 2 (x 0)(x 3)(x 5) x x 6x = (5 (2))(5 0)(5 3) 70

Por lo tanto, 3 3 3 x 8x2 + 15x x 6x2 x + 30 x 3x2 10x p3 (x) = (3) +1 +4 70 30 30 3 2 x x 6x +9 70 1 (15x3 57x2 + 246x + 210) = 210 Observe la gura 5.3.

82

CAPTULO 6. INTERPOLACIN Y APROXIMACIN

Figura 5.3: Polinomio p3 interpolando 4 puntos, segn Lagrange

6.2.

Aproximacin: Caso Discreto

En la seccin anterior vimos algunos mtodos para interpolar una funcin f en los puntos x0, x1 , ..., xn , por una funcin g previamente denida. Pero interpolar no es aconsejable cuando: 1. Es necesario obtener un valor aproximado de la funcin en algn punto fuera del intervalo de la tabla de denicin, a esto se conoce con el nombre de extrapolacin. 2. Los valores en la tabla de datos son resultados de algn experimento fsico, los cuales pueden contener errores inherentes, que en general no son previsibles. Surge la necesidad de buscar una funcin que sea una buena aproximacin para los valores deseados y, que en cierto modo, resuelva nuestro problema con un cierto margen de seguridad. As, sean ahora los puntos (x1 , f (x1 )), ..., (xm , f (xm )), donde x1 , ..., xm [a, b] . El problema que intentamos resolver consiste en lo siguiente: Una vez elegida n funciones g1 , ..., gn : R 7 R continuas en [a, b] , obtener n constantes 1 , ..., n tales que la funcin (x) = 1 g1 (x) + ... + n gn (x) (6.3)

se aproxime lo mejor posible a f. Las funciones gi , i = 1, ..., n, pueden asumir formas no lineales, tales como g1 (x) = ex , g2 (x) = sen x, g3 (x) = x2 + log x, etc. Pero, Con qu criterio elegir las funciones gi , i = 1, ..., n? consiste en lo siguiente: El procedimiento comn

1. Gracar los puntos (x1 , f (x1 )), ..., (xm , f (xm )) en el plano cartesiano. El grco resultante se denomina diagrama de dispersin. Por medio de este diagrama es posible identicar la curva, o una combinacin lineal de curvas, que mejor se aproxime a nuestro problema. 2. Basndose en fundamentos tericos del experimento que nos otorgaron las observaciones.

6.2. APROXIMACIN: CASO DISCRETO Ejemplo 6.3 (Eleccin de gi basado en el diagrama de dispersin) Considere los siguientes datos, resultados de alguna observacin: x 2 1, 5 0, 5 0 1 1, 2 2, 5 2, 65 2, 8 3, 1 f (x) 3 1 1 1, 5 2 4 6 7 10 12

83

(6.4)

El diagrama de dispersin est dado en la gura 5.4. Claramente, una forma conveniente es elegir g1 (x) = x2 y g2 (x) = 1, entonces encontrar los valores de 1 y 2 que determinan la funcin (x) = 1 x2 + 2 .

Figura 5.4: Diagrama de dispersin Ejemplo 6.4 (Eleccin de gi utilizando fundamentos tericos) Si consideramos una experiencia donde fueron medidos varios valores de corriente elctrica que pasan por una resistencia sometida a varias tensiones, podemos elegir gi adecuadamente, colocando esos valores inicialmente en el grco 5.5:

Figura 5.5: Eleccin de gi mediante fundamento terico En este caso existe un fundamento terico relacionado a la corriente elctrica con la tensin: V = RI . Es decir, V (voltaje) es una funcin lineal de I (intensidad). Luego, como busca aproximar la tensin relacionada a esos datos, resulta prudente elegir gi (I ) = I, de tal modo que: (I ) = 1 g1 (I ) = 1 I

84

CAPTULO 6. INTERPOLACIN Y APROXIMACIN

Note que existe una sola incgnita, 1 . Al resolver este problema, el valor que tome 1 ser a su vez una aproximacin a la resistencia R. En resumen, al determinar 1 habremos ajustado ese diagrama a una curva lineal, todo eso basado en un fundamento terico. Una vez que en (6.3) fueron elegidas las funciones gi , i = 1, 2, ..., n, resta calcular las constantes que denen (x) = 1 g1 (x) + ... + n gn (x). A continuacin veremos una de las formas ms comunes.

6.2.1.

Mtodo de Mnimos cuadrados

Dados los puntos conocidos (x1 , f (x1 )), ..., (xm , f (xm )), donde las funciones g1 , ..., gn : R 7 R fueron elegidas de la forma adecuada. En experimentos de este tipo, por lo general, se tiene que m n (m mucho mayor que n). Nuestro objetivo es encontrar los valores de 1 , ..., n de modo que (x) = 1 g1 (x) + ... + n gn (x) se aproxime lo mejor posible a la supuesta funcin f. Anlogamente a lo que vimos cuando estudiamos interpolacin, desearamos ahora que: (x1 ) = 1 g1 (x1 ) + ... + n gn (x1 ) = f (x1 ) . . . (xm ) = 1 g1 (xm ) + ... + n gn (xm ) = f (xm )

(6.5)

Este sistema de m las y 1 , ..., n variables, podemos nuevamente verlo matricialmente as: A = b donde g1 (x1 ) . ... . A= . g1 (xm ) gn (x1 ) . . , . gn (xm ) 1 . = . . n f (x1 ) . . y b= . f (xm ) (6.6)

Note que el sistema (6.6) por lo general es incompatible, es decir, que no tiene solucin. Esto sucede justamente porque m (el nmero de observaciones) es mucho mayor que n (las incgnitas). Gustaramos ahora de encontrar Rn , tal que A b ( A prximo de b). El mtodo de Mnimos Cuadrados consiste en calcular del siguiente modo: = m n kA bk2 n
R

(6.7)

donde kk es la norma eucldea. Observe que si i representa la i-la de A, entonces ai bi representa (xi ) f (xi ), esta diferencia es denida como el i-desvio: di = (xi ) f (xi ) Adems,
m m X X 2 (ai bi ) = d2 kA bk = i 2 i=1 i=1

6.2. APROXIMACIN: CASO DISCRETO

85

P 2 Es decir, al minimizar kA bk2 estamos minimizando m i=1 di , la sumatoria de los desvos al cuadrado, de ah el nombre del mtodo (vea la gura 5.6)

Figura 5.6: Valor absoluto de la desviacin di Para minimizar la funcin q : Rn 7 R, denida por la regla q () = kA bk2 , utilizamos algunos resultados del clculo diferencial. Despus de derivar parcialmente q , igualando a cero, tenemos el siguiente sistema de n incgnitas y n ecuaciones: q () = 0 1 . . . q () = 0 n Esto es equivalente resolver q () = 0, donde q() = At A At b As, si A tiene rango completo, entonces At A es inversible. Luego, existe una nica solucin para el sistema q () = 0, la cual la obtenemos de resolver: At A = At b o = (At A)1 At b (6.9) (6.8)

Note que At A es una matriz cuadrada e inversible desde que A tenga rango completo. En estas condiciones, la matriz At A es simtrica y denida positiva, as, el sistema (6.8) puede resolverse utilizando la descomposicin de Cholesky. Ejemplo 6.5 Considere los siguientes datos vistos en (6.4), cuyo diagrama de dispersin asociado est en la gura 5.4. Despus de elegir g1 (x) = x2 y g2 = 1, tal que (x) = 1 x2 + 2 ,

86

CAPTULO 6. INTERPOLACIN Y APROXIMACIN

procedemos a encontrar 1 y 2 exigiendo que:

(2) (1, 5) (0, 5) (0) (1) (1, 2) (2, 5) (2, 65) (2, 8) (3, 1)

= = = = = = = = = =

1 (2)2 + 2 = 3 1 (1, 5)2 + 2 = 1 1 (0, 5)2 + 2 = 1 1 (0)2 + 2 = 1, 5 1 (1)2 + 2 = 2 1 (1, 2)2 + 2 = 4 1 (2, 5)2 + 2 = 6 1 (2, 65)2 + 2 = 7 1 (2, 8)2 + 2 = 10 1 (3, 1)2 + 2 = 12

y obteniendo el sistema lineal A = b. donde: y = 1 2

A=

4, 0000 2, 2500 0, 2500 0, 0000 1, 0000 1, 4400 6, 2500 7, 0225 7, 8400 9, 6100

1 1 1 1 1 1 1 1 1 1

b=

3, 0000 1, 0000 1, 0000 1, 5000 2, 0000 4, 0000 6, 0000 7, 0000 10, 0000 12, 0000

Resolviendo por mnimos cuadrados, es decir, resolviendo At A = At b : tenemos 1 2 = 1, 22449857738665 0, 606667482559783

266, 3943 39, 6625 39, 6625 10

1 2

302, 1375 42, 5000

La curva que se aproxima, de la mejor manera, es (x) = 1, 224498x2 0, 606667, el grco

6.2. APROXIMACIN: CASO DISCRETO est dado en la gura 5.7.

87

Figura 5.7: aproximando el diagrama de dispersin

88

CAPTULO 6. INTERPOLACIN Y APROXIMACIN

Captulo 7 Integracin numrica

Sabemos que si una funcin f : R 7 R es continua en el intervalo [a, b] , entonces tal funcin tiene antiderivada F en [a, b] , es decir F 0 (x) = f (x), y Z x [a, b]

f (x)dx = F (b) F (a)

Sin lugar a dudas, esto constituy uno de los ms grandes descubrimientos dentro del clculo innitesimal. Pero es casi improbable calcular la integral de una funcin arbitraria f . En este captulo veremos algunos de los ms populares mtodos numricos para calcular integrales denidas. La idea bsica de los mtodos numricos consiste en establecer una particin uniforme {x0 , x1 , ..., xn1 , xn } del intervalo de integracin [a, b] : a = x0 , x1 , ..., xn1 , xn = b (7.1)

donde x0 < x1 < xn y h = x1 x0 = x2 x1 = = xn xn1 = (b a)/n. Posteriormente, sustituir la funcin f por polinomios que la aproximen razonablemente en el intervalo [a, b] . As, el problema se resumira a la integracin de polinomios, lo cual es trivial.

7.1.

Mtodo de los Trapecios

Si usamos la frmula de Lagrange para aproximar f por un polinomio p1 de grado uno en x0 y x1 , tenemos Z x1 Z x1 Z x1 x x1 x x0 f (x0 ) + f (x1 ) dx f (x)dx p1 (x)dx = h h x0 x0 x0 h = (f (x0 ) + f (x1 )) = T 2 Observe que, cuando f sobre [x0 , x1 ] est por encima del eje de las abscisas, T puede ser interpretada como el rea del trapecio mostrado en la gura 6.1, la cual es una aproximacin 89

90

CAPTULO 7. INTEGRACIN NUMRICA

para el valor de la integral de f en [x0 , x1 ] . Mientras ms pequeo sea h, mejor ser la aproximacin.

Figura 6.1: Podemos repetir ese procedimiento para cada uno de los subintervalos generados por la particin uniforme, as Z b n1 Z xk+1 n1 n1 X X hX h (f (xk ) + f (xk+1 )) = f (x)dx = f (x)dx (f (xk ) + f (xk+1 )) 2 2 a x k k=0 k=0 k=0

Note adems que xk = x0 + kh = a + kb, esto nos lleva al siguiente algoritmo.

Algoritmo 7.1 (Regla de los Trapecios) Dada la funcin f : R 7 R continua en [a, b] . Paso inicial: Denir n y calcular h = (b a)/n Paso principal: Calcular hX I= {f (a + kh) + f (a + (k + 1)h)} 2 k=0
n1

I es una aproximacin a

Una funcin en Matlab quedara as: function I=trapecio(a,b,n) h=(b-a)/n; I=0; for k=1:n I=I+f(a+(k-1)*h)+f(a+k*h); end I=I*h/2 Ejemplo 7.1 Usando la funcin de arriba, calcule y verique que: R2 1. Para n = 10000, se tiene 1 x2 dx 3, 00000004499998 R2 2 2. Para n = 10000, se tiene 0 ex dx 0, 882081390518215 R7 3. Para n = 10000, se tiene 3 sen(7x3 + sen(ln(x4 + x2 + 5)))dx 0, 681133991768069

Rb
a

f (x)dx

7.2. MTODO DE SIMPSON

91

7.2.

Mtodo de Simpson

Consideremos la misma particin introducida en (7.1). Podemos usar de nuevo la frmula de Lagrange para obtener la frmula de integracin resultante de la interpolacin de f por un polinomio de grado 2. Sea p2 el polinomio que interpola f en los puntos x0 , x1 = x0 + h y x2 = x0 + 2h (vea la gura 6.2)

Figura 6.2: As, p2 (x) = Luego, Z (x x1 )(x x2 ) (x x0 )(x x2 ) (x x0 )(x x1 ) f (x0 ) + f (x1 ) + f (x2 ) (h)(2h) (h)(h) (2h)(h) Z
x2

x2

x0

Z f (x0 ) x2 f (x)dx p2 (x)dx = (x x1 )(x x2 )dx 2h2 x0 x0 Z Z f (x1 ) x2 f (x2 ) x2 2 (x x0 )(x x2 )dx + (x x0 )(x x1 )dx h 2h2 x0 x0 Z
x2

integrando, tenemos f (x)dx h (f (x0 ) + 4f (x1 ) + f (x2 )) = S 3

x0

Podemos tomar de tres puntos de la particin y aplicar repetidamente esta regla, para eso requerimos que la particin sea denida por n par. Luego: Z
b

f (x)dx =

xn

f (x)dx =

x0

+ [f (x2 ) + 4f (x3 ) + f (x4 )] + + [f (xn2 ) + 4f (xn1 ) + f (xn )]} Ms an, Z

b

n/2 Z X k=1

x2k

x2k2

f (x)dx

h {[f (x0 ) + 4f (x1 ) + f (x2 )] 3

f (x)dx

h {[f (x0 ) + f (xn )] + 4 [f (x1 ) + f (x3 ) + + f (xn1 )] 3 +2 [f (x2 ) + f (x4 ) + + f (xn2 )]}

Esto nos lleva al siguiente algoritmo.

92

CAPTULO 7. INTEGRACIN NUMRICA

Algoritmo 7.2 (Regla de Sipmson) Dada la funcin f : R 7 R continua en [a, b]. Paso inicial: Denir n y calcular h = (b a)/n Paso principal: Calcular S = h {[f (a) + f (b)] + 4 [f (a + h) + f (a + 3h) + + f (a + (n 1)h)] 3 +2 [f (a + 2h) + f (a + 4h) + + f (a + (n 2)h)]} Rb
a

S es una aproximacin para

f (x)dx.

Una funcin en Matlab para este caso, estara dada del siguiente modo: function S=simpson(a,b,n) S1=0; S2=0; h=(b-a)/n; for k=1:2:n-1 S1=S1+f(a+k*h); end for k=2:2:n-2 S2=S2+f(a+k*h); end S=(h/3)*((f(a)+f(b))+4*S1+2*S2); Ejemplo 7.2 Implementamos el mtodo de Simpson y lo usamos para calcular la misma inR2 2 tegral 0 e x dx, del ejemplo anterior. Para el caso de apenas n = 100, obtuvimos Z
2 0

ex dx 0, 882081390111187

Observacin 7.1 En la prctica, el mtodo de Simpson es ms eciente que el mtodo del trapecio, una manera de medir analticamente esto es por medio del estudio de error cometido. Problema 7.1 Haga un estudio de la estimativa de los errores cometidos por los mtodos del Trapecio y Simpson.

Bibliografa
[1] Steven C. ,Mtodos Numricos para ingenieros .quinta edicin. McGraw-Hill, 2006 [2] Jaan K. Numerical Methods in Engineering with MATLAB. 2005 [3] John H. Mtodos numricos con Matlab. tercera edicin. Prentice 2000 [4] Hern M. Matlab 7.para ciencias e ingeniera. primera edicin. Megabyte. 2005 [5] Burden R. Anlisis numrico. 7 edicin. Cengage

93