Descarga la aplicación para disfrutar aún más
Vista previa del material en texto
Problemas de Métodos Numéricos Miguel Alemán Flores, Luis Alvarez León y Javier Sánchez Pérez Departamento de Informática y Sistemas Universidad de Las Palmas Campus de Tafira 35017 Las Palmas, España Tfl: 45.87.10/08 Email: {maleman/lalvarez/jsanchez}@dis.ulpgc.es Contents 1 INTRODUCCION. 1 2 ARITMETICAS DE PRECISION FINITA Y FUENTES DE ERRORES NUMERI- COS. 1 3 CALCULO DE CEROS DE UNA FUN- CION 4 4 INTERPOLACION DE FUNCIONES I 6 5 ANALISIS NUMERICO MATRICIAL I 9 6 DIFERENCIACION E INTEGRACION NUMERICA 14 7 ANALISIS NUMERICO MATRICIAL II 23 8 INTERPOLACION DE FUNCIONES II 36 INTRODUCCION. El presente documento es el libro de problemas donde se encuentran resueltos todos los problemas presentes en el libro de Métodos Numéricos publicado por los mismos au- tores. Nunca se insistirá lo suficiente sobre la necesidad de hacer problemas para comprender correctamente cualquier teoría y sus aplicaciones. Además la manera de afrontar el estudio de los problemas debe ser bien distinta a la forma de estudiar teoría. Primero se debe intentar hacer los problemas sin mirar en absoluto la solución y después de reflexionar e intentar resolverlo de diferentes formas, muchas de las cuales nos llevarán a callejones sin salida, se mirará la solución. Es un hecho fácilmente constatable, que se aprende mucho más de un problema que no se ha conseguido resolver, pero al que se ha dedicado suficiente esfuerzo, que de un problema del cual se mira directamente la solución sin ninguna fase de reflexión previa. Además se tiende a olvidar con facilidad la técnica de resolución de un problema sobre el cual no se ha reflexionado sufi- cientemente. De todo ello se deduce que el estudio cor- recto de los problemas de una asignatura va reñido con las prisas de última hora que suelen asaltar a los estudiantes cuando se acercan los exámenes, puesto que el esfuerzo de reflexión que requieren precisa de un trabajo diario y continuado, difícilmente compatible con las prisas de úl- tima hora. Resulta inquietante observar como en muchas ocasiones la realización de problemas se aborda bajo un espíritu de aprender rápidamente 4 técnicas básicas, que muchas veces ni se entienden, y a partir de ahí intentar re- producir esas técnicas, de forma absolutamente mecánica, en problemas análogos. El problema de esta actitud, es que aunque a corto plazo puede dar lugar a resultados posi- tivos, aprobando asignaturas con un conocimiento mínimo e insuficiente, a la larga, tiene efectos catastróficos sobre la formación del alumno, a través de una disminución im- portante de la capacidad de razonamiento y del sentido crítico. ARITMETICAS DE PRECISION FINITA Y FUENTES DE ERRORES NUMERICOS. Problema 1 Demostrar que al representar el número real 0.1 como 0.1 = 2e ∞X n=1 an 2n el número de elementos no-nulos an es infinito. Solución: Supongamos que para algún t finito y e entero se tiene: 0.1 = 2e tX n=1 an 2n despejando en esta igualdad obtenemos 2t−e = 10 tX n=1 an2 t−n ahora bien, como el número m = Pt n=1 an2 t−n es entero, de la desigualdad anterior obtenemos 2t−e = 5 · 2m pero esta igualdad implica que el número 2t−e es divisible por 5 lo cual es imposible. 1 Problema 2 Representar el número 0.0 703 125 como 0.0 703 125 = 2e ∞X n=1 an 2n Solución: En primer lugar tenemos que encontrar un en- tero e tal que 1 2 ≤ 0.0 703 125 · 2−e < 1 para e = −3 obtenemos 0.0 703 125 · 23 = 0. 562 5 ahora tenemos que escribir el número 0.5625 como 0.5625 = 1 2 + ∞X n=2 an 2n los an se calculan de la siguiente forma 0.5625 < 1 2 + 1 22 = 0.75⇒ a2 = 0 0.5625 < 1 2 + 1 23 = 0.625⇒ a3 = 0 0.5625 = 1 2 + 1 24 = 0.5625⇒ a4 = 1 por tanto 0.0 703 125 = 2−3 µ 1 2 + 1 24 ¶ en términos binarios, este numero se escribiría con e = −3 y la mantisa viene dada por la secuencia 1, 0, 0, 1, 0, 0, .... (si no almacenamos el primer término a1 porque siempre es 1, la mantisa sería 0, 0, 1, 0, 0, ....) Problema 3 (1 puntos) Calcular los valores positivos mínimo y máximo que puede tomar un número real en una aritmética de precisión finita en función de t, emin y emax. Solución: Los valores positivos mínimo y máximo son xmin = 2 emin−1 xmax = 2 emax tX n=1 1 2n = 2emax 1 2 − 1 2t+1 1 2 = 2emax µ 1− 1 2t ¶ Problema 4 Calcular todos los números reales que se pueden construir tomando 5 bits de la forma siguiente: 1 bit para el signo, 2 bits para la mantisa (es decir t = 3, puesto que a1 = 1 sólo se almacenan a2 y a3, y 2 bits para el exponente e, tomando como rango de e = −1, 0, 1, 2. Representar dichos números sobre una recta. Solución: Los valores posibles positivos se representan en la siguiente tabla e = −1 1 22 , 1 22 + 1 24 , 1 22 + 1 23 , 1 22 + 1 23 + 1 24 e = 0 1 2 , 1 2 + 1 23 , 1 2 + 1 22 , 1 2 + 1 22 + 1 23 e = 1 1, 1 + 1 22 , 1 + 1 2 , 1 + 1 2 + 1 22 e = 2 2, 2 + 1 2 , 2 + 1, 2 + 1 + 1 2 los valores negativos son los mismos cambiados de signo. Simplificando las fracciones nos queda e = −1 0.25, 0.3125, 0.375, 0.437 5 e = 0 0.5, 0.625, 0.75, 0.875 e = 1 1, 1.25, 1.5, 1.75 e = 2 2, 2.5, 3, 3.5 Si representamos los números positivos sobre una recta obtenemos 0 1 2 3 4 -2 -1 0 1 2 x y Problema 5 Dada una aritmética de precisión finita cualquiera, calcular la distancia que hay entre el número 1 y su inmediato superior (es decir el número que va después de 1), y la distancia entre el número 1, y su inmediato in- ferior. Solución: El número 1 en una aritmética de precisión finita se escribe como 1 = 2 µ 1 2 ¶ el número inmediato superior a 1 en la aritmética es 2 µ 1 2 + 1 2t ¶ = 1 + 1 2t−1 y el número inmediato inferior a 1 viene dado por 1 2 + ..+ 1 2t = 1 2 − 1 2t+1 1− 12 = 1− 1 2t 2 Problema 6 Se considera una aritmética de 16 bits donde se dedican 1 bit al signo, 9 bits a la mantisa (t = 10) y 6 bits al exponente ( emin = −30 emax = 31). Escribir, si es posible, los siguientes números en esta aritmética: 1. 2, y los números más cercanos a 2 por arriba y por debajo. Solución: 2 = 22 µ 1 2 ¶ Si guiente = 22 µ 1 2 + 1 210 ¶ Anterior = 2 à 10X i=1 1 2i ! = 2 µ 1 2 − 1 211 1 2 ¶ 2. El cero, el infinito y Na. Solución: 0 = 2−31 µ 1 2 ¶ ∞ = 232 µ 1 2 ¶ NaN = 232 µ 1 2 + 1 22 ¶ 3. Los números positivos más grande y más pequeño de la aritmética (teniendo en cuenta las excepciones) Solución: Mayor = 231 à 10X i=1 1 2i ! = 231 µ 1 2 − 1 211 1 2 ¶ Menor = 2−31 µ 1 210 ¶ 4. 19 . Solución: No se puede escribir de forma exacta. Si suponemos 1 9 = 2e à tX i=1 ai 2i ! =⇒ 1 = 9 · 2e à tX i=1 ai 2i ! =⇒ 2t−e = 32 à tX i=1 ai2 t−i ! =⇒ 2t−e = 32m donde m es un número entero. Ahora bien esta igual- dad es imposible porque resultaría que 3 divide a 2. 5. 2 ¡ 1 2 − 1 210 ¢ . Solución: 2 ¡ 1 2 − 1 210 ¢ = 20 ¡ 1 2 + 1 22 + 1 23 + 1 24 + 1 25 + 1 26 + 1 27 + 1 28 + 1 29 ¢ Problema 7 Sean A = 2 ¡ 1 2 + 1 23 + 1 25 ¢ B = 23 ¡ 1 2 + 1 26 + 1 27 ¢ . Calcular B +A y B −A Solución: B +A = 23 ¡ 1 2 + 1 23 + 1 24 ¢ 1. B −A = 22 ¡ 1 2 + 1 23 + 1 24 + 1 25 ¢ Problema 8 Sean emin, emax, los valores mínimo y máx- imo del exponente e. Demostrar que si emin < e < emax, entonces los números: 2e à tX n=1 an 2n ± 1 2t ! pertenecen al conjunto A de números reales generados por la aritmética de precisión finita. Solución: Que los números pertenecen a la aritmética significa que existe un conjunto de valores binarios a0i y un entero e0 tal que 2e à tX n=1 an 2n ± 1 2t ! = 2e 0 tX n=1 a0n 2n Consideremos primero el caso de sumar 1/2t. Si ak = 1 para todo k, entonces 2e à tX n=1 1 2n + 1 2t ! = 2e+1 1 2 Si por el contrario existe un k0 tal que ak0 = 0, y tal que ak = 1 para todo k0 < k ≤ t entonces basta tomar a0k = ak si 1 ≤ k < k0, a0k0 = 1 y a 0 k = 0 si k0 < k ≤ t Consideremos ahora el caso de restar 1/2t. Si el único elemento ak distinto de 0 es a1, entonces 2e µ 1 2 − 1 2t ¶ = 2e−1 tX n=1 1 2n Si por el contrario existe un k0 > 1 tal que ak0 = 1, y tal queak = 0 para todo k0 < k ≤ t entonces basta tomar a0k = ak si 1 ≤ k < k0, a0k0 = 0 y a 0 k = 1 si k0 < k ≤ t. Problema 9 Dado un número ez = 2ePtn=1 an2n , en una aritmética de precisión finita. Calcular el número inmedi- atamente inferior y superior a él en dicha aritmética. Solución: Si el número es de la forma ez = 2e 1 2 entonces el inmediato superior es ez + 2e 1 2t y el inmediato inferior es 2e−1 tX n=1 1 2n 3 para cualquier otro número ez, el inmediato superior e in- ferior son ez ± 2e 1 2t Problema 10 (1 puntos) Calcular las raíces del poli- nomio P (x) = x2 − 2x+ 0.01 evitando los errores de can- celación. Solución: x1 = 2 + √ 4− 0.04 2 = 1.995 x2 = 0.01 1.995 Problema 11 Escribir el pseudocódigo para implementar el cálculo de las raíces reales de ax2+bx+c = 0 evitando los errores de cancelación y teniendo en cuenta las diferentes opciones que aparecen cuando a 6= 0 y a = 0. Solución: Algoritmo Calculo raiz polinomio ax2 + bx+ c = 0 variables reales a,b,c leer(a,b,c) si (a==0 ) entonces si (b==0 ) entonces PRINT ’EL POLINOMIO ES CONSTANTE’ parar finsi PRINT ’EL POLINOMIO ES DE GRADO 1.’ PRINT ’LA RAIZ ES ’,−c/b parar finsi d=b*b-4*a*c si (d< 0 ) entonces PRINT EL POLINOMIO NO TIENE RAICES REALES parar finsi si (b> 0) entonces x1=(-b-SQRT(d))/(2*a) además x1=(-b+SQRT(d)/(2*a) finsi x2=c/(x1*a) PRINT *,’LAS RAICES SON: ’,x1,x2 fin algoritmo CALCULO DE CEROS DE UNA FUNCION Problema 12 Calcular 2 iteraciones del algoritmo de la bisección para buscar un cero de la función f(x) = x2 − 2 en el intervalo [−2, 0] Solución: x = 0 + (−2) 2 = −1 f(−2) > 0, f(0) < 0, f(−1) < 0 Nuevo Intervalo = [−2,−1] x = −1 + (−2) 2 = −1.5 f(−2) > 0, f(−1) < 0, f(−1.5) > 0 Nuevo Intervalo = [−1.5,−1] Problema 13 Escribir el pseudocódigo del algoritmo el método de la bisección Solución: Algoritmo: Método de la bisección variables reales x,a,b,tol leer(a,b,tol) si (a> b) entonces PRINT ’INTERVALO INCORRECTO’ parar finsi si (F(a)*F(b)> 0) entonces PRINT ’NO HAY CAMBIO DE SIGNO EN EL INTERVALO’ parar finsi mientras (F(x)!=0 Y (b-a)>tol) x=(a+b)/2 si((F(a)*F(x))<0) entonces b=x además A=X finsi fin mientras _PRINT ’LA RAIZ ES’ x fin algoritmo real F(real x) real a←− cos(x) + x ∗ x− 6 devolver a fin función Problema 14 Calcular 2 iteraciones del algoritmo de la regula-falsi para buscar un cero de la función f(x) = x2−2 en el intervalo [0, 2] 4 Solución: x = 0− 2 f(2)− f(0)f(0) = 1 f(2) > 0, f(0) < 0, f(1) < 0 Nuevo Intervalo = [1, 2] x = 1− 1 f(2)− f(1)f(1) = 4 3 f(2) > 0, f(1) < 0, f( 4 3 ) < 0 Nuevo Intervalo = [ 4 3 , 2] Problema 15 Escribir el pseudocódigo del algoritmo del método de la Regula-falsi Solución: Algoritmo: Método de la Regula-falsi variables reales x,a,b,tol leer(a,b,tol) si (a> b) entonces PRINT ’INTERVALO INCORRECTO’ parar finsi si (F(a)*F(b)> 0) entonces PRINT ’NO HAY CAMBIO DE SIGNO EN EL INTERVALO’ parar finsi mientras (F(x)!=0 Y (b-a)>tol) x=a-F(a)*(b-a)/(F(b)-F(a)) si((F(a)*F(x))<0) entonces b=x además a=x finsi fin mientras _PRINT ’LA RAIZ ES’ x fin algoritmo Problema 16 Calcular una iteración del método de Newton-Raphson para calcular un cero de la función f(x) = x3 − 3 partiendo de x0 = 1. Solución: x1 = 1− −2 3 = 5 3 Problema 17 (1 punto Calcular una iteración del método de la secante para calcular un cero de la función f(x) = x3 − 3 partiendo de x0 = 0, x1 = 1 Solución: x1 = 1− −2³ −2−(−3) 1−0 ´ = 3 Problema 18 Escribir pseudocódigo del algoritmo del método de la Secante utilizando reales de doble precisión. Los datos de entrada son las aproximaciones iniciales x0, y x1, El número máximo de iteraciones N max, y la toler- ancia TOL para determinar la igualdad de dos números. Solución: Algoritmo: Método de la secante variables reales x0,x1,x2,tol variable entera Nmax leer(a,b,tol,Nmax) si (x0==x1) entonces PRINT ’LAS DOS APROXIMACIONES INI- CIALES COINCIDEN’ parar finsi para k←− 1 hasta Nmax hacer si(ABS(x1-x0)< tol) entonces PRINT *,’LA RAIZ DE LA FUNCION ES: ’,x1 parar finsi si(F(x1)==F(x0)) entonces PRINT *,’METODO NO CONVERGE’ parar finsi x2=x1-F(x1)*(x1-x0)/(F(x1)-F(x0)) x0=x1 x1=x2 fin para PRINT *,’NUMERO MAXIMO DE ITERACIONES EXCEDIDO’ fin algoritmo Problema 19 Calcular una iteración del método de Muller para calcular un cero de la función f(x) = x3 − 3 partiendo de x0 = 1 (Calculando las derivadas de la fun- ción de forma exacta) y quedándonos con la raíz más cer- cana a x0. Solución: −2 + 3(x− 1) + 3(x− 1)2 = 0 x1 = 1 + −3 + √ 33 6 Problema 20 Dado el polinomio P (x) = 2x3+3x2+4x+ 5. Evaluar el polinomio y su derivada en el punto x = 2, utilizando el algoritmo de Horner 5 Solución: P (x) = ((2x+ 3)x+ 4)x+ 5 P (2) = ((7)2 + 4)2 + 5 P (2) = (18)2 + 5 = 41 P 0(x) = (2x+ 7)x+ 18 P 0(2) = (4 + 7)2 + 18 = 40 Problema 21 Calcular el número máximo de raíces pos- itivas y negativas del polinomio x5− 35x3+30x2+124x− 120, y localizarlas en un intervalo. Solución: Teniendo en cuenta que 1 + maxk=0,..,n−1 | ak | | an | = 125 las raíces del polinomio están en el intervalo [−125, 125]. Para calcular el número máximo de raíces positivas mi- ramos los cambios de signo de los coeficientes, en este caso los signos son: +−++− por tanto el número de raíces positivas es 1 ó 3. Para esti- mar el número de raíces negativas cambiamos x por −x y miramos los signos de los coeficientes que en este caso son: −++−− por tanto el número de raíces negativas son 0 ó 2. Problema 22 Aislar en intervalos las raíces del poli- nomio P (x) = 20x3 − 45x2 + 30x− 1. Solución: Teniendo en cuenta que en este caso 1 + maxk=0,..,n−1 | ak | | an | = 1 + 45 20 = 65 20 todas las raíces están en el intervalo [−6520 , 65 20 ]. Para ais- lar las raíces calculamos los ceros de la derivada P 0(x) = 60x2 − 90x + 30, dichas raíces son 1 y 1/2. Por otro lado tenemos P (−65 20 ) = −1260. 4 P ( 1 2 ) = 21 4 P (1) = 4 P ( 65 20 ) = 307. 75 por tanto hay una única raíz en el intervalo [−6520 , 1 2 ]. Problema 23 Aislar en intervalos las raíces del poli- nomio P (x) = 2x3 + 3x2 − 12x+ 1 Solución: P 0(x) = 6x2 + 6x− 12 raı́ces x = 1,−2 Intervalo Inicial [−7, 7] P (−7) = −454 P (−2) = 21 P (1) = −6 P (7) = 750. Intervalos donde están las raíces: [-7,-2] [-2,1] [1,7] INTERPOLACION DE FUNCIONES I Problema 24 Calcular el polinomio interpolador de La- grange P3(x) de la función f(x) = sen(x) en los puntos 0, π2 , π y 3π 2 . Solución: Puesto que sen(0) = sen(π) = 0 sólo necesi- tamos los polinomios base de Lagrange centrados en π2 y 3π 2 . Pπ 2 (x) = x (x− π) ¡ x− 3π2 ¢ π 2 ¡ π 2 − π ¢ ¡ π 2 − 3π 2 ¢ P 3π 2 (x) = x (x− π) ¡ x− π2 ¢ 3π 2 ¡ 3π 2 − π ¢ ¡ 3π 2 − π 2 ¢ Por tanto el polinomio interpolador es P (x) = Pπ 2 (x)− P 3π 2 (x) Problema 25 Calcular la expresión del error interpo- lación al aproximar la función f(x) = sen(x) en el in- tervalo [0, 2π] interpolando en los puntos 0, π2 , π, 3π 2 . y aco- tarlo superiormente. Solución: El error de interpolación viene dada por la expresión: f(x)− PN(x) = sen(ξ) 4! x ³ x− π 2 ´ (x− π) µ x− 3π 2 ¶ el valor máximo del sen(ξ) es 1. Por otro lado el valor donde alcanza el máximo el polinomio del error en [0, 2π] es x = 2π, por tanto la cota del error que obtenemos es |f(x)− PN (x)| ≤ 1 4! 2π ³ 2π − π 2 ´ (2π − π) µ 2π − 3π 2 ¶ Problema 26 Calcular el error máximo de interpolación en el intervalo [0, 1] al interpolar la función cos(x) en los puntos dados por los polinomios de Chebyshev tomando N = 5. 6 Solución: Según las fórmulas vistas en teoría el error viene dado por la expresión: | f(x)− PN (x) |≤ maxx∈[a,b] ¯̄ fN+1)(ξ) ¯̄ (N + 1)!2N µ b− a 2 ¶N+1 en nuestro caso como N = 5 y la derivada sexta de cos(x) es − cos(x) cuyo máximo en valor absoluto es 1, obtenemos | f(x)− PN (x) |≤ 1 6!25 µ 1 2 ¶6 = 6. 78× 10−7 Problema 27 Interpolar la función f(x) = 10x2+1 en los puntos x0 = −2, x1 = −1, x2 = 1, x3 = 2 utilizando las diferencias de Newton y evaluar el polinomio en x = 0 utilizando el algoritmo de Horner. Solución: −2→2 3 -1 → 5 -1 0 0 1 → 5 -1 -3 2 → 2 P (x) = 2 + 3(x+ 2)− 1(x+ 2)(x+ 1) + 0(x+ 2)(x+ 1)(x− 1) = (−1(x+ 1) + 3)(x+ 2) + 2 P (0) = (−1(0 + 1) + 3)(0 + 2) + 2 = 6 Nota: Quitar paréntesis en P(x) y aplicar Horner sobre el polinomio resultante no es lo que pide el problema y por lo tanto está mal Problema 28 Calcular el polinomio interpolador de La- grange P3(x) de la función f(x) = sen(x) en los puntos 0, π2 , π y 3π 2 utilizando las diferencias divididas de Newton. Solución: Las diferencias divididas son: f [0] = 0, f [π2 ] = 1, f [π] = 0, f [3π2 ] = −1, f [0, π 2 ] = 2 π f [ π 2 , π] = − 2 π f [π, 3π 2 ] = − 2 π f [0, π 2 , π] = − 4 π2 f [ π 2 , π, 3π 2 ] = 0 f [0, π 2 , π, 3π 2 ] = 8 3π3 por tanto, el polinomio interpolador es P (x) = 2 π x− 4 π2 x ³ x− π 2 ´ + 8 3π3 x ³ x− π 2 ´ (x− π) Problema 29 Calcular el polinomio interpolador de La- grange P3(x) de la función f(x) = 2x en los puntos 0, 1, 3, 4 utilizando las diferencias divididas de Newton. Expresar el polinomio tomando en primer lugar x0 = 0, x1 = 1, x2 = 3 y x3 = 4, y en segundo lugar x0 = 4, x1 = 3, x2 = 1, y x3 = 0. Solución: En el primer caso, las diferencias divididas son f [x0] = 1, f [x1] = 2, f [x2] = 8, f [x3] = 16. f [x0, x1] = 1 f [x1, x2] = 3 f [x2, x3] = 8 f [x0, x1, x2] = 2 3 f [x1, x2, x3] = 5 3 f [x0, x1, x2, x3] = 1 4 y el polinomio interpolador es: P (x) = 1 + x+ 2 3 x(x− 1) + 1 4 x(x− 1)(x− 3) Si tomamos ahora los puntos en orden inverso: f [x0] = 16, f [x1] = 8, f [x2] = 2, f [x3] = 1. f [x0, x1] = 8 f [x1, x2] = 3 f [x2, x3] = 1 f [x0, x1, x2] = 5 3 f [x1, x2, x3] = 2 3 f [x0, x1, x2, x3] = 1 4 El polinomio interpolador es: P (x) = 16+8(x−4)+5 3 (x−4)(x−3)+1 4 (x−4)(x−3)(x−1) como puede observarse, al cambiar el orden de los pun- tos de interpolación, el polinomio de Lagrange expresado a través de las diferencias divididas cambia totalmente, salvo el último coeficiente 14 que es el mismo en ámbos ca- sos pues como se había demostrado en teoría el valor de f [x0, x1, x2, x3] no depende del orden en que se toman los puntos de interpolación. Problema 30 Dada una función f(x), y una secuen- cia de valores xn, aproximar f(x) por la parábola que pasa por los puntos (xn−1, f(xn−1)) , (xn−2, f(xn−2)) y (xn−3, f(xn−3)), calcular posteriormente las derivadas del polinomio, y comprobar que coinciden con las fórmu- las dadas en el método de Muller para el cálculo de las derivadas f 00(xn−1) y f 0(xn−1). 7 Solución: Si utilizamos las diferencias divididas para in- terpolar obtenemos f [xn−1] = f(xn−1) f [xn−1, xn−2] = f(xn−1)− f(xn−2) xn−1 − xn−2 f [xn−2, xn−3] = f(xn−2)− f(xn−3) xn−2 − xn−3 f [xn−1, xn−2, xn−3] = f [xn−1, xn−2]− f [xn−2, xn−3] xn−1 − xn−3 El polinomio interpolador es P (x) = f(xn−1) + f [xn−1, xn−2](x− xn−1)+ f [xn−1, xn−2, xn−3](x− xn−1)(x− xn−2) por tanto P 00(xn−1) = 2f [xn−1, xn−2, xn−3] P 0(xn−1) = f [xn−1, xn−2] + f [xn−1, xn−2, xn−3](xn−1 − xn−2) que corresponde a las fórmulas utilizadas por el método de Muller. Problema 31 Aproximar la función sen(x) en el inter- valo [0, π4 ] utilizando el desarrollo de Taylor, y calcular el valor de n a partir del cual la aproximación es la mejor posible dentro de una aritmética de 32 bits. Solución: El desarrollo de Taylor en 0 del sen(x) viene dado por: sen(x) u Pn(x) = x− x3 3! + x5 5! + ....+ (−1)n x 2n+1 (2n+ 1)! y el error máximo cometido por el desarrollo de Taylor en un punto x ∈ [0, π4 ] es | Pn(x)− sen(x) |≤ sen( π 4 ) (x) 2n+2 (2n+ 2)! donde ξ ∈ [0, π4 ]. Para que la aproximación sea la mejor dentro de una aritmética de 32 bits tiene que cumplirse | Pn(x)− sen(x) | sen(x) ≤ 2−24 = 5. 96× 10−8 por otro lado, en el intervalo [0, π4 ] se verifica sen ¡ π 4 ¢ π 4 x ≤ sen(x) por tanto: | Pn(x)− sen(x) | sen(x) ≤ ¡ π 4 ¢2n+2 (2n+ 2)! para n = 4 se tiene que¡ π 4 ¢2n+2 (2n+ 2)! = 2. 46× 10−8 por tanto n = 4 determina la mejor aproximación en una aritmética de 32 bits. Problema 32 Demostrar que utilizando relaciones trigonométricas es posible calcular las funciones sen(x) y cos(x) para cualquier x (en radianes), utilizando únicamente su valor en el intervalo [0, π8 ]. Solución: En teoría se demostró como se pueden definir el sen(x) y cos(x) para cualquier valor de x a partir de su definición en [0, π4 ], por tanto, en este problema sólo tenemos que definir las funciones trigonométricas en [0, π4 ] a partir de su definición en [0, π8 ]. Basta tener en cuenta las relaciones: cos[0,π4 ](x) = ( cos[0,π8 ](x) si x ≤ π 8 cos2[0,π8 ] ¡ x 2 ¢ − sin2[0,π8 ] ¡ x 2 ¢ si x > π8 sen[0,π4 ](x) = ½ sen[0,π8 ](x) si x ≤ π 8 2 cos[0,π8 ] ¡ x 2 ¢ sin[0,π8 ] ¡ x 2 ¢ si x > π8 Problema 33 Calcular los polinomios necesarios para in- terpolar las funciones trigonométricas cos(x) y sen(x) en el intervalo [0, π8 ] en una aritmética de 32 bits Solución: En primer lugar, la función cos(x) la desarrol- lamos por serie de Taylor como cos(x) u Pn(x) = 1− x2 2! + x4 4! + ....+ (−1)n x 2n (2n)! y el error máximo cometido por el desarrollo de Taylor en un punto x ∈ [0, π8 ] es | Pn(x)− cos(x) |≤ sen( π 8 ) (x) 2n+1 (2n+ 1)! donde ξ ∈ [0, π8 ]. Para que la aproximación sea la mejor dentro de una aritmética de 32 bits tiene que cumplirse | Pn(x)− cos(x) | cos(x) ≤ 2−24 = 5. 96× 10−8 por tanto: | Pn(x)− cos(x) | cos(x) ≤ tan(π 8 ) ¡ π 8 ¢2n+1 (2n+ 1)! para n = 3 se tiene que tan( π 8 ) ¡ π 8 ¢2n+1 (2n+ 1)! = 1. 18 × 10−7 con lo cual ya estamos muy cerca de la precisión óptima. Para n = 4 tan( π 8 ) ¡ π 8 ¢2n+1 (2n+ 1)! = 2. 53× 10−10 por tanto n = 4 determina la mejor aproximación en una aritmética de 32 bits. 8 Análogamente, para la función sen(x) tenemos sen(x) u Pn(x) = x− x3 3! + x5 5! + ....+ (−1)n x 2n+1 (2n+ 1)! y el error máximo cometido por el desarrollo de Taylor en un punto x ∈ [0, π4 ] es | Pn(x)− sen(x) |≤ sen( π 8 ) (x)2n+2 (2n+ 2)! donde ξ ∈ [0, π8 ]. Para que la aproximación sea la mejor dentro de una aritmética de 32 bits tiene que cumplirse | Pn(x)− sen(x) | sen(x) ≤ 2−24 = 5. 96× 10−8 por otro lado, en el intervalo [0, π8 ] se verifica sen ¡ π 8 ¢ π 8 x ≤ sen(x) por tanto: | Pn(x)− sen(x) | sen(x) ≤ ¡ π 8 ¢2n+2 (2n+ 2)! para n = 3 se tiene que¡ π 8 ¢2n+2 (2n+ 2)! = 1. 402 679 863× 10−8 por tanto n = 3 determina la mejor aproximación en una aritmética de 32 bits. Problema 34 (1 puntos) Como se puede obtener la fun- ción yx, donde x, y son números reales, utilizando las fun- ciones ex y ln(x). Solución: Se utiliza la equivalencia yx = ex ln y ANALISIS NUMERICO MATRICIAL I Problema 35 Calcular el número de operaciones básicas (sumas, restas, multiplicaciones y divisiones) en función de la dimensión N necesarias para realizar un remonte para resolver un sistema A0u = b0 donde A0 es una matriz triangular superior. Solución: Escribimos la matriz A0 de la siguiente manera,⎛⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 · · · a1,n−2 a1,n−1 a1n 0 a22 · · · a2,n−2 a2,n−1 a2n ... ... . . . ... ... ... 0 0 0 an−2,n−2 an−2,n−1 an−2,n 0 0 0 0 an−1,n−1 an−1,n 0 0 0 0 0 an,n ⎞⎟⎟⎟⎟⎟⎟⎟⎠ En el remonte se empiezan a calcular los ui de abajo hacia arriba. Las operaciones que se realizan vienen dadas por: un = bn ann un−1 = bn−1−an−1,nun an−1,n−1 un−2 = bn−2−(an−2,nun+an−2,n−1un−1) an−2,n−2 un−3 = bn−3−(an−3,nun+an−3,n−1un−1+an−3,n−2un−2) an−3,n−3 ... En la siguiente tabla se muestra el número de opera- ciones que se realizan en cada iteración: Sumas Multiplic. Divisiones Total n− 1 n− 1 1 2n− 1 ... ... ... ... 3 3 1 7 2 2 1 5 1 1 1 3 0 0 1 1 A partir de esta tabla podemos calcular el total de operaciones sumando por columnas: Sumas = 0 + 1 + 2 + 3 + . . .+ n− 1 = (n−1)n2 Multiplicac. = 0 + 1 + 2 + 3 + . . .+ n− 1 = (n−1)n2 Divisiones = 1 + 1 + 1 + 1 + . . .+ 1 = n Total = 1 + 3 + 5 + 7 + . . .+ 2n− 1 = = Sumas+Multiplicac.+Divisiones = = (n−1)n2 + (n−1)n 2 + n = n 2 El orden del algoritmo es entonces O(n2). Problema 36 Resolver por el método de Gauss el sistemaµ −1 2 2 −1 ¶µ x y ¶ = µ 3 0 ¶ Solución:µ −1 22 −1 ¶µ x y ¶ = µ 3 0 ¶ →µ 2 −1 −1 2 ¶µ x y ¶ = µ 0 3 ¶ →µ 2 −1 0 32 ¶µ x y ¶ = µ 0 3 ¶ → 32y = y → y = 2 → x = 22 = 1 9 Problema 37 Calcular el número de operaciones básicas necesarias para descomponer el sistema Au = b en el sis- tema A0u = b0 utilizando el método de Gauss, y teniendo en cuenta la siguiente relación M−1X k=1 k2 = 1 3 M3 − 1 2 M2 + 1 6 M Solución: A = ⎛⎜⎜⎜⎝ a11 a12 · · · a1n a21 a22 · · · a2n ... ... . . . ... an,1 an,2 · · · an,n ⎞⎟⎟⎟⎠ En cada iteración se realizan las siguientes opera- ciones: Para cada iteración (i): Para cada fila (j) ∗ ³ aii aii ´ ∗aj1 − ai1 ³ aji aii ´ . . . ajn − ain ³ aji aii ´ En la primera iteración, este proceso se repite N − 1 veces (para las N − 1 j-filas inferiores). En la segunda, se repite N −2 veces, y así sucesivamente hasta la penúltima fila, en donde sólo se realiza una vez. Iteración Fila Division. Multiplic. Sumas 1a 2a 3a ... na 1 1 ... 1 n n ... n n n ... n 2a 3a 4a ... na 1 1 ... 1 n− 1 n− 1 ... n− 1 n− 1 n− 1 ... n− 1 ... ... ... ... ... (n− 1)a na 1 2 2 A continuación obtenemos el total de operaciones en cada iteración sumando por columnas: 1a Iteración: Divisiones = 1 + 1 + . . .+ 1 = n− 1 Multiplicac. = n+ n+ . . .+ n = n(n− 1) Sumas = n+ n+ . . .+ n = n(n− 1) 2a Iteración: Divisiones = 1 + 1 + . . .+ 1 = n− 2 Multiplicac. = (n− 1) + (n− 1) + . . .+ (n− 1) = =(n− 1)(n− 2) Sumas = (n−1)+(n−1)+. . .+(n−1) = (n−1)(n−2) ... (n-1)a Iteración: Divisiones = 1 Multiplicac. = 2 Sumas = 2 Total operaciones1: Divisiones = (n−1)+(n−2)+(n−3)+. . .+1 = n(n−1)2 Multiplicac. = n(n− 1) + (n− 1)(n− 2) + . . .+ 2 = =((n−1)+1)(n−1)+((n−2)+1)(n−2)+ . . . = =(n− 1)2 + (n− 1) + (n− 2)2 + (n− 2) + . . . = =2n 3−3n2+n 6 + (n−1)n 2 = n3−n 3 Sumas = n(n− 1) + (n− 1)(n− 2) + . . .+ 2 = n3−n3 Total=Sumas+Multiplicac.+Divisiones = = n 3−n 3 + n(n−1) 2 = 2 3n 3 + 12n 2 − 76n El orden del algoritmo es entonces O(2n33 ). Problema 38 Escribir el pseudocódigo del algoritmo de la funcion IDESCENSO(A, b, u, N) que resuelve un sis- tema donde A es una matriz triangular inferior, b es el vector de términos independientes, u el vector solución, N es la dimensión del sistema La función devuelve 0 si ter- mina correctamente y 1 en caso contrario. Nota Impor- tante: Las líneas de código tienen que ir todas numeradas y no pueden superar las 12 lineas de instrucciones como máximo. Solución: 01 IDESCENSO(matriz real A,vector real b,vector real u,entero N) 02 para variable entera I ← 0 hasta N-1 hacer 03 si(A(I, I) == 0) entonces 04 devolver 1 05 finsi 06 u(I) = b(I) 07 para variable entera J ← 0 hasta I − 1 hacer 08 u(I) = u(I)−A(I, J) ∗ u(J) 11 + 22 + 32 + ...+ (n− 1)2 = 2n 3−3n2+n 6 1 0 09 fin mientras 10 u(I) = u(I)/A(I, I) 11 fin mientras 12 devolver 0 Problema 39 Resolver por el método de Gauss el sigu- iente sistema de ecuaciones⎛⎝ 0 −1 2−1 2 −1 2 −1 0 ⎞⎠⎛⎝ u1u2 u3 ⎞⎠ = ⎛⎝ 10 1 ⎞⎠ Solución: Pasos en la descomposición por Gauss: 1. Intercambiamos la tercera fila con la primera:⎛⎝ 0 −1 2 1−1 2 −1 0 2 −1 0 1 ⎞⎠−−−−→pivoteo ⎛⎝ 2 −1 0 1−1 2 −1 0 0 −1 2 1 ⎞⎠ 2. Hacemos ceros en la primera columna³ filaj − fila1 · aj1a11 ; j > 1 ´ :⎛⎝ 2 −1 0 1−1 2 −1 0 0 −1 2 1 ⎞⎠ −−−→ceros ⎛⎝ 2 −1 0 10 3 −2 1 0 −1 2 1 ⎞⎠ 3. Hacemos ceros en la segunda columna³ filaj − fila2 · aj1a11 ; j > 2 ´ :⎛⎝ 2 −1 0 10 3 −2 1 0 −1 2 1 ⎞⎠ −−−→ceros ⎛⎝ 2 −1 0 10 3 −2 1 0 0 4 4 ⎞⎠ 4. Realizamos el remonte, y obtenemos como solución: u3 = 4 4 = 1 u2 = 1−2u3 −1 = −1 −1 = 1 u1 = 1+u2 2 = 2 2 = 1 Problema 40 Demostrar que si A = B ·Bt (B triangular inferior) y |B| 6= 0, entonces A es simétrica y definida positiva Solución:Tenemos que demostrar, por una parte, que At = A (A simétrica) y, por otra, que x̄tAx̄ > 0 (A definida positiva2). 1. Simétrica: At = (B ·Bt)t = (B ·Bt)t = (Bt)tBt = B ·Bt = A 2Matriz definida positiva: ∀x̄ 6= 0 =⇒ x̄tAx̄ > 0. Esta es la definición formal. De forma práctica, se comprueba que los menores principales de la matriz sean positivos. También se cumple si todos sus autovalores son positivos: x̄tAx̄ = x̄tλx̄ = λx̄tx̄ > 0. 2. Definida positiva: Como |B| 6= 0, si Bx̄ = 0 =⇒ x̄ = 0 Una matriz se dice definida positiva si se cumple que ∀x̄ 6= 0, x̄tAx̄ > 0 =⇒ =⇒ x̄tAx̄ = x̄tBBtx̄ = (Btx̄)Btx̄ = = ȳtȳ = P y2i > 0 Problema 41 Descomponer la siguiente matriz A por el método de Cholesky A = ⎛⎝ 1 1 41 5 6 4 6 26 ⎞⎠ Solución: La descomposición por el método de Cholesky tiene la forma siguiente: A = B ·Bt, donde la matriz B es triangular inferior. B = ⎛⎝ b11 0 0b21 b22 0 b31 b32 b33 ⎞⎠ Bt = ⎛⎝ b11 b21 b310 b22 b32 0 0 b33 ⎞⎠ Cálculo de los elementos de la matriz B : A = B ·Bt = = ⎛⎝ b11 0 0b21 b22 0 b31 b32 b33 ⎞⎠⎛⎝ b11 b21 b310 b22 b32 0 0 b33 ⎞⎠ = = ⎛⎝ b211 b11b21 b11b31b11b21 b221 + b222 b21b31 + b22b32 b11b31 b21b31 + b22b32 b 2 31 + b 2 32 + b 2 33 ⎞⎠ Igualamos los elementos de la matriz anterior con los elementos de la matriz A y se obtienen los siguientes re- sultados: b211 = 1 b11 = 1 b11b21 = 1 b21 = 1 b11 = 1 b11b31 = 4 b31 = 4 b11 = 4 1 1 b221 + b 2 22 = 5 b22 = ± p (5− b221) = p (4) = 2 b21b31 + b22b32=6 b32 = 6−b21b31 b22 = 6−42 = 1 b231 + b 2 32 + b 2 33 = 26 b33 = p (26− b231 − b232) = p (26− 16− 12) = 3 La descomposición queda de la siguiente manera: A = B ·Bt = = ⎛⎝ 1 0 01 2 0 4 1 3 ⎞⎠⎛⎝ 1 1 40 2 1 0 0 3 ⎞⎠ Problema 42 Calcular el número de operaciones nece- sarias para resolver un sistema por el método de Cholesky. Solución: Las operaciones que se realizan en cada it- eración vienen dadas por: Iteración Operaciones i = 1 j = 1 : b11 = √ a11 j = 2 : b21 = a21 b11 ... j = n : bn1 = an1 b11 i = 2 j = 2 : b22 = p a22 − b221 j = 3 : b32 = a32−b21b31 b22 ... j = n : bn2 = an2−b21bn1 b22 ... ... i = n j = n : bnn = q ann − ¡ b2n1 + . . .+ b 2 n,n−1 ¢ En la siguiente tabla se muestra de forma esquemati- zada, el número de operaciones en cada iteración: Iteración Sumas Multiplic. Divisiones i = 1 0 0 ... 0 0 0 ... 0 0 1 ... 1 n−1 i = 2 1 1 ... 1 n−1 1 1 ... 1 n−1 0 1 ... 1 n−2 i = 3 2 2 ... 2 2(n−2) 2 2 ... 2 2(n−2) 0 1 ... 1 n−3 ... ... ... ... i = n n− 1 n− 1 0 El total de operaciones se obtiene sumando los totales parciales de la tabla anterior: Sumas =Multiplicac. = = (n− 1)+ 2 (n− 2)+ 3(n− 3)+ . . .+(n− 1) = = Pn−1 i=1 i(n− i) = n 3−n 6 Divisiones = (n− 1) + (n− 2) + (n− 3) + . . .+ 1 = = Pn−1 i=1 i = n(n−1) 2 El resultado final es: Total=Sumas+Multiplicac.+Divisiones = = 2n 3−n 6 + n(n−1) 2 = 1 3n 3 − 56n+ 1 2n 2 El orden del algoritmo es O ³ n3 3 ´ Problema 43 Demostrar que a partir de un método para resolver sistemas de ecuaciones se puede construir de forma inmediata un método para calcular la inversa A−1 de una matriz A. Solución: AA−1 = Id = ⎛⎜⎜⎜⎝ 1 0 · · · 0 0 1 · · · 0 ... ... . . . ... 0 0 · · · 1 ⎞⎟⎟⎟⎠ Si expresamos la matriz inversa de la siguiente man- era: 12 A ⎛⎜⎜⎜⎝ c11 c12 · · · c1n c21 c22 · · · c2n ... ... . . . ... cn1 cn2 · · · cnn ⎞⎟⎟⎟⎠ = ⎛⎜⎜⎜⎝ 1 0 · · · 0 0 1 · · · 0 ... ... . . . ... 0 0 · · · 1 ⎞⎟⎟⎟⎠ , se pueden calcular las columnas de esa matriz a partir de N sistemas de ecuaciones de la siguiente forma: A ⎛⎜⎜⎜⎝ c11 c21 ... cn1 ⎞⎟⎟⎟⎠ = ⎛⎜⎜⎜⎝ 1 0 ... 0 ⎞⎟⎟⎟⎠ A ⎛⎜⎜⎜⎝ c12 c22 ... cn2 ⎞⎟⎟⎟⎠ = ⎛⎜⎜⎜⎝ 0 1 ... 0 ⎞⎟⎟⎟⎠ ... A ⎛⎜⎜⎜⎝ c1n c2n ... cnn ⎞⎟⎟⎟⎠ = ⎛⎜⎜⎜⎝ 0 0 ... 1 ⎞⎟⎟⎟⎠ , c.q.d. Problema 44 Demostrar el algoritmo de Crout para de- scomponer matrices tridiagonales. Solución: Consideremos la matriz tridiagonal siguiente: A = ⎛⎜⎜⎜⎜⎜⎝ a1 b1 0 · · · 0 c1 a2 b2 · · · 0 0 c2 a3 b3 0 ... ... ... . . . bn−1 0 0 0 cn−1 an ⎞⎟⎟⎟⎟⎟⎠ La descomposición por el método de Crout genera dos matrices de la forma: A = ⎛⎜⎜⎝ l1 0 . 0 m1 l2 . 0 0 . . 0 0 . mn−1 ln ⎞⎟⎟⎠ ⎛⎜⎜⎝ 1 u1 . 0 0 1 . 0 0 . . un−1 0 . 0 1 ⎞⎟⎟⎠ = = ⎛⎜⎜⎝ l1 l1u1 0 . 0 m1 m1u1 + l2 l2u2 . 0 0 . . . ln−1un−1 0 0 . mn−1 mn−1un−1 + ln ⎞⎟⎟⎠ Igualando ambas matrices y despejando los elementos li, ui y mi, l1u1 = b1 l1 = a1,u1 = b1 l1 ,m1= c1 m1u1 + l2 = a2 l2u2 = b2 l2 = a2 −m1u1,u2 = b2 l2 ,m2 = c2 ... mn−2un−2 + ln−1 = an−1 ln−1un−1 = bn−1 ln−1 = an−1 −mn−2un−2,un−1 = bn−1 ln−1 ,mn−1 = cn−1 mn−1un−1 + ln = an ln = an −mn−1un−1 El algoritmo queda de la siguiente manera: l1 = a1 u1 = b1 l1 Para i = 2, . . . , n− 1 mi−1 = ci−1 li = ai −mi−1ui−1 ui = bi li Fin Para mn−1 = cn−1 ln = an −mn−1un−1 Problema 45 Resolver utilizando el método de Crout el siguiente sistema de ecuaciones⎛⎝ 2 4 0−1 0 4 0 −1 0 ⎞⎠⎛⎝ xy z ⎞⎠ = ⎛⎝ 63 −1 ⎞⎠ Solución: Aplicando el algoritmo del problema anterior, obtenemos los siguientes resultados: i=1 l1 = 2 u1 = 4 2 = 2 i=2 m1 = −1 l2 = 0− 2 (−1) = 2 u2 = 4 2 = 2 i=3 m2 = −1 l3 = 0− 2 (−1) = 2 1 3 Sustituyendo estos valores en las matrices de Crout, la descomposición queda: A = L · U = ⎛⎝ 2 0 0−1 2 0 0 −1 2 ⎞⎠ · ⎛⎝ 1 2 00 1 2 0 0 1 ⎞⎠ Para resolver el sistema, se tiene en cuenta lo sigu- iente: Ax = b LUx = b (Ux = y) y nos queda un sistema de la forma: Ly = b Calculamos el valor de y a partir del sistema anterior:⎛⎝ 2 0 0−1 2 0 0 −1 2 ⎞⎠⎛⎝ y1y2 y3 ⎞⎠ = ⎛⎝ 63 −1 ⎞⎠ , aplicando un algoritmo de descenso,⎛⎝ y1y2 y3 ⎞⎠ = ⎛⎝ 623+y1 2−1+y2 2 ⎞⎠ = ⎛⎝ 33 1 ⎞⎠ Calculamos el vector x por remonte: Ux = y⎛⎝ 1 2 00 1 2 0 0 1 ⎞⎠⎛⎝ x1x2 x3 ⎞⎠ = ⎛⎝ 33 1 ⎞⎠ ⎛⎝ x1x2 x3 ⎞⎠ = ⎛⎝ 3− 2x23− 2x3 1 ⎞⎠ = ⎛⎝ 11 1 ⎞⎠ quedándonos la solución final x = ¡ 1 1 1 ¢ Problema 46 Calcular el número de operaciones nece- sarias para resolver un sistema tridiagonal por el método de Crout. Solución: Las operaciones que se realizan en cada it- eración vienen dadas por: Iteración Operaciones i = 1 l1 = a1;u1 = b1 l1 i = 2 m1 = c1; l2 = a2 −m1u1;u2 = b2l2 ... ... i = n ln = an −mn−1un−1 En la siguiente tabla se muestra el número de op- eraciones en cada iteración: Iteración Sumas Multiplic. Divisiones i = 1 0 0 1 i = 2 1 1 1 i = 3 1 1 1 ... ... ... ... i = n 1 1 0 El total de operaciones se obtiene de la tabla anterior como: Sumas =Multiplicac. = Divisiones = = 1 + 1 + . . .+ 1 = (n− 1) Total=Sumas+Multiplicac.+Divisiones = = 3 (n− 1) El orden del algoritmo es O (3n) DIFERENCIACION E INTEGRACION NUMERICA Problema 47 Calcular analíticamente y numéricamente la matriz gradiente en el punto (1, 1) (utilizar h = 0.1) de la función: f(x, y) = ½ x2 + y2 − 1 x− y Solución: Analíticamente ∇f(x, y) = µ 2x 2y 1 −1 ¶ →∇f(1, 1) = µ 2 2 1 −1 ¶ Numéricamente, si llamamos f1(x, y) = x2 + y2 − 1 y f2(x, y) = x−y, entonces aplicando las máscaras vistas en teoría tenemos f1(x, y) = x− y ∂f1(1,1) ∂x = 1 4(0.1) ¡ 2− √ 2 ¢ (f1(1 + 0.1, 1− 0.1)− f1(1− 0.1, 1− 0.1))+ 1 4(0.1) ¡ 2− √ 2 ¢ (f1(1 + 0.1, 1 + 0.1)− f1(1− 0.1, 1 + 0.1))+ 1 4(0.1)2 ¡√ 2− 1 ¢ (f1(1 + 0.1, 1)− f1(1− 0.1, 1)) = 0.585 79 + 0.585 79 + 0.828 43 = 2.0 ∂f1(1,1) ∂y = 1 4(0.1) ¡ 2− √ 2 ¢ (f1(1− 0.1, 1 + 0.1)− f1(1− 0.1, 1− 0.1))+ 1 4(0.1) ¡ 2− √ 2 ¢ (f1(1 + 0.1, 1 + 0.1)− f1(1 + 0.1, 1− 0.1))+ 1 4(0.1)2 ¡√ 2− 1 ¢ (f1(1, 1 + 0.1)− f1(1, 1− 0.1)) = 1 4 0.585 79 + 0.585 79 + 0.828 43 = 2.0 De la misma forma, para f2(x, y) tenemos ∂f2(1,1) ∂x = 1 4(0.1) ¡ 2− √ 2 ¢ (f2(1 + 0.1, 1− 0.1)− f2(1− 0.1, 1− 0.1))+ 1 4(0.1) ¡ 2− √ 2 ¢ (f2(1 + 0.1, 1 + 0.1)− f2(1− 0.1, 1 + 0.1))+ 1 4(0.1)2 ¡√ 2− 1 ¢ (f2(1 + 0.1, 1)− f2(1− 0.1, 1)) = 0.292 89 + 0.292 89 + 0.414 21 = 1 ∂f2(1,1) ∂y = 1 4(0.1) ¡ 2− √ 2 ¢ (f2(1− 0.1, 1 + 0.1)− f2(1− 0.1, 1− 0.1))+ 1 4(0.1) ¡ 2− √ 2 ¢ (f2(1 + 0.1, 1 + 0.1)− f2(1 + 0.1, 1− 0.1))+ 1 4(0.1)2 ¡√ 2− 1 ¢ (f2(1, 1 + 0.1)− f2(1, 1− 0.1)) = 0.292 89 + 0.292 89 + 0.414 21 = 1 Con lo cual, en este caso la matriz gradiente calculada numéricamente coincide con la calculada analíticamente. → ∇f(1, 1) = µ 2 2 1 −1 ¶ Problema 48 Dados 3 puntos distintos xl, xi, xr demostrar que la fórmula: f 0(xi) ≈ (xi − xl) f(xr)−f(xi)xr−xi + (xr − xi) f(xi)−f(xl) xi−xl xr − xl aproxima la derivada de f 0(xi) con un orden de aproxi- mación de 2. Solución: Evaluamos el desarrollo de Taylor de la función en los puntos xr, xl : f (xl) = f(xi) + f 0(xi)(xl − xi)+ +f 00(xi) 2! (xl − xi)2 + f 000(xi) 3! (xl − xi)3 f (xr) = f(xi) + f 0(xi)(xr − xi)+ +f 00(xi) 2! (xr − xi)2 + f 000(xi) 3! (xr − xi)3 (xr − xi) f(xl)−f(xi)(xl−xi) = (xr − xi)[f 0(xi)+ +f 00(xi) 2! (xl − xi) + f 000(xi) 3! (xl − xi)2] (xi − xl)f(xi)−f(xr)(xi−xr) = (xi − xl)[f 0(xi)+ +f 00(xi) 2! (xr − xi) + f 000(xi) 3! (xr − xi)2] Sumamos las expresiones anteriores y nos queda: (xr − xi) f(xl)−f(xi)(xl−xi) + (xi − xl) f(xr)−f(xi) (xr−xi) = = (xr − xi)f 0(xi) + (xi − xl)f 0(xi)+ +f 00(xi) 2! (xi − xl)(xr − xi)+ +f 00(xi) 2! (xr − xi)(xl − xi)+ +f 000(xi) 3! (xr − xi)(xl − xi)2+ +f 000(xi) 3! (xi − xl)(xr − xi)2 Agrupamos por las derivadas de la función: (xr − xi)f(xl)−f(xi)(xl−xi) + (xi − xl) f(xr)−f(xi) (xr−xi) = = (xr − xl)f 0(xi) + f 00(xi) 2! · 0+ +f 000(xi) 3! (xr − xi)(xl − xi) ((xl − xi)− (xr − xi)) El término de la tercera derivada nos da el orden de la fórmula: (xr − xl)f 0(xi) = (xr − xi)f(xl)−f(xi)(xl−xi) + +(xi − xl) f(xr)−f(xi)(xr−xi) +O(h 3) f 0(xi) ≈ (xi−xl) f(xr)−f(xi)xr−xi +(xr−xi) f(xi)−f(xl) xi−xl xr−xl +O(h 2) Problema 49 Dados 3 puntos distintos xl, xi, xr calcular el polinomio de Lagrange que interpola a f(x) en esos 3 puntos, calcular la derivada de ese polinomio en xi y com- probar que da la misma fórmula que la presentada en el problema anterior. Solución: El polinomio de Lagrange es: f (x) = (x−xr)(x−xl)(xi−xr)(xi−xl)f (xi) + (x−xi)(x−xl) (xr−xi)(xr−xl)f (xr)+ + (x−xi)(x−xr)(xl−xi)(xl−xr)f (xl) Derivamos la expresión anterior y obtenemos: f 0 (x) = (x−xl)+(x−xr)(xi−xr)(xi−xl)f (xi) + (x−xl)+(x−xi) (xr−xi)(xr−xl)f (xr)+ + (x−xr)+(x−xi)(xl−xi)(xl−xr)f (xl) Evaluamos la derivada en el punto xi y desarrollamos hasta obtener el resultado: f 0 (xi) = (xi−xl)+(xi−xr) (xi−xr)(xi−xl) f (xi)+ + (xi−xl)+(xi−xi)(xr−xi)(xr−xl) f (xr) + (xi−xr)+(xi−xi) (xl−xi)(xl−xr) f (xl) f 0 (xi) = f(xi) (xi−xr) + f(xi) (xi−xl) + (xi−xl)f(xr) (xr−xi)(xr−xl)+ + (xi−xr)f(xl)(xl−xi)(xl−xr) extraemos el factor (xr − xl), (xr − xl) f 0 (xi) = (xr−xl)f(xi)(xi−xr) + (xr−xl)f(xi) (xi−xl) + 1 5 + (xi−xl)f(xr)(xr−xi) − (xi−xr)f(xl) (xl−xi) (xr − xl) f 0 (xi) = −xrf(xi)+xlf(xi)(xr−xi) + xrf(xi)−xlf(xi) (xi−xl) + + (xi−xl)f(xr)(xr−xi) − (xr−xi)f(xl) (xi−xl) agrupamos términos, (xr − xl) f 0 (xi) = ³ (xi−xl)f(xr) (xr−xi) − (xi−xl)f(xi) (xr−xi) ´ + + ³ (xr−xi)f(xi) (xi−xl) − (xr−xi)f(xl) (xi−xl) ´ + xif(xi)(xr−xi) − xrf(xi)(xr−xi) + xif(xi) (xi−xl) − xlf(xi) (xi−xl) (xr − xl) f 0 (xi) = (xi−xl)(f(xr)−f(xi))(xr−xi) + + (xr−xi)(f(xi)−f(xl))(xi−xl) + xif(xi) (xr−xi) − xrf(xi) (xr−xi)+ + xif(xi)(xi−xl) − xlf(xi) (xi−xl) (xr − xl) f 0 (xi) = (xi−xl)(f(xr)−f(xi))(xr−xi) + + (xr−xi)(f(xi)−f(xl))(xi−xl) + +xif(xi)(xi−xl)−xrf(xi)(xi−xl)(xr−xi)(xi−xl) + +xif(xi)(xr−xi)−xlf(xi)(xr−xi)(xr−xi)(xi−xl) simplificando, f 0 (xi) = (xi−xl)(f(xr)−f(xi)) (xr−xi) + (xr−xi)(f(xi)−f(xl)) (xi−xl) (xr−xl) Problema 50 Calcular una aproximación de la derivada tercera f 000(xi) de una función f(x) en un punto xi, uti- lizando f(xi), f(xi + h), f(xi − h), f(xi − 2h) Solución: a→ f(xi+ h) = f(xi)+ hf 0(xi)+ h 2 2 f 00(xi)+ h3 6 f 000(xi)+ O(h4) 1. b → f(xi − h) = f(xi) − hf 0(xi) + h 2 2 f 00(xi) − h3 6 f 000(xi) +O(h 4) c → f(xi − 2h) = f(xi) − 2hf 0(xi) + 2h2f 00(xi) − 4h3 3 f 000(xi) +O(h 4) Sistema: ⎛⎝ a− b− 2c = 0a 2 + b 2 + 2c = 0 a 6 − b 6 − 4c 3 = 1 ⎞⎠ Solución: a = 1, b = 3, c = −1. f 000(xi) = af(xi+h)+bf(xi−h)+cf(xi−2h)−(a+b+c)f(xi) h3 = f(xi+h)+3f(xi−h)−f(xi−2h)−3f(xi) h3 +O(h) Problema 51 Dados 3 puntos. Demostrar que la fórmula f 00(xi) ≈ 2 f(xr)−f(xi) xr−xi − f(xi)−f(xl) xi−xl xr − xl aproxima la derivada segunda de f(x) en xi con un orden de aproximación de 1. Solución: Desarrollo de Taylor de la función en el punto xi y evaluación en xr y xl: f (xr) ≈ f (xi) + f 0 (xi) (xr − xi)+ +f 00(xi) 2! (xr − xi) 2 + f 000(xi) 3! (xr − xi) 3 f (xl) ≈ f (xi) + f 0 (xi) (xl − xi)+ +f 00(xi) 2! (xl − xi) 2 + f 000(xi) 3! (xl − xi) 3 Extraemos en ambas ecuaciones:f(xr)−f(xi) (xr−xi) ≈ f 0 (xi) + f 00(xi) 2! (xr − xi)+ +f 000(xi) 3! (xr − xi) 2 f(xl)−f(xi) (xl−xi) ≈ f 0 (xi) + f 00(xi) 2! (xl − xi)+ +f 000(xi) 3! (xl − xi) 2 Restamos las expresiones anteriores: f(xr)−f(xi) (xr−xi) − f(xi)−f(xl) (xi−xl) ≈ f 00(xi) 2! (xr − xl)+ +f 000(xi) 3! ³ (xr − xi)2 − (xl − xi)2 ´ Despejamos la segunda derivada y obtenemos: f 00 (xi) ≈ 2 f(xr)−f(xi) (xr−xi) − f(xi)−f(xl) (xi−xl) xr−xl − −2 f000(xr) 3! ((xr−xi) 2−(xl−xi)2) xr−xl f 00 (xi) ≈ 2 f(xr)−f(xi) (xr−xi) − f(xi)−f(xl) (xi−xl) xr−xl +O(h) Problema 52 Considerar en el problema anterior que xl = xi−h, y xr = xi+h. Deducir como queda la fórmula anterior para aproximar la derivada segunda, y demostrar que en este caso el orden de aproximación es 2. Solución: Sustituyendo xl = xi − h, y xr = xi + h, tenemos: f 00 (xi) ≈ 2 f(xr)−f(xi) (xi+h−xi) − f(xi)−f(xl) (xi−xi+h) xi+h−xi+h = = f(xr)−f(xi) h − f(xi)−f(xl) h h = 1 6 = f(xr)−f(xi)−f(xi)−f(xl)h2 − f 000(xi) 3! (h− h) + O(h2) = = f(xr)−2f(xi)−f(xl)h2 +O(h 2) La aproximación de la segunda derivada queda de la forma, f 00 (xi) ≈ f (xr)− 2f (xi)− f (xl) h2 Problema 53 Dados 3 puntos xl < xi < xr calcular el polinomio de Lagrange que interpola a f(x) en esos 3 pun- tos, calcular la derivada segunda de ese polinomio en xi y comprobar que da la misma fórmula que utilizando los desarrollos de Taylor. Solución: Por las diferencias divididas de Newton obten- emos lo siguiente: xl → f (xl) xi → f (xi) À f(xi)−f(xl) xi−xl xr → f (xr) À f(xr)−f(xi) xr−xi + f(xr)−f(xi) xr−xi − f(xi)−f(xl) xi−xl xr − xl Polinomio de Lagrange: P (x) ' f(xl) + f(xi)−f(xl)xi−xl (x− xl)+ + f(xr)−f(xi) xr−xi − f(xi)−f(xl)xi−xl xr−xl (x− xl) (x− xi) Derivamos el polinomio: P 0 (x) ' f(xi)−f(xl)xi−xl + f(xr)−f(xi) xr−xi − f(xi)−f(xl)xi−xl xr−xl (x− xl)+ + f(xr)−f(xi) xr−xi − f(xi)−f(xl)xi−xl xr−xl (x− xi) Calculamos la segunda derivada, obteniendo: P 00 (x) ' 2 f(xr)−f(xi) xr−xi − f(xi)−f(xl)xi−xl xr−xl , c.q.d. Problema 54 Calcular una aproximación de la derivada primera y segunda de una función f(x) en x = 0, teniendo en cuenta que f(0) = 1, f(1) = 0, f(4) = 9 Solución: f 0(xi) ≈ (xi−xl) f(xr)−f(xi)xr−xi +(xr−xi) f(xi)−f(xl) xi−xl xr−xl = = (0−1) f(4)−f(0)4−0 +(4−0) f(0)−f(1) 0−1 4−1 = = − 9−14 +4 1−0 −1 3 = −2−4 3 = −63 = −2 f 00 (xi) ≈ 2 f(xr)−f(xi) (xr−xi) − f(xi)−f(xl) (xi−xl) xr−xl = = 2 9−1 (4−0)− 1−0 (0−1) 4−1 = 2 2+1 3 = 2 Problema 55 Demostrar, utilizando el desarrollo de Taylor, que las siguientes expresiones son discretizaciones del laplaciano: ∆F = Fi+1,j+1 + Fi−1,j+1 + Fi−1,j−1 + Fi+1,j−1 − 4Fi,j 2h2 ∆F = Fi+1,j + Fi−1,j + Fi,j+1 + Fi,j−1 − 4Fi,j h2 Solución: A partir del desarrollo de Taylor de la función F , se obtiene lo siguiente: F (x+ h, y + h) = F (x, y) + hFx + hFy+ +12h 2 (Fxx + 2Fxy + Fyy) F (x− h, y − h) = F (x, y)− hFx − hFy+ +12h 2 (Fxx + 2Fxy + Fyy) F (x− h, y + h) = F (x, y)− hFx + hFy+ +12h 2 (Fxx − 2Fxy + Fyy) F (x+ h, y − h) = F (x, y) + hFx − hFy+ +12h 2 (Fxx − 2Fxy + Fyy) Sumamos estas cuatro ecuaciones, F (x+ h, y + h)+F (x− h, y − h)+F (x− h, y + h)+ +F (x+ h, y − h) = 4F (x, y) + 2h2 (Fxx + Fyy) Fxx + Fyy = = F (x+h,y+h)+F (x−h,y−h)+F (x−h,y+h)+F (x+h,y−h)−4F (x,y)2h2 , discretizando se obtiene el resultado esperado, ∆F = Fi+1,j+1 + Fi−1,j+1 + Fi−1,j−1 + Fi+1,j−1 − 4Fi,j 2h2 Para demostrar la segunda igualdad, tomamos las siguientes ecuaciones: F (x+ h, y) = F (x, y) + hFx + h2 2 Fxx F (x− h, y) = F (x, y)− hFx + h 2 2 Fxx F (x, y + h) = F (x, y) + hFy + h2 2 Fyy F (x, y − h) = F (x, y)− hFy + h 2 2 Fyy Sumamos estas expresiones y obtenemos: F (x+ h, y) + F (x− h, y) + F (x, y + h)+ +F (x, y − h) = 4F (x, y) + h2Fxx + h2Fyy Fxx + Fyy = F (x+h,y)+F (x−h,y)+F (x,y+h)+F (x,y−h)−4F (x,y) h2 , discretizando ∆F = Fi+1,j + Fi−1,j + Fi,j+1 + Fi,j−1 − 4Fi,j h2 1 7 Problema 56 Calcular una aproximación del lapla- ciano de una función F (x, y) en el punto (x, y) = (0, 0) conociendo los siguientes valores: F (0, 0) = 0, F (12 , 0) = 1 4 , F (− 1 2 , 0) = 1 4 , F (0, 1 2) = 1 4 , F (0,− 1 2 ) = 1 4 , F (12 , 1 2) = 1 2 , F (− 1 2 ,− 1 2) = 1 2 , F (− 1 2 , 1 2 ) = 1 2 , F ( 1 2 ,− 1 2) = 1 2 Solución: Si representamos estos valores en una tabla, obtenemos lo siguiente: 1 2 1 4 1 2 1 4 0 1 4 1 2 1 4 1 2 El valor de h es 12 . Aproximamos el laplaciano promediando las dos ex- presiones del ejercicio anterior. Si no realizáramos este promediado, no se tendrían en cuenta todos los valores de la función. ∆F = γ Fi+1,j+1+Fi−1,j+1+Fi−1,j−1+Fi+1,j−1−4Fi,j 2h2 + +(1− γ) Fi+1,j+Fi−1,j+Fi,j+1+Fi,j−1−4Fi,jh2 , γ = 23 ∆F (0, 0) = 23 1 2+ 1 2+ 1 2+ 1 2 2 14 + ¡ 1− 23 ¢ 1 4+ 1 4+ 1 4+ 1 4 1 4 = = 83 + 4 3 = 4 Problema 57 Demostrar que las máscaras Fx = 1 4h − ¡ 2− √ 2 ¢ 0 ¡ 2− √ 2 ¢ −2 ¡√ 2− 1 ¢ 0 2 ¡√ 2− 1 ¢ − ¡ 2− √ 2 ¢ 0 ¡ 2− √ 2 ¢ Fy = 1 4h − ¡ 2− √ 2 ¢ −2 ¡√ 2− 1 ¢ − ¡ 2− √ 2 ¢ 0 0 0¡ 2− √ 2 ¢ 2 ¡√ 2− 1 ¢ ¡ 2− √ 2 ¢ dan lugar a una discretización del gradiente tal que su norma euclídea es invariante por rotaciones de 45 grados. Solución: Procedemos de la misma forma que al calcular el valor de γ en el caso del laplaciano. Consideramos una función que tiene los siguientes val- ores en un entorno de un punto (hi0, hj0) : 1 1 1 0 0 0 0 0 0 Calculamos el valor del gradiente en el punto central de la siguiente manera: Fx = γ 0 2h + (1− γ) 0 4h = 0 Fy = (1− γ) −12h + γ −2 4h = − 1 2 γ h − 1 2 1−γ h = − 1 2h ∇1F (hi0, hj0) = (Fx, Fy) = ¡ − 12h , 0 ¢ Rotamos la función anterior 45o : 1 1 0 1 0 0 0 0 0 y calculamos su gradiente: Fx = (1− γ) −12h + γ −1 4h = − 1 2 1−γ h − 1 4 γ h = 1 4 γ−2 h Fy = (1− γ) −12h + γ −1 4h = − 1 2 1−γ h − 1 4 γ h = 1 4 γ−2 h ∇2F (hi0, hj0) = (Fx, Fy) = 14 γ−2 h (1, 1) Calculamos las normas de los gradientes e igualamos: k∇1F (hi0, hj0)k = k∇2F (hi0, hj0)k 1 2h = q 2 ¡ 1 4 γ−2 h ¢2 1 2h = 1 4h √ 2 |γ − 2|( 2√ 2 = − (γ − 2)→ γ = 2− √ 2 2√ 2 = (γ − 2)→ γ = 2 + √ 2 La solución válida es γ = 2− √ 2, ya que el gradiente ∇2F debe ser negativo en sus dos derivadas. Sustituyendo este valor en las expresiones de Fx, Fy tenemos: Fx = (1− γ) Fi+1,j−Fi−1,j2h + +γ Fi+1,j+1−Fi−1,j+1+Fi+1,j−1−Fi−1,j−1 4h = = 2 ¡√ 2− 1 ¢ Fi+1,j−Fi−1,j 4h + + ¡ 2− √ 2 ¢ Fi+1,j+1−Fi−1,j+1+Fi+1,j−1−Fi−1,j−1 4h Fy = (1− γ) Fi,j+1−Fi,j−12h + +γ Fi+1,j+1−Fi+1,j−1+Fi−1,j+1−Fi−1,j−1 4h = = 2 ¡√ 2− 1 ¢ Fi,j+1−Fi,j−1 4h + + ¡ 2− √ 2 ¢ Fi+1,j+1−Fi+1,j−1+Fi−1,j+1−Fi−1,j−1 4h , cuyas máscaras son las que se muestran en el enunciado del problema. Problema 58 Calcular una aproximación del gradi- ente de una función F (x, y) en el punto (x, y) = (0, 0) conociendo los siguientes valores: F (0, 0) = 0, F (12 , 0) = 1 2 , F (− 1 2 , 0) = − 1 2 , F (0, 1 2) = − 1 2 , F (0,− 1 2) = 1 8 1 2 , F ( 1 2 , 1 2) = 0, F (− 1 2 ,− 1 2) = 0, F (− 1 2 , 1 2) = −1, F (12 ,− 1 2) = 1 Solución: Los valores de la función en una tabla quedan de la siguiente manera: 0 12 1−1 2 0 1 2 −1 −12 0 Sustituimos estos valores en las derivadas de la fun- ción: Fx = 2 ¡√ 2− 1 ¢ Fi+1,j−Fi−1,j 4h + + ¡ 2− √ 2 ¢ Fi+1,j+1−Fi−1,j+1+Fi+1,j−1−Fi−1,j−1 4h = = 2 ¡√ 2− 1 ¢ 1 2+ 1 2 4h + ¡ 2− √ 2 ¢ 1+1 4h = = 12 √ 2−1 h + 1 2 2− √ 2 h = 1 2h Fy = 2 ¡√ 2− 1 ¢ Fi,j+1−Fi,j−1 4h + + ¡ 2− √ 2 ¢ Fi+1,j+1−Fi+1,j−1+Fi−1,j+1−Fi−1,j−1 4h = = 2 ¡√ 2− 1 ¢ −1 2 − 1 2 4h + ¡ 2− √ 2 ¢ −1−1 4h = = −12 √ 2−1 h − 1 2 2− √ 2 h = − 1 2h y obtenemos el valor del gradiente: ∇F = µ Fx Fy ¶ = µ 1 2h − 12h ¶ = µ 1 −1 ¶ Este vector nos da la dirección de máximo ascenso, que en este caso será en diagonal hacia arriba a la derecha. Problema 59 Aproximar el valor de la siguiente integral, utilizando las fórmulas de Legendre para n = 2 y n = 3Z 1 −1 ¡ x3 − x4 ¢ dx Cual es el valor exacto de la integral? Solución:R 1 −1 ¡ x3 − x4 ¢ dx ' NP k=0 wkP (xk) P (x) = ¡ x3 − x4 ¢ 1. n = 2 2P k=1 wkP (xk) = = 1 · P (0.5773502692 + 1 · P (−0.5773502692) = = −. 222 22 2. n = 3 3P k=1 wkP (xk)= = 0.55555555555 · P (0.7745966692)+ +0.88888888 · P (0)+ +0.55555555555 · P (−0.7745966692) = = −. 4 El valor exacto de la integral es R 1 −1 ¡ x3 − x4 ¢ dx = −25 = −. 4, que coincide con el valor del segundo caso. La fórmula de integración numérica es exacta hasta el orden 2n − 1, que en el segundo caso es equivalente a 5, con lo que ya se sabía que el valor obtenido sería exacto. Problema 60 Se considera para el intervalo [−1, 1], los puntos x0 = −0.5, x1 = 0 y x2 = 0.5 y los pesos w0 = w1 = w2 = 2/3. Estos puntos y estos pesos se utilizan para aproximar la integral de una función en [−1, 1]. Usar esta fórmula de integración para calcular númericamente la siguiente integral y compararla con el resultado análitico (exacto). Z π 2 −π2 cos(x)dx Solución: 1. R π 2 −π2 cos(x)dx = sin(x)] π 2 −π2 = 1− (−1) = 2 R π 2 −π2 cos(x)dx = R 1 −1 cos( π 2 t) π 2 dt = 2 3 π 2 cos ¡ −π4 ¢ + 2 3 π 2 cos (0) + 2 3 π 2 cos ¡ π 4 ¢ = 13 √ 2π + 13π = 2. 528 2 Problema 61 Encontrar cual sería la fórmula de inte- gración numérica en el intervalo [−1, 1] utilizando un sólo punto de interpolación, y de tal manera que sea exacta para polinomios de grado 1 Solución: La fórmula que usa un único punto se puede expresar como Z 1 −1 f (x) dx ' w0 · f (x0) Vamos a imponer que se exacta para los polinomios f(x) = 1 y f(x) = xZ 1 −1 1dx = 2 = w0 · f (x0) = w0 → w0 = 2 1 9 Z 1 −1 xdx = 0 = w0 · f (x0) = w0 · x0 → x0 = 0 Por lo tanto, la fórmula de integración numérica de Legen- dre es: Z 1 −1 f (x) dx ' 2 · f (0) , Problema 62 A partir de los ceros y de los pesos aso- ciados a los polinomios de Legendre, y dado un intervalo [a, b] cualquiera, encontrar los puntos xk, y los pesos wk que hacen exacta hasta el orden 2N − 1 una fórmula de integración numérica sobre el intervalo [a, b] Solución: Para encontrar los puntos x̂k, y los pesos ŵk, hay que hacer un cambio de variable en la integral:Z b a f (x) dx ' NX k=1 ŵkf (x̂k) Hacemos el siguiente cambio de variable: x = (b−a)t+(b+a)2 dx = b−a2 dt este cambio representa la recta que pasa por los puntos -1, 1 para t = a, b, respectivamente. Z b a f (x) dx = Z 1 −1 f µ (b− a) t+ b+ a 2 ¶ b− a 2 dt Z b a f (x) dx ' NX k=1 w̃k b− a 2 f µ (b− a) x̃k + b+ a 2 ¶ de donde se deduce que los cambios a realizar son de la forma x̂k = (b−a)x̃k+(b+a) 2 , ŵk = (b−a) 2 w̃k Problema 63 Utilizar el resultado del problema anterior para calcular de forma exacta la siguiente integralZ 1 0 ¡ x2 − x3 ¢ dx Solución: El resultado de la integral calculada de forma analítica, da el siguiente resultado:R 1 0 ¡ x2 − x3 ¢ dx = 112 = 8. 333 3× 10−2 Aplicando el método de integración numérica: f (x) = ¡ x2 − x3 ¢ R 1 0 ¡ x2 − x3 ¢ dx = 3P k=1 wkf (xk) = = w1f (x1) + w2f (x2) + w3f (x3) = = ¡ 1−0 2 ¢ (w̃0f ¡ x̃0+1 2 ¢ + w̃1f ¡ x̃1+1 2 ¢ + + w̃2f ¡ x̃2+1 2 ¢ ) = = 12(0.55555556 · f ¡ 0.7745966692+1 2 ¢ + +0.8888888889 · f ¡ 0+1 2 ¢ + +0.55555556 · f ¡−0.7745966692+1 2 ¢ ) = = 8. 333 3× 10−2 Problema 64 Calcular de forma exacta la integralZ ∞ −∞ ¡ x3 − x2 ¢ e−x 2 dx utilizando los polinomios de Hermite. Solución: De forma analítica la integral da como resul- tado:R∞ −∞ ¡ x3 − x2 ¢ e−x 2 dx = −12 √ π = −. 886 23 Utilizando el método de integración numérica: f (x) = ¡ x3 − x2 ¢ R∞ −∞ ¡ x3 − x2 ¢ e−x 2 dx = 2P k=1 wkf (xk) = w1f (x1) + w2f (x2) = = 0.8862269255 · f (−0.707106781)+ +0.8862269255 · f (0.707106781) = = −. 886 23 Problema 65 Aproximar, utilizando dos puntos de aproximación, el valor de la integral:Z ∞ −∞ 1 1 + x2 dx Solución:R∞ −∞ 1 1+x2 dx = R∞ −∞ ex 2 1+x2 e −x2dx = π f(x) = e x2 1+x2R∞ −∞ 1 1+x2 dx ' w1f(x1) + w2f(x2) = = 0.8862269255 · f (−0.707106781)+ +0.8862269255 · f (0.707106781) = = 1. 948 2 2 0 Problema 66 Calcular de forma exacta la integralZ ∞ 0 ¡ x3 − x2 ¢ e−xdx utilizando los polinomios de Laguerre. Solución:R∞ 0 ¡ x3 − x2 ¢ e−xdx = 4 R∞ 0 ¡ x3 − x2 ¢ e−xdx = 1P k=0 wkf (xk) = = w0f (x0) + w1f (x1) = = 0.8535533903 · f (0.585786438)+ +0.1464466093 · f (3.414213562) = = 4.0 Problema 67 Calcular una fórmula de aproximación numérica de la integral siguienteZ ∞ a f(x)e−xdx donde a es un número real cualquiera Solución: Para calcular esta integral realizamos un cam- bio de variable R∞ 0 f(t)e−tdx = ½ t = x− a dt = dx ¾ = = R∞ a f(x− a)e−x+adx = ea R∞ a f(x− a)e−x R∞ 0 f(t)e−tdx = NP k=0 w̃kf (x̃k) ea R∞ a f(x− a)e−x = ea NP k=1 wkf (xk − a) Para que estas dos igualdades sean equivalentes, basta hacer: xk = x̃k + a wk = e −aw̃k Problema 68 Aproximar, por el método de Simpson, la integral Z 1 −1 ¡ x3 − x4 ¢ dx utilizando únicamente el valor de la función en los puntos: −1,−12 , 0, 1 2 y 1. Solución:R 1 −1 ¡ x3 − x4 ¢ dx = − 25 = −. 4 Aplicamos el método de Simpson: f(x) = ¡ x3 − x4 ¢ R 1 −1 ¡ x3 − x4 ¢ dx = = R 0 −1 ¡ x3 − x4 ¢ dx+ R 1 0 ¡ x3 − x4 ¢ dx ' ' f(xk+1)+f(xk)+4f xk+1+xk 2 6 (xk+1 − xk) #0 −1 + + f(xk+1)+f(xk)+4f xk+1+xk 2 6 (xk+1 − xk) #1 0 = ' f(0)+f(−1)+4f( −1+0 2 ) 6 (0 + 1)+ + f(1)+f(0)+4f( 1+02 ) 6 (1− 0) = = − 512 = −. 416 67 Problema 69 Deducir la fórmula de integración numérica sobre el rectángulo [−1, 1]x[−1, 1] resultante de aplicar la integración numérica en una variable en los intervalos [−1, 1], y [−1, 1]. Solución:R 1 −1 R 1 −1 F (x, y) dydx = R 1 −1 NP k=1 w̃kF (x̃k, y) dy = = NP k=1 w̃k R 1 −1 F (x̃k, y) dy = NP k=1 w̃k à NP j=1 w̃jF (x̃k, ỹk) ! = NP k,j=1 W̃kF (x̃k, ỹk) , donde W̃k = w̃kw̃j w̃k = R −1 −1 i 6=k(x−x̃i) i6=k(x̃k−x̃i) w̃j = R −1 −1 i 6=k(y−ỹi) i6=k(ỹk−ỹi) y los x̃k e ỹk son los ceros del polinomio de Legendre. Problema 70 Deducir la fórmula de integración numérica sobre un rectángulo [a, b]x[c, d] resultante de aplicar la integración numérica en una variable en los intervalos [a, b], y [c, d]. 2 1 Solución: R b a R d c F (x, y) dydx = R d c NP k=1 w̃kF (x̃k, y) dy = = NP k=1 w̃k R d c F (x̃k, y) dy = NP k=1 w̃k à NP j=1 w̃0jF (x̃k, ỹk) ! = NP k,j=1 w̃kw̃ 0 jF (x̃k, ỹk) , ahora bien, teniendo en cuenta los resultados obtenidos al integrar en una variable tenemos que : x̃k = (b−a)xk+(b+a) 2 w̃k = (b−a) 2 wk ỹk = (d−c)yk+(d+c) 2 w̃0j = (d−c) 2 wk donde wk son los pesos al integrar en una variable en el intervalo [−1, 1]. Problema 71 Calcular de forma exacta la integralZ 1 −1 Z 1 −1 x2y2dxdy utilizando integración numérica. Solución: El resultado de la integral es:R 1 −1 R 1 −1 x 2y2dxdy = 49 = . 444 44 Utilizando la fórmula de integración numérica: P (x) = x2 P (y) = y2R 1 −1 R 1 −1 x 2y2dxdy = = R 1 −1 x 2dx R 1 −1 y 2dy = = 2P k=1 w̃kP (x̃k) 2P k=1 w̃jP (ỹk) = = (w̃1P (x̃1) + w̃2P (x̃2)) · · (w̃1P (ỹ1) + w̃2P (ỹ2)) = = (P (0.5773502692) + P (−0.5773502692)) · · (P (0.5773502692) + P (−0.5773502692)) = = . 666 67 · . 666 67 = . 444 45 Problema 72 Calcular una aproximación numérica de la integral Z ∞ −∞ Z 2 0 x 1 + ey2 dxdy utilizando la evaluación de F (x, y) en 4 puntos. Solución: Si calculamos el resultado de la integral de forma analíta, nos queda,R∞ −∞ R 2 0 x 1+ey2 dxdy = R∞ −∞ 2 1+ey2 dy = = R∞ −∞ 2 1+ey2 dy = 2. 144 3 R 2 0 xdx = 1P k=0 wkP (xk), realizando un cambio de variables, y utilizando el poli- nomio de Legendre de segundo orden, xk = (b−a)x̃k 2 + (b+a) 2 = x̃k + 1, 1. wk = (b−a) 2 w̃k = w̃k, tenemos: P (x) = x R 2 0 xdx = = w1P (x1) + w2P (x2) = = (0.5773502692 + 1) + (−0.5773502692 + 1) = = 2.0 R∞ −∞ 1 1+ey2 dy = R∞ −∞ 1 e−y2+1 e−y 2 dy = 2P k=1 w̃jP (ỹk), por Hermite, 1. P (y) = 1 e−y2+1 1P k=0 wjP (yk) = = w1P (y1) + w2P (y2) = = 0.8862269255 · P (−0.707106781)+ +0.8862269255 · P (0.707106781) = = 1. 103 3 2 2 El resultado de la aproximación numérica es,R∞ −∞ R 2 0 x 1+ey2 dxdy = 2.0 · 1. 103 3 = 2. 206 6 Problema 73 Se considera el triángulo T de vértices (0, 0), (1, 0) y (0, 1). Deducir cual debe ser el punto (x0, y0) y el peso w0 para que la fórmula de integración numérica:Z T F (x, y)dxdy ≈ F (x0,y0)w0 sea exacta para polinomios de grado 1 en x e y. Es decir P (x, y) = ax+ by + c Solución: Calculamos la integral de forma analítica:R 1 0 R 1−x 0 (ax+ by + c) dydx = 16a+ 1 6b+ 1 2c Igualamos el valor de la integral con la fórmula de integración numérica: 1 6a+ 1 6b+ 1 2c = w0 (x0a+ y0b+ c) Calculamos w0, x0 e y0 dando valores a a, b, c a = b = 0, c = 1; 12c = w0c→ w0 = 1 2 a = c = 0, b = 1; 16b = w0y0b→ y0 = 1 3 b = c = 0, a = 1; 16a = w0x0a→ x0 = 1 3 , luego para los valores w0 = 12 , x0 = 1 3 , y0 = 1 3 la fórmula de integración es exacta. Problema 74 Calcular una aproximación numérica de la integral Z Ω x2ydxdy donde Ω es el triángulo de vértices (0, 0), (2, 0) y (0, 2) utilizando 1 punto, 3 puntos, y 4 puntos Solución: El cálculo de la integral de forma analítica nos da: R Ω x2ydxdy = R 2 0 R 2−x 0 ¡ x2y ¢ dydx = 815 = . 533 33 Utilizando las fórmulas de integración numérica: F (x, y) = x2y El área del triángulo→ Area(T ) = 12 ¯̄̄̄ ¯̄ 1 1 12 0 0 0 2 0 ¯̄̄̄ ¯̄ = 2 1. Para 1 punto: R Ω x2ydxdy = = F (23 , 2 3)Area(T ) = = ¡ 2 3 ¢2 2 32 = 16 27 = = . 592 59 2. Para 3 puntos: R Ω x2ydxdy = = 13Area(T ) ¡ F ( 22 , 0) + F ( 2 2 , 2 2) + F (0, 2 2) ¢ = = 23 · 1 = 2 3 = = . 666 67 3. Para 4 puntos: R Ω x2ydxdy = = Area(T )[2548 ¡ F ( 410 , 4 10) + F ( 12 10 , 4 10) + F ( 4 10 , 12 10 ) ¢ − −2748F ( 2 3 , 2 3)] = 8 15 = = . 533 33 ANALISIS NUMERICO MATRICIAL II Problema 75 (4 puntos) Tomar N = 2 y demostrar que la norma k x k2 verifica las propiedades de la defini- ción de norma kxkp = p q |x1|p + |x2|p Solución: En esta demostración vamos a generalizar para cualquier p. Al final particularizamos para p = 2 con el fin de hacer que la demostración sea más sencilla. Las propiedades que debe verificar, para cumplir con la defición de norma, son: 1. kxkp = 0⇐⇒ x = 0; p p |x1|p + |x2|p = 0 =⇒ |x1|p + |x2|p = 0, la suma, en valor absoluto, de elementos distintos de cero da un valor positivo mayor que cero, con lo que para que se cumpla esta condición, se tiene que cumplir que x1 = x2 = 0, c.q.d.. 2 3 2. kλxkp = |λ| kxkp ,∀λ ∈ K y x ∈ E; kλxkp = p p |λx1|p + |λx2|p kλxkp = p p |λ|p |x1|p + |λ|p |x2|p kλxkp = p p |λ|p (|x1|p + |x2|p) kλxkp = |λ| p p |x1|p + |x2|p kλxkp = |λ| kxkp , c.q.d. 3. kx+ ykp ≤ kxkp + kykp ,∀x, y ∈ E; p p |x1 + y1|p + |x2 + y2|p ≤ kxkp + kykp =⇒ =⇒ |x1 + y1|p + |x2 + y2|p ≤ ≤ ³ p p |x1|p + |x2|p + p p |y1|p + |y2|p ´p Para p = 2 tenemos: |x1 + y1|2 + |x2 + y2|2 ≤ ≤ µq |x1|2 + |x2|2 + q |y1|2 + |y2|2 ¶2 =⇒ =⇒ x21 + 2x1y1 + y21 + x22 + 2x2y2 + y22 ≤ ≤ x21+x22+2 p (x21 + x 2 2) p (y21 + y 2 2)+y 2 1+y 2 2 =⇒ =⇒ x1y1 + x2y2 ≤ p (x21 + x 2 2) p (y21 + y 2 2) =⇒ =⇒ x21y21 + 2x1y1x2y2 + x22y22 ≤ ≤ x21y21 + x21y22 + x22y21 + x22y22 =⇒ =⇒ 2x1y1x2y2 ≤ x21y22 + x22y21 =⇒ =⇒ 0 ≤ x21y22 + 2x1y1x2y2 + x22y21 =⇒ =⇒ 0 ≤ (x1y2 + x2y1)2, que siempre es cierto, con lo que queda demostrado. Problema 76 Demostrar que Limp→∞ kxkp = maxi |xi| Solución: Limp→∞ kxkp = Limp→∞ µ p qPN i=1 |xi| p ¶ Extraemos el máximo componente de x, xmax. Limp→∞ µ p qPN i=1 |xi| p ¶ = = Limp→∞ µ p r |xmax|p PN i=1 ³ |xi| |xmax| ´p¶ = = Limp→∞ µ |xmax| p rPN i=1 ³ |xi| |xmax| ´p¶ = = |xmax|Limp→∞ µ p rPN i=1 ³ |xi| |xmax| ´p¶ = = |xmax|Limp→∞ ³PN i=1 ³ |xi| |xmax| ´p´1/p Todos los elementos |xi||xmax| son menores o iguales que 1, con lo que Limp→∞ ³ |xi| |xmax| ´p = ½ 1 si xi = xmax 0 si xi 6= xmax , entonces |xmax|Limp→∞ ³PN i=1 ³ |xi| |xmax| ´p´1/p = = |xmax|Limp→∞ (0 + . . .+ 0 + 1 + . . .+ 1)1/p = = |xmax|, c.q.d. Problema 77 Tomar N = 2, y dibujar el lugar ge- ométrico de los vectores x = (x1, x2) que verifican 1. kxk1 < 1 2. kxk2 < 1 3. kxk∞ < 1 Solución: En las gráficas 1, 2 y 3 se muestran los lugares geométricos de las normas 1, 2 e infinito, respectivamente. 1. kxk1 < 1 =⇒ |x|+ |y| < 1 =⇒ y < 1− x Esta ecuación representa, como borde, una recta de pendiente negativa. Tal y como se ve en la figura 1, el lugar geométrico está contenido en un rombo. 2. kxk2 < 1 =⇒ p (x2 + y2) < 1 =⇒ ¡ x2 + y2 ¢ < 1 Esta es la ecuación de un círculo de radio menor que 1 y centro el origen. En la figura 2 se muestra el lugar geométrico. 24 x y 1 1 1 1 Figure 1: Lugar geométrico de kxk1 3. kxk∞ < 1 =⇒ max(x, y) < 1 Esto representa una recta de valor constante (x, y) menor que 1. En la figura 3 se puede ver el lugar geométrico. Problema 78 Tomar N = 2 y demostrar la siguiente de- sigualdad k x k∞≤k x k2≤k x k1 Solución: Esta desigualdad es equivalente a lo siguiente: max(|x1| , |x2|) ≤ p (x21 + x 2 2) ≤ |x1|+ |x2| 1. max(|x1| , |x2|) ≤ p (x21 + x 2 2)⇐⇒ ⇐⇒ |xmax| ≤ p (x21 + x 2 2)⇐⇒ ⇐⇒ x2max ≤ x21 + x22 Esta desigualdad siempre es cierta ya que xmax es o bien x1 o bien x2. 2. p (x21 + x 2 2) ≤ |x1|+ |x2|⇐⇒ ⇐⇒ ¡ x21 + x 2 2 ¢ ≤ (|x1|+ |x2|)2 ⇐⇒ ⇐⇒ x21 + x22 ≤ |x1| 2 + 2 |x1| |x2|+ |x2|2 ⇐⇒ ⇐⇒ x21 + x22 ≤ x21 + 2 |x1| |x2|+ x22 ⇐⇒ ⇐⇒ 0 ≤ 2 |x1| |x2| x y 1 1 1 1 Figure 2: Lugar geométrico de kxk2 Esto siempre es cierto ya que el producto de valores positivos siempre es positivo (o igual a cero si algún xi es cero). 3. max(|x1| , |x2|) ≤ |x1| + |x2| . Es trivial (propiedad transitiva). De estas demostraciones se deduce que las distintas normas coinciden cuando el vector x descansa sobre uno de los ejes de coordenadas. Problema 79 Demostrar que si A,B son dos matrices de dimensión NxN, entonces para cualquier norma de matri- ces subordinada a una norma vectorial se verifica k AB k≤k A k · k B k Solución: supx kABxk kxk = supx kABxk kxk kBxk kBxk , supx kABxk kBxk kBxk kxk ≤ supx kBxk kxk · supx kABxk kBxk supx kABxk kBxk kBxk kxk ≤ kBk · kAk , entonces supx kABxk kxk ≤ kBk · kAk , c.q.d. Problema 80 Demostrar que los autovalores de A son los ceros del polinomio característico P (λ). 2 5 x y 1 1 1 1 Figure 3: Lugar geométrico de kxk∞ Solución: Definición de autovalor de una matriz A: xi 6= 0 ∈ E, λi ∈ CÁAxi = λixi Axi = λixi =⇒ (A− λiId)xi = 0 como xi 6= 0, entonces |A− λiId| = 0 =⇒ P (λ) = 0, c.q.d. Problema 81 Calcular los autovectores de la matriz⎛⎝ 1 1 01 1 0 0 0 2 ⎞⎠ y determinar una base ortonormal de R3 de autovectores de A. Solución: Calculamos los autovalores de A : |A− λiId| = 0¯̄̄̄ ¯̄ 1− λ 1 01 1− λ 0 0 0 2− λ ¯̄̄̄ ¯̄ = ((1− λ)2 − 1)(2− λ) = 0 λ1 = 0 λ2 = 2 λ3 = 2 Calculamos los autovectores de A : 1. λ1 = 0⎛⎝ 1 1 01 1 0 0 0 2 ⎞⎠⎛⎝ x1x2 x3 ⎞⎠ = ⎛⎝ 00 0 ⎞⎠ ⎧⎨⎩ x1 + x2 = 0x1 + x2 = 0 2x3 = 0 ⎧⎨⎩ x1 = −x2 x3 = 0 x̄1 = ⎛⎝ 1√2−1√ 2 0 ⎞⎠ 2. λ2, λ3 = 2⎛⎝ −1 1 01 −1 0 0 0 0 ⎞⎠⎛⎝ x1x2 x3 ⎞⎠ = ⎛⎝ 00 0 ⎞⎠ ⎧⎨⎩ −x1 + x2 = 0x1 − x2 = 0 x3 libre ⎧⎨⎩ x1 = x2 x3 libre x̄2 = ⎛⎝ 1√21√ 2 0 ⎞⎠ , x̄3 = ⎛⎝ 00 1 ⎞⎠ La matriz, B = ⎛⎝ 1√2 1√2 0−1√ 2 1√ 2 0 0 0 1 ⎞⎠ contiene los autovectores de A que forman una base ortog- onal en R3. Problema 82 Calcular las normas 2, 1 e infinito de la matriz A = µ 1 0 1 1 ¶ Solución: 1. kAk2 = p ρ (ATA) kAk2 = s ρ µ 2 1 1 1 ¶ = q 3+1 √ 5 2 = 1. 618 ¯̄̄̄ 2− λ 1 1 1− λ ¯̄̄̄ = 1− 3λ+ λ2 = 0, λ = 32 ± 1 2 √ 5 2. kAk1 = maxj ³P2 i=1 |aij | ´ kAk1 = max(1, 2) = 2 3. kAk∞ = maxi ³P2 j=1 |aij | ´ = kAk∞ = max(2, 1) = 2 2 6 Problema 83 Demostrar la siguiente igualdad: ρ(AtA) = ρ(AAt) Solución: Veamos que los polinomios característicos coinciden : |AtA− λiId| = |At|−1 |AtA− λiId| |At| = = ¯̄̄ (At) −1 AtAAt − λi (At)−1 IdAt ¯̄̄ = = ¯̄̄ AAt − λi (At)−1At ¯̄̄ = = |AAt − λiId| Problema 84 Demostrar que si los autovectores de una matriz A de dimensión NxN forman una base ortonormal de RN , entonces para la norma 2 se cumple: χ(A) = kAk2 · °°A−1°° 2 = maxi{|λi|} mini{|λi|} Solución: Al ser una base de autovectores ortonormal, la norma kAk2 = ρ(A) = maxi {|λi|} Los autovalores de A−1 vienen dados por: Axi = λixi =⇒ =⇒ A−1Axi = A−1λixi =⇒ =⇒ 1λixi = A −1xi =⇒ =⇒ A−1xi = λ0ixi, donde λ0i = 1 λi , es decir, los autovalores de A−1 son los inversos de los de A y sus autovectores son los mismos, luego la norma de °°A−1°° 2 = ρ(A−1) °°A−1°°2 = max i ©¯̄ λ0i ¯̄ª = max i ½¯̄̄̄ 1 λi ¯̄̄̄¾ = 1 mini {|λi|} , entonces χ(A) = kAk2 · °°A−1°° 2 χ(A) = maxi {|λi|} · 1mini{|λi|} χ(A) = maxi{|λi|}mini{|λi|} , c.q.d. Problema 85 Calcular el condicionamiento para la norma 2, de las siguientes matrices: 1. A = ⎛⎝ 2 2 −22 1 1 −2 1 1 ⎞⎠ 2. A = ⎛⎝ 2 −1 0−1 2 −1 0 −1 2 ⎞⎠ Solución: El condicionamiento de una matriz χ(A) = kAk2 · °°A−1°° 2 . Calculamos los autovalores de ambas ma- trices: 1. ¯̄̄̄ ¯̄ 2− λ 2 −22 1− λ 1 −2 1 1− λ ¯̄̄̄ ¯̄ = (2− λ) £(1− λ)2 − 1¤ −2 [2(1− λ) + 2] −2 [2 + 2(1− λ)] = (2− λ) £ λ2 − 2λ ¤ − 8 (2− λ) = (2− λ) £ λ2 − 2λ+ 8 ¤ de donde obtenemos λ1 = 2 λ2 = −2 λ3 = 4 Esta matriz es simétrica, luego posee una base orto- normal 3 de autovectores, con lo que el condi- cionamiento de A se puede calcular como: χ(A) = kAk2 · °°A−1°° 2 = maxi{| λi |} mini{| λi |} = 4 2 = 2 2. ¯̄̄̄ ¯̄ 2− λ −1 0−1 2− λ −1 0 −1 2− λ ¯̄̄̄ ¯̄ = 4− 10λ+ 6λ2 − λ3 = 0 λ1 = 2 λ1 = 2 + √ 2 λ1 = 2− √ 2 También es una matriz simétrica, con lo que sus au- tovectores forman una base ortonormal y su condi- cionamiento es: χ(A) = kAk2· °°A−1°° 2 = maxi{| λi |} mini{| λi |} = 2 + √ 2 2− √ 2 = 3+2 √ 2 Problema 86 Sean las matrices A,R. Demostrar que la matriz A, y la matriz B = R−1AR poseen los mismos autovalores Solución: Bxi = λixi =⇒ =⇒ ¡ R−1AR ¢ xi = λixi =⇒ =⇒ RR−1ARxi = Rλixi =⇒ 3Vectores ortonormales: dos vectores son ortonormales si cumplen lo siguiente, xTi xj = 0 si i 6= j 1 si i = j 2 7 =⇒ ARxi = λiRxi =⇒ =⇒ Ayi = λiyi, de donde se deduce que los autovalores son los mismos y los autovectores están relacionados por la siguiente igualdad: yi = Rxi, c.q.d. Problema 87 Se considera la matriz A = µ 1 1 1 1 ¶ calcular el ángulo α tal que la matriz R = µ cosα sinα − sinα cosα ¶ verifique que la matriz B = R−1AR sea diagonal. Solución: Realizamos el cálculo de la matriz B : B = R−1AR = = µ cosα − sinα sinα cosα ¶µ 1 1 1 1 ¶µ cosα sinα − sinα cosα ¶ = = µ −2 cosα sinα+ 1 2 cos2 α− 1 2 cos2 α− 1 2 cosα sinα+ 1 ¶ = = µ b1 0 0 b2 ¶ Se tiene que cumplir que los elementos que están fuera de la diagonal sean iguales a cero, 2 cos2 α− 1 = 0 cosα = ± r 1 2 De esta igualdad se obtiene el valor del ángulo α : α = π 4 , α = 3π 4 La matriz de rotación queda como sigue: R1 = µ cos π4 sin π 4 − sin π4 cos π 4 ¶ = µ 1 2 √ 2 12 √ 2 −12 √ 2 12 √ 2 ¶ R2 = µ cos 3π4 sin 3π 4 − sin 3π4 cos 3π 4 ¶ = µ −12 √ 2 12 √ 2 −12 √ 2 −12 √ 2 ¶ Calculamos los elementos de la diagonal: b1 = −2 cosα sinα+ 1 b1 = 0, b1 = 2 b2 = 2 cosα sinα+ 1 b2 = 2, b2 = 0, luego las soluciones posibles son: B1 = µ 0 0 0 2 ¶ , B2 = µ 2 0 0 0 ¶ Problema 88 Demostrar las siguientes igualdades trigonométricas tan(α) = − cot(2α) + sign(cot(2α)) q 1 + cot2(2α) donde α ∈ ¡ −π4 , π 4 ¢ , sign(x) = 1 si x ≥ 0 y sign(x) = −1 si x < 0, cosα = 1p 1 + tan2(α) sinα = tan(α) cosα cot(2α) = − tan(α) + sin(2α) 2 sin2(α) Solución: 1. cot(2α) = cos(2α)sin(2α) = cos2(α)−sin2(α) 2 sin(α) cos(α) = 1−tan2(α) 2 tan(α) 2 tan(α) cot(2α) = 1− tan2(α) realizando el cambio de variable x = tan(α), tenemos x2 + 2 cot(2α)x− 1 = 0 x = tan(α) = −2 cot(2α)± √ 4 cot2(2α)+4 2 = = − cot(2α)± p 1 + cot2(2α) tan(α) = ½ − cot(2α) + p 1 + cot2(2α) si α ≥ 0 − cot(2α)− p 1 + cot2(2α) si α < 0 El segundo término es siempre mayor que el primero, con lo que es éste el que va a determinar el signo de la ecuación. Como sign (tan(α)) = sign (cot(α)) , podemos expre- sar la anterior igualdad de la siguiente forma: tan(α) = − cot(2α) + sign(cot(2α)) q 1 + cot2(2α) 2. cosα = 1√ 1+tan2(α) = 1 1+ sin2(α) cos2(α) = = 1 cos2(α)+sin2(α) cos2(α) p cos2(α) = cosα 3. sinα = tan(α) cosα = = sin(α)cos(α) cosα = sinα 2 8 4. cot(2α) = − tan(α)+sin(2α) 2 sin2(α) = = − sin(α) cos(α) +2 sin(α) cos(α) 2 sin2(α) = = − sin(α)+2 sin(α) cos(α) cos(α) cos(α) 2 sin2(α) = = sin(α)(−1+2 cos(α) cos(α)) 2 sin2(α) cos(α) = (2 cos2(α)−1) 2 sin(α) cos(α) = = cos(2α)sin(2α) = cot(2α) Problema 89 Dentro del método de Jacobi para el cálculo de autovalores demostrar las igualdades a0pq = 0 a0pp = app − tan(α)apq a0qq = aqq + tan(α)apq a0pj = apj cosα− aqj sinα j 6= p, q a0qj = apj sinα+ aqj cosα j 6= p, q Solución: En el método de Jacobi se persigue construir una matriz diagonal a partir de una matriz A cualquiera, aplicándole transformaciones de la forma B = R−1AR. Según se ha demostrado en problemas anteriores, los autovalores de B y de A coinciden, con lo que si se consigue encontrar la matriz R que cumpla con la ecuación anterior, entonces habremos encontrado los autovalores de A. La matriz R es una matriz de rotación y se calcula el ángulo, α, de la misma, transformando los valores de A que están fuera de la diagonal en ceros. Vamos a expresar las matrices de la siguiente manera: A = ⎛⎜⎜⎜⎜⎝ a11 ap1 ai1 aq1 an1 ap1 app apj apq apn ai1 apj aij aqj ani aq1 apq aqj aqq aqn an1 apn ani aqn ann ⎞⎟⎟⎟⎟⎠ Rpq (α) = ⎛⎜⎜⎜⎜⎝ 1 0 0 0 0 0 cosα . sinα 0 0 . 1 . 0 0 − sinα . cosα 0 0 0 0 0 1 ⎞⎟⎟⎟⎟⎠ B = R−1AR = = ⎛⎜⎜⎜⎜⎝ 1 0 0 0 0 0 cosα 0 − sinα 0 0 0 1 0 0 0 sinα 0 cosα 0 0 0 0 0 1 ⎞⎟⎟⎟⎟⎠ · · ⎛⎜⎜⎜⎜⎝ a11 ap1 ai1 aq1 an1 ap1 app apj apq apn ai1 apj aij aqj ani aq1 apq aqj aqq aqn an1 apn ani aqn ann ⎞⎟⎟⎟⎟⎠ · · ⎛⎜⎜⎜⎜⎝ 1 0 0 0 0 0 cosα 0 sinα 0 0 0 1 0 0 0 − sinα 0 cosα 0 0 0 0 0 1 ⎞⎟⎟⎟⎟⎠ = = ⎛⎜⎜⎜⎜⎝ a11 ap1 cosα− aq1 sinα ap1 cosα− aq1 sinα app cos2 α+ aqq sin2 α− apq sin 2α ai1 apj cosα− aqj sinα ap1 sinα+ aq1 cosα (app−aqq) 2 sin 2α+ apq cos 2α an1 apn cosα− aqn sinα ai1 apj cosα− aqj sinα aij apj sinα+ aqj cosα ani ap1 sinα+ aq1 cosα an1 (app−aqq) 2 sin 2α+ apq cos 2α apn cosα− aqn sinα apj sinα+ aqj cosα ani app sin 2 α+ aqq cos 2 α+ apq sin 2α apn sinα+ aqn cosα apn sinα+ aqn cosα ann ⎞⎟⎟⎟⎟⎠ De esta matriz se deducen las siguientes igualdades: a0pq = (app−aqq) 2 sin 2α+ apq cos 2α a0pp = app cos 2 α+ aqq sin 2 α− apq sin 2α a0qq = app sin 2 α+ aqq cos 2 α+ apq sin 2α a0pj = apj cosα− aqj sinα j 6= p, q aqj = apj sinα+ aqj cosα j 6= p, q En donde se iguala a0pq a cero para calcular el ángulo de rotación: a0pq = 0 = (app−aqq) 2 sin 2α+ apq cos 2α tan(2α) = 2apq (aqq − app) aqq = app + 2apq tan(2α) Las dos últimas igualdades se obtienen directamente de la matriz final. Para obtener a0pp y a 0 qq, se opera de la siguiente manera: 1. a0pp = app cos 2 α+ aqq sin 2 α− apq sin 2α = = app cos 2 α+ ³ app + 2apq tan(2α) ´ sin2 α− −apq sin 2α = app cos2 α+ 2 9 + ³ app sin(2α)+2apq cos(2α) 2 sinα cosα ´ sin2 α− apq sin 2α = = app cos 2 α+ + ³ app2 sinα cosα+2apq cos 2 α−2apq sin2 α 2 cosα ´ sinα− −2apq sinα cosα = app cos2 α+ app sin2 α+ +apq cosα sinα− apq tanα+ apq sinα cosα− −2apq sinα cosα = app − apq tanα 2. a0qq = app sin 2 α+ aqq cos 2 α+ apq sin 2α = = ³ aqq − 2apqtan(2α) ´ sin2 α+ aqq cos 2 α+ +apq sin 2α = = ³ aqq2 sinα cosα−2apq cos2 α+2apq sin2 α 2 sinα cosα ´ sin2 α+ +aqq cos 2 α+ apq sin 2α = (aqq sinα− apq cosα+ + apq cosα−apq cosα) sinα+aqq cos2 α+apq sin 2α = = aqq sin 2 α+ aqq cos 2 α− apq cosα sinα+ +apq tanα− apq cosα sinα+ 2apq cosα sinα = = aqq + tan(α)apq Problema 90 Utilizar el método de Jacobi para aproxi- mar los autovalores y autovectores de la matriz: A = ⎛⎝ 2 0 10 1 0 1 0 1 ⎞⎠ Solución: R (α) = ⎛⎝ cosα 0 sinα0 1 0 − sinα 0 cosα ⎞⎠ (app−aqq) 2 sin 2α+ apq cos 2α = 0 tan(2α) = 2apq (aqq − app) tan(2α) = 2 (1− 2) α = arctan (−2) 2 α = −1 2 arctan 2 = −. 553 57 a11 = 2− tan(α) a11 = 2 + tan µ 1 2 arctan 2 ¶ = 2. 618 a33 = 1 + tan(α) a33 = 1− tan µ 1 2 arctan 2 ¶ = . 381 97 a21 = a32 = 0 B = R−1AR = ⎛⎝ 2. 618 0 00 1 0 0 0 0. 381 97 ⎞⎠ Los autovalores son los elementos de la diagonal (2. 618, 1, 0. 381 97) . Como α = −. 553 57, la matriz R (α) = ⎛⎝ cosα 0 sinα0 1 0 − sinα 0 cosα ⎞⎠ = ⎛⎝ 0.850 65 0 −0.525 730 1 0 0.525 73 0 0.850 65 ⎞⎠ por tanto, en este caso, como con una única matriz de rotación conseguimos transformar A en una matriz diago- nal, tendremos que los autovectores de A son simplemente los vectores columnasde R(α). Es decir el autovector del autovalor 2. 618 es ( 0.850 65 0 0.525 73 ), el autovector del auto- valor 1 es ( 0 1 0 ), el autovector del autovalor 0. 381 97 es ( −0.525 73 0 0.850 65 ). Problema 91 Aplicar el método de la potencia para aproximar el autovalor máximo, y el autovector asociado, de las siguientes matrices, dando 3 pasos en el método, hasta calcular u4 y partiendo de u1 = (1, 1) A = µ 2 1 0 1 ¶ A = µ −3 0 1 1 ¶ Solución: En este problema vamos a utilizar la norma euclídea aunque cualquier otra norma también sería válida. La norma infinito, por ejemplo, simplificaría los cálculos ya que es inmediato obtener el máximo de un vector. 1. A = µ 2 1 0 1 ¶ u2 = A u 1 ku1k = 3 0 = µ 2 1 0 1 ¶Ã 1√ 2 1√ 2 ! = = µ 3 2 √ 2 1 2 √ 2 ¶ , °°u2°° = √5 = 2. 236 1 u3 = A u 2 ku2k = = µ 2 1 0 1 ¶Ã 3 2 √ 5 √ 2 1 2 √ 5 √ 2 ! = = µ 7 10 √ 5 √ 2 1 10 √ 5 √ 2 ¶ , °°u3°° = √5 = 2. 236 1 u4 = A u 3 ku3k = = µ 2 1 0 1 ¶µ 7 10 √ 2 1 10 √ 2 ¶ = = µ 3 2 √ 2 1 10 √ 2 ¶ , °°u4°° = 1 5 √ 113 = 2. 126 El autovalor máximo aproximado es λ = 2. 126 y su autovector asociado es: xλ = µ 15 226 √ 113 √ 2 1 226 √ 113 √ 2 ¶ = µ . 997 79 6. 651 9× 10−2 ¶ 2. A = µ −3 0 1 1 ¶ u2 = A u 1 ku1k = = µ −3 0 1 1 ¶Ã 1√ 2 1√ 2 ! = = µ − 32 √ 2√ 2 ¶ , °°u2°° = 1 2 √ 26 = 2. 549 5 u3 = A u 2 ku2k = = µ −3 0 1 1 ¶µ − 326 √ 26 √ 2 1 13 √ 26 √ 2 ¶ = = µ 9 26 √ 26 √ 2 − 126 √ 26 √ 2 ¶ , °°u3°° = 1 13 √ 1066 = 2. 511 5 u4 = A u 3 ku3k = = µ −3 0 1 1 ¶µ 9 2132 √ 1066 √ 26 √ 2 − 12132 √ 1066 √ 26 √ 2 ¶ = = µ − 272132 √ 1066 √ 26 √ 2 2 533 √ 1066 √ 26 √ 2 ¶ , °°u4°° = 1 82 √ 65 026 = 3. 109 8 El autovalor máximo aproximado es λ = −3. 109 8, con signo negativo ya que sign ¡ u4, u3 ®¢ = −1 y su autovector asociado es xλ = à − 27·82 2132 √ 65 026 √ 1066 √ 26 √ 2 2·82 533 √ 65 026 √ 1066 √ 26 √ 2 ! = µ −. 958 8 . 284 09 ¶ , con signo positivo ya que ¡ sign ¡ u4, u3 ®¢¢n = (−1)4 = 1. Problema 92 Calcular el autovalor mayor y el autovec- tor correspondiente de la matriz µ 2 −1 −1 1 ¶ utilizando el método de la potencia, dando 2 iteraciones del método a partir de u1 = (1, 1) y tomando como norma kuk = maxi |ui| Solución: 1. ku1k = 1→ u2 = µ 2 −1 −1 1 ¶µ 1 1 ¶ = µ 1 0 ¶ ku2k = 1→ u3 = µ 2 −1 −1 1 ¶µ 1 0 ¶ = µ 2 −1 ¶ Producto escalar (u2, u3) = 2 > 0. → autovalor máx- imo = ku3k = 2 Autovector asociado normalizado u3ku3k = µ 1 −1/2 ¶ 3 1 Problema 93 Utilizar el método de la potencia inversa para aproximar el autovalor más pequeño de la matriz A = µ −2 1 0 3 ¶ llegar hasta u3 partiendo de u = (1, 1). Solución: Aun = u n−1 kun−1kµ −2 1 0 3 ¶ u2 = à 1√ 2 1√ 2 ! u2 = µ −16 √ 2 1 6 √ 2 ¶ , °°u2°° = 1 3 = . 333 33 µ −2 1 0 3 ¶ u3 = µ −12 √ 2 1 2 √ 2 ¶ u3 = µ 1 3 √ 2 1 6 √ 2 ¶ , °°u3°° = 1 6 √ 10 = . 527 05 °°u3°° es el autovalor máximo de A−1, con lo que el autovalor mínimo de A es λmin = −1ku3k = − 6 10 √ 10 = −1. 897 4, con signo negativo ya que sign ¡ u3, u2 ®¢ = −1. Problema 94 Calcular el autovalor y autovector más cer- cano a 2 de la matriz⎛⎝ 0 −1 00 3 −1 0 0 −1 ⎞⎠ para ello calcular dos iteraciones del método de la potencia inversa partiendo de u1 = (1, 1, 1). Solución: A0 = A− 2Id = ⎛⎝ −2 −1 00 1 −1 0 0 −3 ⎞⎠ Vamos a utilizar la norma infinito con el fin de sim- plificar los cálculos. A0un = u n−1 kun−1k⎛⎝ −2 −1 00 1 −1 0 0 −3 ⎞⎠u2 = ⎛⎝ 11 1 ⎞⎠ u2 = ⎛⎝ −562 3 −13 ⎞⎠ ,°°u2°° = 5 6 = . 833 33 ⎛⎝ −2 −1 00 1 −1 0 0 −3 ⎞⎠u3 = 65 ⎛⎝ −562 3 −13 ⎞⎠ u3 = ⎛⎝ 13014 15 2 15 ⎞⎠ ,°°u2°° = 14 15 El autovalor máximo de (A− 2Id)−1 es λmax = 1415 con signo positivo (sign ¡ u3, u2 ®¢ = 1) (A− 2Id)−1 x̄ = λmaxx̄ Para calcular el autovalor más cercano a 2, realizamos las siguientes operaciones: 1 λmax x̄ = (A− 2Id) x̄ µ A− 2Id− 1 λmax Id ¶ x̄ = 0µ A− µ 2 + 1 λmax ¶ Id ¶ x̄ = (A− λproxId) x̄ = 0 λprox = µ 2 + 1 λmax ¶ = µ 2 + 1 14 15 ¶ = 43 14 λprox = 3. 071 4 Problema 95 Calcular 3 iteraciones del método de Ja- cobi para resolver el sistema⎛⎝ 1 −1 0−1 2 0 0 −1 3 ⎞⎠⎛⎝ xy z ⎞⎠ = ⎛⎝ −13 1 ⎞⎠ partiendo de u1 = (0, 0, 0) Solución: Despejamos la diagonal para plantear el método iterativo : xn = yn−1 + 1 yn = xn−1 + 3 2 zn = yn−1 + 1 3 haciendo iteraciones obtenemos 1. u2 = ⎛⎝ −13 2 1 3 ⎞⎠ 2. u3 = ⎛⎝ 121 5 6 ⎞⎠ 3 2 3. u4 = ⎛⎝ 07 4 2 3 ⎞⎠ Problema 96 Calcular una base ortogonal de autovec- tores de la matriz ⎛⎝ 1 0 10 2 0 1 0 1 ⎞⎠, Solución: 1. Autovectores y autovalores: ⎧⎨⎩ ⎛⎝ −10 1 ⎞⎠⎫⎬⎭ ↔ 0, ⎧⎨⎩ ⎛⎝ 10 1 ⎞⎠ , ⎛⎝ 01 0 ⎞⎠⎫⎬⎭↔ 2 Problema 97 Calcular 3 iteraciones del método de Gauss-Seidel para resolver el sistema⎛⎝ 1 −1 0−1 2 0 0 −1 3 ⎞⎠⎛⎝ xy z ⎞⎠ = ⎛⎝ −13 1 ⎞⎠ partiendo de u1 = (0, 0, 0) Solución: Despejamos la diagonal para plantear el método iterativo, teniendo en cuenta además que vamos actualizando los valores según los cálculamos: xn = yn−1 + 1 yn = xn + 3 2 zn = yn + 1 3 haciendo iteraciones partiendo de (0, 0, 0) : 1. x1 = −1 y1 = 3−1 2 = 1 z1 = 1+1 3 = 2 3 2. x2 = 0 y2 = 3 2 z3 = 5 6 3. x3 = −1 + 32 = 1 2 y3 = 3+ 12 2 = 7 4 z3 = 1+ 74 3 = 11 12 Problema 98 Una variante del método de Gauss-Seidel es tomar M = (D + U)−1 (−L), y c = (D + U)−1 b. indicar en este caso que diferencias de implementación habría con respecto al caso anterior. Solución: El método es igual que en el problema anterior, excepto que en este caso los cálculos se realizarían de abajo para arriba, es decir, primero se calcularía z, se sustituiría su valor en la ecuación de y y, por último, estos dos valores se sustituirían en la primera ecuación. Problema 99 Calcular 3 iteraciones del método de rela- jación para resolver el sistema⎛⎝ 1 −1 0−1 2 0 0 −1 3 ⎞⎠⎛⎝ xy z ⎞⎠ = ⎛⎝ −13 1 ⎞⎠ , partiendo de u1 = (0, 0, 0). Calcular previamente el parámetro de relajación óptimo Solución: x− y = −1 −x+ 2y = 3 −y + 3z = 1 Cálculo del wopt: Al ser A tridiagonal, el wopt se puede calcular como wopt = 2 1 + q 1− ρ (MJ)2 MJ es la matriz del método de Jacobi que se obtiene despejando la diagonal en el sistema x = y − 1 y = x2 + 3 2 z = y3 + 1 3 = ⎛⎝ 0 1 01 2 0 0 0 13 0 ⎞⎠⎛⎝ xy z ⎞⎠+ ⎛⎝ −13 2 1 3 ⎞⎠ Los autovalores de MJ : ³ 0, 1√ 2 ,− 1√ 2 ´ , luego ρ (MJ) 2 = 12 wopt = 2 1+ √ 1−ρ(MJ)2 = 2 1+ √ 1− 12 wopt = 2 1+ 12 √ 2 = 4 2+ √ 2 wopt = 1. 171 6 Iteraciones del sistema: xn = w (yn−1 − 1) + (1− w)xn−1 yn = w 3+xn 2 + (1− w) yn−1 zn = w 1+yn 3 + (1− w) zn−1 u2 = ⎛⎝ −ww 3−w2 w 1+1. 071 13 ⎞⎠ = ⎛⎝ −1. 171 61. 071 1 . 808 83 ⎞⎠ 3 3 u3 = ⎛⎝ w (1. 071 1− 1)− (1− w) 1. 171 6w 3+. 284 352 + (1− w) 1. 071 1 w 1+1. 740 23 + (1− w) . 808 83 ⎞⎠ = = ⎛⎝ . 284 351. 740 2 . 931 34 ⎞⎠ u4 = ⎛⎝ w (1. 740 2− 1) + (1− w) . 284 35w 3+. 818 422 + (1− w) 1. 740 2 w 1+1. 938 23 + (1− w) . 931 34 ⎞⎠ = = ⎛⎝ . 818 421. 938 2 . 987 65 ⎞⎠ Problema 100 Demostrar que si una matriz A verifica que por filas o columnas su suma es siempre igual a 0, en- tonces el determinante de A es cero, y por tanto el sistema asociado a A no tiene solución. Solución: Si |A| = 0, entonces la matriz A no es invertible y el sistema no tiene solución. 1. Vamos a demostrar que si la suma por filas de A es igual a cero, entonces |A| = 0 :Pn j aij = 0, esto es equivalente a lo siguiente: A ⎛⎜⎜⎜⎝ 1 1 ... 1 ⎞⎟⎟⎟⎠ = ⎛⎜⎜⎜⎝ 0 0 ... 0 ⎞⎟⎟⎟⎠ = 0 · ⎛⎜⎜⎜⎝ 1 1 ... 1 ⎞⎟⎟⎟⎠ Esto significa que la matriz A posee un autovalor igual a cero (λ = 0). El determinante de una matriz es igual al producto de sus autovalores: |A| = nY i=1 λi = λ1 · . . . 0 . . . · λn = 0 2. Para demostrar que |A| = 0 cuando la suma por columnas es cero, basta saber que |A| = ¯̄ AT ¯̄ y aplicar el argumento anterior a AT Problema 101 Dado un sistema iterativo un =Mun−1 + c Demostrar que aunque el radio espectral de M sea mayor que 1, si u1 y c son combinaciones lineales de autovectores deM correspondientes a autovalores de módulo menor que 1, entonces el método converge. Solución: Sean xi los autovectores deM correspondientes a autovalores menores
Compartir