Problemas Resueltos de Analisis Numericos I

Problemas de Análisis Numérico 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 Ta…ra 35017 Las Palmas, España T‡: 45.87.10/08 Email: {maleman/lalvarez/jsanchez}@dis.ulpgc.es

Contents 1 INTRODUCCION.

prisas de última hora que suelen asaltar a los estudiantes cuando se acercan los exámenes, puesto que el esfuerzo de re‡exión que requieren precisa de un trabajo diario y continuado, difícilmente compatible con las prisas de última 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 reproducir esas técnicas, de forma absolutamente mecánica, en problemas análogos. El problema de esta aptitud, es que aunque a corto plazo puede dar lugar a resultados positivos, aprobando asignaturas con un conocimiento mínimo e insu…ciente, a la larga, tiene efectos catastró…cos sobre la formación del alumno, a través de una disminución importante de la capacidad de razonamiento y del sentido crítico.

1

2 ARITMETICAS DE PRECISION FINITA Y FUENTES DE ERRORES NUMERICOS. 1 3 CALCULO DE CEROS DE UNA FUNCION 4 4 INTERPOLACION DE FUNCIONES I

6

5 ANALISIS NUMERICO MATRICIAL I

9

6 DIFERENCIACION NUMERICA

E

INTEGRACION

14

7 ANALISIS NUMERICO MATRICIAL II

23

8 INTERPOLACION DE FUNCIONES II

37

ARITMETICAS DE PRECISION FINITA Y FUENTES DE ERRORES NUMERICOS. Problema 1 Demostrar que al representar el número real 0:1 como 1 X an 0:1 = 2e 2n n=1

INTRODUCCION. El presente documento es el libro de problemas donde se encuentran resueltos todos los problemas presentes en el libro de Análisis Numérico publicado por los mismos autores. Nunca se insistirá lo su…ciente 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 re‡exionar 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 su…ciente esfuerzo, que de un problema del cual se mira directamente la solución sin ninguna fase de re‡exió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 re‡exionado su…cientemente. De todo ello se deduce que el estudio correcto de los problemas de una asignatura va reñido con las

el número de elementos no-nulos an es in…nito.

Solución: Supongamos que para algún t …nito y e entero se tiene: t X an e 0:1 = 2 2n n=1 despejando en esta igualdad obtenemos t¡e

2

= 10

t X

an 2t¡n

n=1

P ahora bien, como el número m = tn=1 an 2t¡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

Solución: Los valores posibles positivos se representan en la siguiente tabla

1 X an

n=1

2n

Solución: En primer lugar tenemos que encontrar un entero e tal que

e=0 e=1

1 · 0:0 703 125 ¢ 2¡e < 1 2

e=2

para e = ¡3 obtenemos

los valores negativos son los mismos cambiados de signo. Simpli…cando las fracciones nos queda

3

0:0 703 125 ¢ 2 = 0: 562 5

ahora tenemos que escribir el número 0:5625 como 0:5625 =

e = ¡1 e=0 e=1 e=2

1 X

1 an + 2 n=2 2n

los an se calculan de la siguiente forma

0:5625 < 0:5625 < 0:5625 =

0:0 703 125 = 2¡3

µ

1 1 + 2 24

2 1.75 1.5 1.25 y 1 0.75 0.5 0.25 0 -0.25 -0.5 -0.75 -1 -1.25 -1.5 -1.75 -2

¶

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; ::::)

xmax

1 2

¡

1

2t+1 1 2

0.250.5 0.751

1.251.5 1.752

2.252.5 2.753

3.253.5 3.754

Problema 5 Dada una aritmética de precisión …nita 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 inferior.

Solución: Los valores positivos mínimo y máximo son = 2emin ¡1 t X 1 = 2emax = 2emax n 2 n=1

0

x

Problema 3 Calcular los valores positivos mínimo y máximo que puede tomar un número real en una aritmética de precisión …nita en función de t; emin y emax :

xmin

0:25; 0:3125; 0:375; 0:437 5 0:5; 0:625; 0:75; 0:875 1; 1:25; 1:5; 1:75 2; 2:5; 3; 3:5

Si representamos los números positivos sobre una recta obtenemos

1 1 + 2 = 0:75 ) a2 = 0 2 2 1 1 + 3 = 0:625 ) a3 = 0 2 2 1 1 + 4 = 0:5625 ) a4 = 1 2 2

por tanto

1 1 1 1 1 1 1 1 ; + ; + ; + + 22 22 24 22 23 22 23 24 1 1 1 1 1 1 1 1 ; + ; + ; + + 2 2 23 2 22 2 22 23 1 1 1 1 1; 1 + 2 ; 1 + ; 1 + + 2 2 2 2 2 1 1 2; 2 + ; 2 + 1; 2 + 1 + 2 2

e = ¡1

Solución: El número 1 en una aritmética de precisión …nita se escribe como µ ¶ 1 1=2 2

¶ µ 1 = 2emax 1 ¡ t 2

el número inmediato superior a 1 en la aritmética es µ ¶ 1 1 1 2 + t = 1 + t¡1 2 2 2

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.

y el número inmediato inferior a 1 viene dado por 1 1 + :: + t = 2 2

2

1 2

1 ¡ 2t+1 1 =1¡ t 1 2 1¡ 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:

2

1 23

+

1 24

+

1 23

+

1 24

¢ +

1 25

Solución: Que los números pertenecen a la aritmética signi…ca que existe un conjunto de valores binarios a0i y un entero e0 tal que à t ! t X an X 1 a0n e e0 2 § = 2 2n 2t 2n n=1 n=1 Consideremos primero el caso de sumar 1=2t . Si ak = 1 para todo k; entonces ! à t X 1 1 1 e + t = 2e+1 2 n 2 2 2 n=1 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 a0k = 0 si k0 < k · t Consideremos ahora el caso de restar 1=2t : Si el único elemento ak distinto de 0 es a1 ; entonces e

2

µ

1 1 ¡ 2 2t

¶

= 2e¡1

Solución: No se puede escribir de forma exacta. Si suponemos à t ! à t ! X ai X ai 1 = 2e =) 1 = 9 ¢ 2e =) i 9 2 2i i=1 i=1 à t ! X t¡e 2 t¡i =) 2t¡e = 32 m 2 = 3 ai 2

P Problema 9 Dado un número ze = 2e tn=1 a2nn ; en una aritmética de precisión …nita. Calcular el número inmediatamente inferior y superior a él en dicha aritmética.

i=1

Solución: Si el número es de la forma

donde m es un número entero. Ahora bien esta igualdad es imposible porque resultaría que 3 divide a 2:

ze = 2e

1 2

entonces el inmediato superior es +

1 25

+

1 26

¡ Problema 7 Sean A = 2 12 + 213 + ¡ ¢ 23 21 + 216 + 217 : Calcular B + A y B ¡ A

+ 1 25

¢

1 27

+ B

1 28

+ =

t X 1 n 2 n=1

Si por el contrario existe un k0 > 1 tal que ak0 = 1, y tal que ak = 0 para todo k0 < k · t entonces basta tomar a0k = ak si 1 · k < k0 ; a0k0 = 0 y a0k = 1 si k0 < k · t:

1 9:

1 24

¢

pertenecen al conjunto A de números reales generados por la aritmética de precisión …nita.

3. Los números positivos más grande y más pequeño de la aritmética (teniendo en cuenta las excepciones) Solución: Ã 10 ! µ1 1 ¶ X 1 31 31 2 ¡ 211 M ayor = 2 = 2 1 2i 2 i=1 µ ¶ 1 M enor = 2¡31 210

Solución:

¡1

+

Problema 8 Sean emin , emax ; los valores mínimo y máximo del exponente e: Demostrar que si emin < e < emax ; entonces los números: Ã t ! X an 1 e § t 2 2n 2 n=1

2. El cero, el in…nito y Na. Solución: µ ¶ ¡31 1 0 = 2 2 µ ¶ 1 1 = 232 2 ¶ µ 1 32 1 N aN = 2 + 2 22

¡ ¢ 5. 2 21 ¡ 2110 : Solución: ¡ ¢ ¡ 2 12 ¡ 2110 = 20 12 + 212 + 213 +

2

1. B ¡ A = 22

1. 2; y los números más cercanos a 2 por arriba y por debajo. Solución: µ ¶ 2 1 2 = 2 2 µ ¶ 1 2 1 Si guiente = 2 + 2 210 Ã 10 ! µ1 1 ¶ X 1 2 ¡ 211 Anterior = 2 = 2 1 2i 2 i=1

4.

¡1

B + A = 23

1 29

¢

ze + 2e

y el inmediato inferior es

2e¡1

3

1 2t

t X 1 n 2 n=1

CALCULO DE CEROS DE UNA FUNCION

para cualquier otro número ze; el inmediato superior e inferior son 1 ze § 2e t 2

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]

Problema 10 Calcular las raíces del polinomio P (x) = x2 ¡ 2x + 0:01 evitando los errores de cancelación.

Solución:

Solución:

x1 x2

= =

2+

p

x =

4 ¡ 0:04 = 1:995 2

f(¡2) > N uevo Intervalo =

0:01 1:995

x = f(¡2) > N uevo Intervalo =

Problema 11 Escribir el código en fortran 77 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:

0 + (¡2) = ¡1 2 0; f(0) < 0; f (¡1) < 0 [¡2; ¡1] ¡1 + (¡2) = ¡1:5 2 0; f(¡1) < 0; f(¡1:5) > 0 [¡1:5; ¡1]

Solución:

Problema 13 Escribir el código en fortran 77 para implementar el método de la bisección

PRINT *, ’CALCULO DE LAS RAICES DE P(X)=A*X ^2+BX+C’ PRINT *, ’INTRODUZCA EL VALOR DE A’ READ *,A PRINT *, ’INTRODUZCA EL VALOR DE B’ READ *,B PRINT *,’INTRODUZCA EL VALOR DE C’ READ *,C IF (A.EQ.0 ) THEN IF (B.EQ.0 ) THEN PRINT *, ’EL POLINOMIO ES CONSTANTE’ EXIT ENDIF PRINT *, ’EL POLINOMIO ES DE GRADO 1.’ PRINT *, ’LA RAIZ ES X=’,¡C=B EXIT ENDIF D=B*B-4*A*C IF (D.LT.0 )THEN PRINT *, EL POLINOMIO NO TIENE RAICES REALES EXIT ENDIF IF (B.GT.0) THEN X1=(-B-SQRT(D))/(2*A) ELSE X1=(-B+SQRT(D)/(2*A) ENDIF X2=C/(X1*A) PRINT *,’LAS RAICES SON: ’,X1,X2 END

Solución: PRINT *,’METODO DE LA BISECCION’ PRINT *,’INTRODUCIR EL EXTREMO IZQUIERDO DEL INTERVALO’ READ *,A PRINT *,’INTRODUCIR EL EXTREMO DERECHO DEL INTERVALO’ READ *,B IF (A.GT.B) THEN PRINT *, ’INTERVALO INCORRECTO’ EXIT ENDIF PRINT *,’INTRODUCIR LA PRECISION ’ READ *,TOL IF (F(A)*F(B).GT.0) THEN PRINT *,’NO HAY CAMBIO DE SIGNO EN EL INTERVALO’ EXIT ENDIF 1 X=(A+B)/2 IF (F(X).EQ.0).OR.((B-A).LT.TOL) THEN PRINT *,’LA RAIZ DE LA FUNCION ES: ’,X EXIT ENDIF IF((F(A)*F(X)).LE.0) THEN B=X ELSE A=X ENDIF GOTO 1 END FUNCTION F(X) F=COS(X) 4

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:

END

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]

Solución: x1 = 1 ¡

Solución: 2 f (0) = 1 f (2) ¡ f(0) f (2) > 0; f(0) < 0; f (1) < 0 N uevo Intervalo = [1; 2] 1 4 x = 1¡ f (1) = f (2) ¡ f(1) 3 4 f (2) > 0; f(1) < 0; f ( ) < 0 3 4 N uevo Intervalo = [ ; 2] 3 Problema 15 Escribir el código en fortran 77 para implementar el método de la Regula-falsi

¡2 5 = 3 3

x = 0¡

Problema 17 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 ¡ ³

Solución:

¡2

¡2¡(¡3) 1¡0

´ =3

Problema 18 Escribir un programa en fortran 77 que implemente el 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 tolerancia T OL para determinar la igualdad de dos números.

PRINT *,’METODO DE LA REGULA FALSI’ PRINT *,’INTRODUCIR EL EXTREMO IZQUIERDO DEL INTERVALO’ READ *,A PRINT *,’INTRODUCIR EL EXTREMO DERECHO DEL INTERVALO’ READ *,B IF (A.GT.B) THEN PRINT *, ’INTERVALO INCORRECTO’ EXIT ENDIF PRINT *,’INTRODUCIR LA PRECISION ’ READ *,TOL IF (F(A)*F(B).GT.0) THEN PRINT *,’NO HAY CAMBIO DE SIGNO EN EL INTERVALO’ EXIT ENDIF 1 X=A-F(A)*(B-A)/(F(B)-F(A)) IF (F(X).EQ.0).OR.((B-A).LT.TOL) THEN PRINT *,’LA RAIZ DE LA FUNCION ES: ’,X EXIT ENDIF IF(F(A)*F(X)).LE.0 THEN B=X ELSE A=X ENDIF GOTO 1 END

Solución: IMPLICIT DOUBLE PRECISION (X) PRINT *,’METODO DE LA SECANTE’ PRINT *,’INTRODUCIR X0’ READ *,X0 PRINT *,’INTRODUCIR X1’ READ *,X1 IF (X0.EQ.X1) THEN PRINT *, ’LAS DOS APROXIMACIONES INICIALES COINCIDEN’ EXIT ENDIF PRINT *,’INTRODUCIR LA PRECISION ’ READ *,XTOL PRINT *,’INTRODUCIR NUMERO MAXIMO DE ITERACIONES’ READ *,NMAX DO 1 K=1,NMAX IF(ABS(X1-X0).LT.TOL) THEN PRINT *,’LA RAIZ DE LA FUNCION ES: ’,X1 EXIT ENDIF IF(XF(X1).EQ.XF(X0)) THEN PRINT *,’METODO NO CONVERGE’ EXIT ENDIF X2=X1-XF(X1)*(X1-X0)/(XF(X1)-XF(X0)) X0=X1 X1=X2

FUNCTION F(X) F=COS(X) END

5

Problema 22 Aislar en intervalos las raíces del polinomio P (x) = 20x3 ¡ 45x2 + 30x ¡ 1.

1 CONTINUE PRINT *,’NUMERO MAXIMO DE ITERACIONES EXCEDIDO’ END

Solución: Teniendo en cuenta que en este caso

FUNCTION XF(X) IMPLICIT DOUBLE PRECISION (X) F=COS(X) END

1+

45 65 maxk=0;::;n¡1 j ak j =1+ = j an j 20 20

65 todas las raíces están en el intervalo [¡ 65 20 ; 20 ]: Para aislar 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

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 función de forma exacta) y quedándonos con la raíz más cercana a x0 :

65 ) 20 1 P( ) 2 P (1) 65 P( ) 20

P (¡

Solución: ¡2 + 3(x ¡ 1) + 3(x ¡ 1)2 = 0 p ¡3 + 33 x1 = 1 + 6

= ¡1260: 4 21 4 = 4

=

= 307: 75

1 por tanto hay una única raíz en el intervalo [¡ 65 20 ; 2 ]:

Problema 23 Aislar en intervalos las raíces del polinomio P (x) = 2x3 + 3x2 ¡ 12x + 1

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

Solución:

Solución:

P 0 (x) = 6x2 + 6x ¡ 12 P (x) P (2) P (2) P 0 (x) P 0 (2)

= = = = =

Intervalo Inicial [¡7; 7] P (¡7) = ¡454 P (¡2) = 21 P (1) = ¡6 P (7) = 750:

((2x + 3)x + 4)x + 5 ((7)2 + 4)2 + 5 (18)2 + 5 = 41 (2x + 7)x + 18 (4 + 7)2 + 18 = 40

Intervalos donde están las raíces: [-7,-2] [-2,1] [1,7]

INTERPOLACION DE FUNCIONES I

Problema 21 Calcular el número máximo de raíces positivas y negativas del polinomio x5 ¡ 35x3 + 30x2 + 124x ¡ 120; y localizarlas en un intervalo.

Problema 24 Calcular el polinomio interpolador de Lagrange P3 (x) de la función f (x) = sen(x) en los puntos 0; ¼2 ; ¼ y 3¼ 2 .

Solución: Teniendo en cuenta que 1+

ra¶{ces x = 1; ¡2

maxk=0;::;n¡1 j ak j = 125 j an j

por tanto el número de raíces positivas es 1 ó 3: Para estimar el número de raíces negativas cambiamos x por ¡x y miramos los signos de los coe…cientes que en este caso son:

Solución: Puesto que sen(0) = sen(¼) = 0 sólo necesitamos los polinomios base de Lagrange centrados en ¼2 y 3¼ 2 : ¡ ¢ x (x ¡ ¼) x ¡ 3¼ 2 ¼ ¢ ¡ ¼ 3¼ ¢ P 2 (x) = ¼ ¡ ¼ ¡ 2 2 2 ¡¼ ¡2 ¢ x (x ¡ ¼) x ¡ ¼2 ¡ ¢ ¡ ¢ P 3¼ (x) = 3¼ 3¼ 3¼ ¼ 2 2 2 ¡¼ 2 ¡ 2

¡++¡¡

P (x) = P ¼2 (x) ¡ P 3¼ (x) 2

las raíces del polinomio están en el intervalo [¡125; 125]: Para calcular el número máximo de raíces positivas miramos los cambios de signo de los coe…cientes, en este caso los signos son: +¡++¡

Por tanto el polinomio interpolador es

por tanto el número de raíces negativas son 0 ó 2:

6

Problema 25 Calcular la expresión del error interpolación al aproximar la función f(x) = sen(x) en el intervalo [0; 2¼] interpolando en los puntos 0; ¼2 ; ¼, 3¼ 2 : y acotarlo superiormente.

Problema 28 Calcular el polinomio interpolador de Lagrange 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;

Solución: El error de interpolación viene dada por la expresión: µ ¶ sen(») ³ ¼´ 3¼ f (x) ¡ PN (x) = x x¡ (x ¡ ¼) x ¡ 4! 2 2

¼ 2 f[0; ] = 2 ¼ ¼ 2 f [ ; ¼] = ¡ 2 ¼ 3¼ 2 f[¼; ] = ¡ 2 ¼ ¼ 4 f[0; ; ¼] = ¡ 2 2 ¼ ¼ 3¼ f [ ; ¼; ] = 0 2 2 ¼ 3¼ 8 f [0; ; ¼; ] = 2 2 3¼ 3 por tanto, el polinomio interpolador es 2 4 ³ ¼´ 8 ³ ¼´ P (x) = x ¡ 2 x x ¡ + 3x x ¡ (x ¡ ¼) ¼ ¼ 2 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 µ ¶ 1 ³ ¼´ 3¼ jf (x) ¡ PN (x)j · 2¼ 2¼ ¡ (2¼ ¡ ¼) 2¼ ¡ 4! 2 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:

Problema 29 Calcular el polinomio interpolador de Lagrange 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: Según las fórmulas vistas en teoría el error viene dado por la expresión: ¯ ¯µ ¶ maxx2[a;b] ¯f N+1) (»)¯ b ¡ a N+1 j f(x) ¡ PN (x) j· (N + 1)!2N 2

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 1 j f(x) ¡ PN (x) j· 6!25

Solución: En el primer caso, las diferencias divididas son f [x0 ] = 1; f[x1 ] = 2; f [x2 ] = 8; f [x3 ] = 16:

µ ¶6 1 = 6: 78 £ 10¡7 2

f [x0 ; x1 ] = 1 f [x1 ; x2 ] = 3 f [x2 ; x3 ] = 8 2 f [x0 ; x1 ; x2 ] = 3 5 f [x1 ; x2 ; x3 ] = 3 1 f[x0 ; x1 ; x2 ; x3 ] = 4 y el polinomio interpolador es:

Problema 27 Interpolar la función f (x) = x210+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 -1 ! 5 1!5

3

2 1 P (x) = 1 + x + x(x ¡ 1) + x(x ¡ 1)(x ¡ 3) 3 4 Si tomamos ahora los puntos en orden inverso: f [x0 ] = 16; f [x1 ] = 8; f[x2 ] = 2; f [x3 ] = 1:

-1 0

0 -1

-3

2!2

f [x0 ; x1 ] = 8 f [x1 ; x2 ] = 3 f [x2 ; x3 ] = 1 5 f [x0 ; x1 ; x2 ] = 3 2 f [x1 ; x2 ; x3 ] = 3 1 f[x0 ; x1 ; x2 ; x3 ] = 4

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

7

donde » 2 [0; ¼4 ]: Para que la aproximación sea la mejor dentro de una aritmética de 32 bits tiene que cumplirse

El polinomio interpolador es: 5 1 P (x) = 16+8(x¡4)+ (x¡4)(x¡3)+ (x¡4)(x¡3)(x¡1) 3 4

j Pn (x) ¡ sen(x) j · 2¡24 = 5: 96 £ 10¡8 sen(x)

como puede observarse, al cambiar el orden de los puntos de interpolación, el polinomio de Lagrange expresado a través de las diferencias divididas cambia totalmente, salvo el último coe…ciente 14 que es el mismo en ámbos casos 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.

por otro lado, en el intervalo [0; ¼4 ] se veri…ca ¡ ¢ sen ¼4 x · sen(x) ¼ 4

por tanto: ¡ ¼ ¢2n+2 j Pn (x) ¡ sen(x) j · 4 sen(x) (2n + 2)!

Problema 30 Dada una función f(x); y una secuencia 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órmulas dadas en el método de Muller para el cálculo de las derivadas f 00 (xn¡1 ) y f 0 (xn¡1 ):

para n = 4 se tiene que

Solución: Si utilizamos las diferencias divididas para interpolar obtenemos f [xn¡1 ] = f (xn¡1 )

por tanto n = 4 determina la mejor aproximación en una aritmética de 32 bits.

¡ ¼ ¢2n+2 4

(2n + 2)!

f (xn¡1 ) ¡ f (xn¡2 ) xn¡1 ¡ xn¡2 f (xn¡2 ) ¡ f (xn¡3 ) xn¡2 ¡ xn¡3 f [xn¡1 ; xn¡2 ] ¡ f [xn¡2 ; xn¡3 ] xn¡1 ¡ xn¡3

f [xn¡1 ; xn¡2 ] = f [xn¡2 ; xn¡3 ] = f [xn¡1 ; xn¡2 ; xn¡3 ] =

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 de…nir el sen(x) y cos(x) para cualquier valor de x a partir de su de…nición en [0; ¼4 ]; por tanto, en este problema sólo tenemos que de…nir las funciones trigonométricas en [0; ¼4 ] a partir de su de…nición en [0; ¼8 ]: Basta tener en cuenta las relaciones: ( ¼ ¼ cos ¡ x ¢[0; 8 ] (x)2 ¡ x ¢ si x · ¼8 cos[0; ¼4 ] (x) = 2 cos[0; ¼ ] 2 ¡ sin[0; ¼8 ] 2 si x > 8

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 )

8

sen[0; ¼4 ] (x) =

que corresponde a las fórmulas utilizadas por el método de Muller.

sen(x) u Pn (x) = x ¡

sen¡ [0;¢¼8 ] (x) ¡ ¢ si x · 2 cos[0; ¼8 ] x2 sin[0; ¼8 ] x2 si x >

¼ 8 ¼ 8

Solución: En primer lugar, la función cos(x) la desarrollamos por serie de Taylor como

Solución: El desarrollo de Taylor en 0 del sen(x) viene dado por: 5

½

Problema 33 Calcular los polinomios necesarios para interpolar las funciones trigonométricas cos(x) y sen(x) en el intervalo [0; ¼8 ] en una aritmética de 32 bits

Problema 31 Aproximar la función sen(x) en el intervalo [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.

3

= 2: 46 £ 10¡8

cos(x) u Pn (x) = 1 ¡

2n+1

x x x + + :::: + (¡1)n 3! 5! (2n + 1)!

x2 x4 x2n + + :::: + (¡1)n 2! 4! (2n)!

y el error máximo cometido por el desarrollo de Taylor en un punto x 2 [0; ¼8 ] es

y el error máximo cometido por el desarrollo de Taylor en un punto x 2 [0; ¼4 ] es

2n+1

¼ (x) j Pn (x) ¡ cos(x) j· sen( ) 8 (2n + 1)!

2n+2

¼ (x) j Pn (x) ¡ sen(x) j· sen( ) 4 (2n + 2)! 8

donde » 2 [0; ¼8 ]: Para que la aproximación sea la mejor dentro de una aritmética de 32 bits tiene que cumplirse

ANALISIS NUMERICO MATRICIAL I

j Pn (x) ¡ cos(x) j · 2¡24 = 5: 96 £ 10¡8 cos(x)

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 A0 u = b0 donde A0 es una matriz triangular superior.

por tanto: ¡ ¢2n+1 j Pn (x) ¡ cos(x) j ¼ ¼8 · tan( ) cos(x) 8 (2n + 1)!

para n = 3 se tiene que ¡ ¢2n+1 ¼ ¼8 = 1: 18 £ 10¡7 tan( ) 8 (2n + 1)!

Solución: Escribimos la matriz A0 de la siguiente manera, 0 1 a11 a12 ¢ ¢ ¢ a1;n¡2 a1;n¡1 a1n B 0 a22 ¢ ¢ ¢ a2;n¡2 a2;n¡1 a2n C B C B .. C .. . . .. . . . . B . C . . . . . B C B 0 0 0 an¡2;n¡2 an¡2;n¡1 an¡2;n C B C @ 0 0 0 0 an¡1;n¡1 an¡1;n A 0 0 0 0 0 an;n

con lo cual ya estamos muy cerca de la precisión óptima. Para n = 4 ¡ ¢2n+1 ¼ ¼8 tan( ) = 2: 53 £ 10¡10 8 (2n + 1)!

por tanto n = 4 determina la mejor aproximación en una aritmética de 32 bits.

En el remonte se empiezan a calcular los ui de abajo hacia arriba. Las operaciones que se realizan vienen dadas por: n un = abnn

Análogamente, para la función sen(x) tenemos sen(x) u Pn (x) = x ¡

x3 x5 x2n+1 + + :::: + (¡1)n 3! 5! (2n + 1)!

y el error máximo cometido por el desarrollo de Taylor en un punto x 2 [0; ¼4 ] es ¼ (x)2n+2 j Pn (x) ¡ sen(x) j· sen( ) 8 (2n + 2)!

donde » 2 [0; ¼8 ]: Para que la aproximación sea la mejor dentro de una aritmética de 32 bits tiene que cumplirse

por otro lado, en el intervalo [0; ¼8 ] se veri…ca ¡ ¢ sen ¼8 x · sen(x) ¼ 8

por tanto:

¡ ¼ ¢2n+2 j Pn (x) ¡ sen(x) j · 8 sen(x) (2n + 2)!

un¡2 =

bn¡2 ¡(an¡2;n un +an¡2;n¡1 un¡1 ) an¡2;n¡2

un¡3 = .. .

bn¡3 ¡(an¡3;n un +an¡3;n¡1 un¡1 +an¡3;n¡2 un¡2 ) an¡3;n¡3

Sumas n¡1 .. .

Multiplic: n¡1 .. .

Divisiones 1 .. .

T otal 2n ¡ 1 .. .

3 2 1 0

3 2 1 0

1 1 1 1

7 5 3 1

A partir de esta tabla podemos calcular el total de operaciones sumando por columnas:

para n = 3 se tiene que ¡ ¼ ¢2n+2 8

bn¡1 ¡an¡1;n un an¡1;n¡1

En la siguiente tabla se muestra el número de operaciones que se realizan en cada iteración:

j Pn (x) ¡ sen(x) j · 2¡24 = 5: 96 £ 10¡8 sen(x)

(2n + 2)!

un¡1 =

= 1: 402 679 863 £ 10¡8

Sumas = 0 + 1 + 2 + 3 + : : : + n ¡ 1 =

por tanto n = 3 determina la mejor aproximación en una aritmética de 32 bits.

(n¡1)n 2

M ultiplicac: = 0 + 1 + 2 + 3 + : : : + n ¡ 1 =

(n¡1)n 2

Divisiones = 1 + 1 + 1 + 1 + : : : + 1 = n

Problema 34 Como se puede obtener la función yx ; donde x; y son números reales, utilizando las funciones ex y ln(x):

Total = 1 + 3 + 5 + 7 + : : : + 2n ¡ 1 = = Sumas + M ultiplicac: + Divisiones =

Solución: Se utiliza la equivalencia

=

yx = ex ln y

(n¡1)n 2

+

(n¡1)n 2

+ n = n2

El orden del algoritmo es entonces O(n2 ): 9

Problema 36 Resolver µ ¶µ ¶por el µ método ¶ de Gauss el sistema ¡1 2 x 3 = 2 ¡1 y 0

A continuación obtenemos el total de operaciones en cada iteración sumando por columnas: 1a Iteración:

Solución: µ ¶µ ¶ µ ¶ ¡1 2 x 3 = ! 2 ¡1 y 0 µ ¶µ ¶ µ ¶ 2 ¡1 x 0 = ! ¡1 2 y 3 µ ¶µ ¶ µ ¶ 2 ¡1 x 0 = ! 23 y = y ! y = 2 ! 0 23 y 3 x = 22 = 1

Divisiones = 1 + 1 + : : : + 1 = n ¡ 1 M ultiplicac: = n + n + : : : + n = n(n ¡ 1) Sumas = n + n + : : : + n = n(n ¡ 1) 2a Iteración: Divisiones = 1 + 1 + : : : + 1 = n ¡ 2 M ultiplicac: = (n ¡ 1) + (n ¡ 1) + : : : + (n ¡ 1) =

Problema 37 Calcular el número de operaciones básicas necesarias para descomponer el sistema Au = b en el sistema A0 u = b0 utilizando el método de Gauss, y teniendo en cuenta la siguiente relación M¡1 X k=1

=(n ¡ 1)(n ¡ 2) Sumas = (n¡1)+(n¡1)+: : :+(n¡1) = (n¡1)(n¡2) .. .

1 1 1 k2 = M 3 ¡ M 2 + M 3 2 6

(n-1)a Iteración: Divisiones = 1

Solución:

0

a11 B a21 B A=B . @ .. an;1

a12 a22 .. .

¢¢¢ ¢¢¢ .. .

a1n a2n .. .

an;2

¢¢¢

an;n

M ultiplicac: = 2

1

Sumas = 2

C C C A

Total operaciones1 : Divisiones = (n¡1)+(n¡2)+(n¡3)+: : :+1 =

En cada iteración se realizan las siguientes operaciones:

M ultiplicac: = n(n ¡ 1) + (n ¡ 1)(n ¡ 2) + : : : + 2 = =((n ¡ 1) + 1)(n ¡ 1) + ((n ¡ 2) + 1)(n ¡ 2) + : : : =

Para cada iteración (i): Para cada …la (j) ³ ´ ¤ aaii ii ¤aj1 ¡ ai1

³

aji aii

=(n ¡ 1)2 + (n ¡ 1) + (n ¡ 2)2 + (n ¡ 2) + : : : = 3

= 2n ´

: : : ajn ¡ ain

³

aji aii

´

1a

2a .. . (n ¡ 1)a

¡3n2 +n 6

+

(n¡1)n 2

=

n3 ¡n 3

Sumas = n(n ¡ 1) + (n ¡ 1)(n ¡ 2) + : : : + 2 =

Fila 2a 3a .. .

Division. 1 1 .. .

Multiplic. n n .. .

Sumas n n .. .

na 3a 4a .. .

1 1 1 .. .

n n¡1 n¡1 .. .

n n¡1 n¡1 .. .

na .. .

1 .. .

n¡1 .. .

n¡1 .. .

na

1

2

2

n3 ¡n 3

Total=Sumas + Multiplicac: + Divisiones =

En la primera iteración, este proceso se repite N ¡ 1 veces (para las N ¡ 1 j-…las inferiores). En la segunda, se repite N ¡ 2 veces, y así sucesivamente hasta la penúltima …la, en donde sólo se realiza una vez. Iteración

n(n¡1) 2

=

n3 ¡n 3

+

n(n¡1) 2

= 23 n3 + 21 n2 ¡ 67 n 3

El orden del algoritmo es entonces O( 2n3 ): Problema 38 Implementar en FORTRAN la funcion IDESCENSO(A; b; u; N; N max) que resuelve un sistema 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 real del sistema y N max la dimensión que se utilizo para reservar la memoria de la matriz A: La función devuelve 0 si termina correctamente y 1 en caso contrario. Nota Importante: Las líneas de código tienen que ir todas numeradas y no pueden superar las 15 líneas como máximo. 11

10

+ 22 + 32 + ::: + (n ¡ 1)2 =

2n3 ¡3n2 +n 6

Solución: 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15

Solución:Tenemos que demostrar, por una parte, que At = A (A simétrica) y, por otra, que x ¹t A¹ x > 0 (A de…nida 2 positiva ).

IDESCEN SO(A; b; u; N; N max) DIMEN SION A(N max; ¤); b(¤); u(¤) DO 13 I = 1; N IF (A(I; I):EQ:0) T HEN IDESCEN SO = 1 RET URN EN DIF u(I) = b(I) DO 11 J = 1; I ¡ 1 u(I) = u(I) ¡ A(I; J) ¤ u(J) CON T INU E u(I) = u(I)=A(I; I) CONT IN U E IDESCEN SO = 0 END

1. Simétrica: t

Como jBj 6= 0; si B¹ x = 0 =) x ¹=0 Una matriz se dice de…nida positiva si se cumple que 8¹ x 6= 0; x ¹t A¹ x > 0 =) =) x ¹t A¹ x=x ¹t BB t x ¹ = (B t x ¹) B t x ¹=

Problema 39 Resolver por el método de Gauss el siguiente sistema de ecuaciones 0 10 1 0 1 0 ¡1 2 u1 1 @ ¡1 2 ¡1 A @ u2 A = @ 0 A 2 ¡1 0 u3 1

1 2 ¡1 0 1 ¡1 2 ¡1 0 A 0 ¡1 2 1

Cálculo de los elementos de la matriz B :

A = B ¢ Bt = 0 b11 0 = @ b21 b22 b31 b32

u1 =

1+u2 2

0

=1 = =

2 2

¡1 ¡1

1 1 4 5 6 A 6 26

donde la matriz B es triangular inferior. 0 1 b11 0 0 B = @ b21 b22 0 A b31 b32 b33 0 1 b11 b21 b31 B t = @ 0 b22 b32 A 0 0 b33

4. Realizamos el remonte, y obtenemos como solución:

1¡2u3 ¡1

la siguiente matriz A por el

A = B ¢ Bt;

3. Hacemos ceros en la segunda columna ³ ´ j1 filaj ¡ f ila2 ¢ aa11 ;j > 2 : 0 1 0 1 2 ¡1 0 1 2 ¡1 0 1 ¡¡! @ 0 3 ¡2 1 A @ 0 3 ¡2 1 A ¡ ceros 0 ¡1 2 1 0 0 4 4

u2 =

yi2 > 0

Solución: La descomposición por el método de Cholesky tiene la forma siguiente:

la primera:

2. Hacemos ceros en la primera columna ³ ´ j1 filaj ¡ f ila1 ¢ aa11 ;j > 1 : 0 1 0 1 2 ¡1 0 1 2 ¡1 0 1 ¡¡! @ 0 3 ¡2 1 A @ ¡1 2 ¡1 0 A ¡ ceros 0 ¡1 2 1 0 ¡1 2 1

4 4

P

Problema 41 Descomponer método de Cholesky 0 1 A=@ 1 4

Solución: Pasos en la descomposición por Gauss:

u3 =

t

2. De…nida positiva:

= y¹t y¹ =

1. Intercambiamos la tercera …la con 0 1 0 0 ¡1 2 1 ¡¡¡! @ ¡1 2 ¡1 0 A ¡ pivoteo @ 2 ¡1 0 1

t

At = (B ¢ B t ) = (B ¢ B t ) = (B t ) B t = B ¢ B t = A

b211 @ b11 b21 = b11 b31

=1

10 0 b11 0 A@ 0 b33 0 b11 b21 b221 + b222 b21 b31 + b22 b32

b21 b22 0

1 b31 b32 A = b33

1 b11 b31 b21 b31 + b22 b32 A b231 + b232 + b233

Igualamos los elementos de la matriz anterior con los elementos de la matriz A y se obtienen los siguientes resultados:

=1

Problema 40 Demostrar que si A = B ¢B t (B triangular inferior) y jBj = 6 0; entonces A es simétrica y de…nida positiva

2 Matriz de…nida positiva: 8¹ x 6= 0 =) x ¹t A¹ x > 0: Esta es la de…nició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 ¹t A¹ x=x ¹t ¸¹ x = ¸¹ xt x ¹ > 0:

11

En la siguiente tabla se muestra de forma esquematizada, el número de operaciones en cada iteración:

b211 = 1 b11 = 1

Iteraci¶ on b11 b21 = 1 b21 =

1 b11

=1

i=1

Sumas 0 0 .. .

M ultiplic: 0 0 .. .

0

0

b11 b31 = 4 b31 =

4 b11

1 n¡1

=4

1 1 .. .

i=2

b221 + b222 = 5 p p b22 = § (5 ¡ b221 ) = (4) = 2 6¡b21 b31 b22

=

6¡4 2

b231 + b232 + b233 = 26 p p b33 = (26 ¡ b231 ¡ b232 ) = (26 ¡ 16 ¡ 12 ) = 3

1

1

n¡1

n¡2

2 2 .. .

2 2 .. .

0 1 .. .

2

2

1

2(n¡2)

n¡3

2(n¡2)

.. . i=n

La descomposición queda de la siguiente manera:

0 1 .. .

1

i=3

=1

1 1 .. .

n¡1

b21 b31 + b22 b32 =6 b32 =

Divisiones 0 1 .. .

.. . n¡1

.. . n¡1

.. . 0

El total de operaciones se obtiene sumando los totales parciales de la tabla anterior:

A = B ¢ Bt = 0 10 1 1 0 0 1 1 4 = @ 1 2 0 A@ 0 2 1 A 4 1 3 0 0 3

Sumas = M ultiplicac: = = (n ¡ 1) + 2 (n ¡ 2) + 3(n ¡ 3) + : : : + (n ¡ 1) = P n3 ¡n = n¡1 i=1 i(n ¡ i) = 6

Problema 42 Calcular el número de operaciones necesarias para resolver un sistema por el método de Cholesky.

Divisiones = (n ¡ 1) + (n ¡ 2) + (n ¡ 3) + : : : + 1 = P n(n¡1) = n¡1 i=1 i = 2

Solución: Las operaciones que se realizan en cada iteración vienen dadas por:

El resultado …nal es:

Total=Sumas + M ultiplicac: + Divisiones = Iteraci¶ on i=1

i=2 .. . i=n

= 2n

Operaciones p j = 1 : b11 = a11 21 j = 2 : b21 = ab11 .. . j = n : bn1 = abn1 p11 j = 2 : b22 = a22 ¡ b221 21 b31 j = 3 : b32 = a32 ¡b b22 .. . 21 bn1 j = n : bn2 = an2 ¡b b22 .. . q ¡ ¢ j = n : bnn = ann ¡ b2n1 + : : : + b2n;n¡1

3

¡n 6

+

n(n¡1) 2

= 13 n3 ¡ 65 n + 12 n2

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:

0

B B AA¡1 = Id = B @

1 0 ¢¢¢ 0 1 ¢¢¢ .. .. . . . . . 0 0 ¢¢¢

0 0 .. . 1

1 C C C A

Si expresamos la matriz inversa de la siguiente manera: 12

0

B B AB @

c11 c21 .. .

c12 c22 .. .

¢¢¢ ¢¢¢ .. .

c1n c2n .. .

cn1

cn2

¢¢¢

cnn

1

0

C B C B C=B A @

1 0 ¢¢¢ 0 1 ¢¢¢ .. .. . . . . . 0 0 ¢¢¢

0 0 .. . 1

1

l2 u2 = b2

C C C; A

l2 = a2 ¡ m1 u1 ; u2 = .. . mn¡2 un¡2 + ln¡1 = an¡1

se pueden calcular las columnas de esa matriz a partir de N sistemas de ecuaciones de la siguiente forma: 0 1 0 1 c11 1 B c21 C B 0 C B C B C AB . C = B . C @ .. A @ .. A

ln¡1 un¡1 = bn¡1 ln¡1 = an¡1 ¡ mn¡2 un¡2 ; un¡1 =

0

cn1

0

1

c12 B c22 B AB . @ .. cn2 .. . 0 c1n B c2n B AB . @ ..

0

C B C B C=B A @ 1

cnn

0

C B C B C=B A @

0 1 .. . 0 0 0 .. . 1

ln = an ¡ mn¡1 un¡1

C C C A

El algoritmo queda de la siguiente manera: l1 = a1 u1 = bl11 Para i = 2; : : : ; n ¡ 1 mi¡1 = ci¡1 li = ai ¡ mi¡1 ui¡1 ui = blii Fin Para mn¡1 = cn¡1 ln = an ¡ mn¡1 un¡1

1

C C C ; c.q.d. A

Problema 44 Demostrar el algoritmo de Crout para descomponer matrices tridiagonales.

0

0

0

Problema 45 Resolver utilizando el método de Crout el siguiente sistema de ecuaciones 0 10 1 0 1 2 4 0 x 6 @ ¡1 0 4 A @ y A = @ 3 A 0 ¡1 0 z ¡1

la matriz tridiagonal siguiente: 1 ¢¢¢ 0 ¢¢¢ 0 C C b3 0 C C C .. A . b n¡1

cn¡1

an

La descomposición por el método de Crout genera dos matrices de la forma: 0 10 1 l1 0 : 0 1 u1 : 0 B m1 l2 B : 0 C 0 C CB 0 1 : C= A=B @ 0 A @ : : 0 0 : : un¡1 A 0 : mn¡1 ln 0 : 0 1 0

l1 B m1 B =@ 0 0

l1 u1 m1 u1 + l2 : 0

bn¡1 ; mn¡1 = cn¡1 ln¡1

mn¡1 un¡1 + ln = an

1

Solución: Consideremos 0 a1 b1 0 B c1 a2 b2 B B A = B 0 c2 a3 B .. .. .. @ . . .

b2 ; m2 = c2 l2

0 l2 u2 : :

: : : mn¡1

1

0 C 0 C A ln¡1 un¡1 mn¡1 un¡1 + ln

Igualando ambas matrices y despejando los elementos li ; ui y mi ,

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

l1 u1 = b1 l1 = a1 ; u1 =

b1 ; m1 = c1 l1

m2 = ¡1 l3 = 0 ¡ 2 (¡1) = 2

m1 u1 + l2 = a2 13

En la siguiente tabla se muestra el número de operaciones en cada iteración:

Sustituyendo estos valores en las matrices de Crout, la descomposición queda: 0 1 0 1 2 0 0 1 2 0 A = L ¢ U = @ ¡1 2 0 A ¢ @ 0 1 2 A 0 ¡1 2 0 0 1 Para resolver el sistema, se tiene en cuenta lo siguiente: Ax = b LU x = b

(U x = y)

Divisiones 1 1 1 .. .

i=n

1

1

0

= 1 + 1 + : : : + 1 = (n ¡ 1) Total=Sumas + M ultiplicac: + Divisiones =

Calculamos el valor de y a partir del sistema anterior: 0 10 1 0 1 2 0 0 y1 6 @ ¡1 2 0 A @ y2 A = @ 3 A ; 0 ¡1 2 ¡1 y3

= 3 (n ¡ 1) El orden del algoritmo es O (3n)

aplicando un algoritmo de descenso, 0

y1 @ y2 A = @ y3

6 2 3+y1 2 ¡1+y2 2

1

0

DIFERENCIACION E INTEGRACION NUMERICA

1

3 A=@ 3 A 1

Problema 47 Dados 3 puntos distintos xl ; xi ; xr demostrar que la fórmula:

Calculamos el vector x por remonte:

(xl ) (xi ) + (xr ¡ xi ) f (xxi )¡f (xi ¡ xl ) f (xxrr)¡f ¡xi i ¡xl f (xi ) ¼ xr ¡ xl

Ux = y 0 10 1 0 1 1 2 0 x1 3 @ 0 1 2 A @ x2 A = @ 3 A 0 0 1 x3 1

1 0 1 0 1 x1 3 ¡ 2x2 1 @ x2 A = @ 3 ¡ 2x3 A = @ 1 A x3 1 1

0

aproxima la derivada de f 0 (xi ) con un orden de aproximación de 2:

0

quedándonos la solución …nal x =

¡

1 1 1

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

Problema 46 Calcular el número de operaciones necesarias para resolver un sistema tridiagonal por el método de Crout.

Operaciones l1 = a1 ; u1 = bl11 m1 = c1 ; l2 = a2 ¡ m1 u1 ; u2 = .. . ln = an ¡ mn¡1 un¡1

00

(xi ) 2! (xl

f 000 (xi ) (xl 3!

¡ xi )2 +

¡ xi )3

f (xr ) = f (xi ) + f 0 (xi )(xr ¡ xi )+ +f

00

(xi ) 2! (xr

¡ xi )2 +

f 000 (xi ) (xr 3!

¡ xi )3

(xi ) (xr ¡ xi ) f (x(xl )¡f = (xr ¡ xi )[f 0 (xi )+ l ¡xi )

Solución: Las operaciones que se realizan en cada iteración vienen dadas por: Iteraci¶ on i=1 i=2 .. . i=n

M ultiplic: 0 1 1 .. .

Sumas = M ultiplicac: = Divisiones =

Ly = b

1

Sumas 0 1 1 .. .

El total de operaciones se obtiene de la tabla anterior como:

y nos queda un sistema de la forma:

0

Iteraci¶ on i=1 i=2 i=3 .. .

+f

00

(xi ) 2! (xl

¡ xi ) +

f 000 (xi ) (xl 3!

¡ xi )2 ]

(xr ) (xi ¡ xl ) f (x(xi )¡f = (xi ¡ xl )[f 0 (xi )+ i ¡xr ) b2 l2

+f

00

(xi ) 2! (xr

¡ xi ) +

f 000 (xi ) (xr 3!

¡ xi )2 ]

Sumamos las expresiones anteriores y nos queda:

14

i ¡xr )f (xl ) + (x(xl ¡x i )(xl ¡xr )

(xi ) (xi ) + (xi ¡ xl ) f (x(xrr)¡f = (xr ¡ xi ) f (x(xl )¡f ¡xi ) l ¡xi )

= (xr ¡ xi )f 0 (xi ) + (xi ¡ xl )f 0 (xi )+ +f

00

+f

00

+f

000

+f

000

(xi ) 2! (xi

extraemos el factor (xr ¡ xl ),

¡ xl )(xr ¡ xi )+

(xi ) 2! (xr

(xr ¡ xl ) f 0 (xi ) =

¡ xi )(xl ¡ xi )+

(xi ) (xr 3! (xi ) (xi 3!

¡xl )f (xr ) + (xi(x ¡ r ¡xi )

¡ xi )(xl ¡ xi )2 +

(xr ¡ xl ) f 0 (xi ) =

¡ xl )(xr ¡ xi )2

+

f

000

(xi ) (xr 3!

(xr ¡ xl ) f 0 (xi ) =

¢ 0+

+

¡ xi )(xl ¡ xi ) ((xl ¡ xi ) ¡ (xr ¡ xi ))

(xi ¡xl )

xr ¡xl

xi f (xi ) + (x ¡ i ¡xl )

+ O(h2 )

(xi ) +

(x¡xi )(x¡xl ) (xr ¡xi )(xr ¡xl ) f

(xi ) +

f 0 (xi ) = (xr ) +

f (xi ) =

(x¡xl )+(x¡xi ) (xr ¡xi )(xr ¡xl ) f

f 0 (xi ) =

f (xi ) (xi ¡xr )

+

f(xi ) (xi ¡xl )

+

xl f (xi ) (xi ¡xl )

(xi ¡xr )+(xi ¡xi ) (xl ¡xi )(xl ¡xr ) f

+

xi f (xi ) (xr ¡xi )

+

(xi ¡xl )(f (xr )¡f (xi )) + (xr ¡xi ) xi f (xi ) (xr ¡xi )

¡

xr f (xi ) (xr ¡xi ) +

xl f (xi ) (xi ¡xl ) (xi ¡xl )(f (xr )¡f (xi )) + (xr ¡xi )

(xi ¡xl )(f (xr )¡f (xi )) (xr ¡xi )(f (xi )¡f (xl )) + (xr ¡xi ) (xi ¡xl ) (xr ¡xl )

Solución: 3 2 a ! f (xi + h) = f (xi ) + hf 0 (xi ) + h2 f 00 (xi ) + h6 f 000 (xi ) + O(h4 )

(xr ) +

1. b ! f(xi ¡ h) = f (xi ) ¡ hf 0 (xi ) + h3 000 4 6 f (xi ) + O(h )

h2 00 2 f (xi )

¡

c ! f (xi ¡ 2h) = f (xi ) ¡ 2hf 0 (xi ) + 2h2 f 00 (xi ) ¡ 4h3 000 4 3 f (xi ) + O(h ) 0 1 a ¡ b ¡ 2c = 0 Sistema: @ a2 + 2b + 2c = 0 A Solución: a = a b 4c ¡ ¡ = 1 6 6 3 1; b = 3; c = ¡1:

(xi ) +

i ¡xl )+(xi ¡xi ) + (x (xr ¡xi )(xr ¡xl ) f (xr ) +

¡

´

´

Problema 49 Calcular una aproximación de la derivada tercera f 000 (xi ) de una función f (x) en un punto xi ; utilizando f(xi ); f(xi + h); f(xi ¡ h); f (xi ¡ 2h)

Evaluamos la derivada en el punto xi y desarrollamos hasta obtener el resultado: (xi ¡xl )+(xi ¡xr ) (xi ¡xr )(xi ¡xl ) f

(xr ¡xi )f (xl ) (xi ¡xl )

(xi ¡xl )f (xi ) (xr ¡xi )

simpli…cando,

r )+(x¡xi ) + (x¡x (xl ¡xi )(xl ¡xr ) f (xl )

0

¡

¡

i )¡xl f (xi )(xr ¡xi ) + xi f (xi )(x(xrr¡x ¡xi )(xi ¡xl )

Derivamos la expresión anterior y obtenemos: (x¡xl )+(x¡xr ) (xi ¡xr )(xi ¡xl ) f

(xi ¡xl )f (xr ) (xr ¡xi )

l )¡xr f (xi )(xi ¡xl ) + xi f (xi )(x(xi r¡x + ¡xi )(xi ¡xl )

i )(x¡xr ) + (x(x¡x f (xl ) l ¡xi )(xl ¡xr )

f 0 (x) =

xr f (xi )¡xl f (xi ) + (xi ¡xl )

)(f (xi )¡f (xl )) + + (xr ¡xi(x i ¡xl )

Solución: El polinomio de Lagrange es: (x¡xr )(x¡xl ) (xi ¡xr )(xi ¡xl ) f

+

(xr ¡xi )f(xl ) (xi ¡xl )

xi f (xi ) (xi ¡xl )

(xr ¡ xl ) f 0 (xi ) =

Problema 48 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 comprobar que da la misma fórmula que la presentada en el problema anterior.

f (x) =

(xi ¡xr )f(xl ) (xl ¡xi )

)(f (xi )¡f (xl )) + (xr ¡xi(x + i ¡xl )

+ O(h2 )

f (xi )¡f (xl ) f (xr )¡f (xi ) +(xr ¡xi ) xr ¡xi xi ¡xl

³

(xr ¡xi )f (xi ) (xi ¡xl )

(xr ¡ xl ) f 0 (xi ) =

(xi ) + (xr ¡ xl )f 0 (xi ) = (xr ¡ xi ) f (x(xl )¡f l ¡xi )

f 0 (xi ) ¼

³

xr f(xi ) + ¡ (x r ¡xi )

El término de la tercera derivada nos da el orden de la fórmula:

(xi ) +(xi ¡ xl ) f(x(xrr)¡f ¡xi )

(xr ¡xl )f (xi ) + (xi ¡xl )

agrupamos términos,

(xi ) (xi ) (xr ¡ xi ) f (x(xl )¡f + (xi ¡ xl ) f (x(xrr)¡f = ¡xi ) l ¡xi )

= (xr ¡ xl )f 0 (xi ) +

+

¡xr f (xi )+xl f (xi ) (xr ¡xi )

¡xl )f (xr ) + (xi(x ¡ r ¡xi )

Agrupamos por las derivadas de la función:

f 00 (xi ) 2!

(xr ¡xl )f (xi ) (xi ¡xr )

(xl )

(xi ¡xl )f (xr ) (xr ¡xi )(xr ¡xl ) +

15

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 ) + O(h) h3

f 000 (xi ) =

=

=

f (xr )¡f (xi ) f (xi )¡f (xl ) ¡ h h

h

f (xr )¡f (xi )¡f (xi )¡f (xl ) h2

= O(h2 ) =

Problema 50 Dados 3 puntos. Demostrar que la fórmula 00

f (xi ) ¼ 2

f (xr )¡f (xi ) xr ¡xi

¡

=

f (xi )¡f (xl ) xi ¡xl

f (xr ) ¼ f (xi ) + f 0 (xi ) (xr ¡ xi ) + (xr ¡ xi )2 +

(xi ) 2!

f 000 (xi ) 3!

(xr ¡ xi )3

00

(xl ¡ xi )2 +

(xi ) 2!

f 000 (xi ) 3!

(xl ¡ xi )3

Extraemos en ambas ecuaciones: f (xr )¡f (xi ) (xr ¡xi )

+f

000

(xi ) 3!

f (xl )¡f (xi ) (xl ¡xi )

+f

f 00 (xi ) 2!

¼ f 0 (xi ) +

000

(xi ) 3!

(xr ¡ xi ) +

(xr ¡ xi )2

(xl ¡ xi )2 f (xi )¡f (xl ) (xi ¡xl )

¡

+

¼

f 00 (xi ) ¼ 2

f (xi ) 2!

P 0 (x) ' +

¡2

+

f (x )¡f (x ) f (xr )¡f (xi ) ¡ xi ¡x l xr ¡xi i l

xr ¡xl

f (x )¡f (x ) f (xr )¡f (xi ) ¡ xi ¡x l xr ¡xi i l

xr ¡xl

(x ¡ xl ) +

(x ¡ xi )

Calculamos la segunda derivada, obteniendo:

¡

P 00 (x) ' 2

xr ¡xl

f (x )¡f (x ) f (xr )¡f (xi ) ¡ xi ¡x l xr ¡xi i l

xr ¡xl

; c.q.d.

Problema 53 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

f (xi )¡f (xl )

¡ (xr ¡xi ) (xi ¡xl ) xr ¡xl

+ O(h)

Problema 51 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: f 0 (xi ) ¼

Solución: Sustituyendo xl = xi ¡ h; y xr = xi + h; tenemos: f (xr )¡f (xi )

f 00 (xi ) ¼ 2

(x ¡ xl ) (x ¡ xi )

((xr ¡xi )2 ¡(xl ¡xi )2 )

f (xr )¡f (xi )

f 00 (xi ) ¼ 2

f (xi )¡f (xl ) xi ¡xl

f (xi )¡f (xl )

(xr ¡xi ) (xi ¡xl ) xr ¡xl

f 000 (xr ) 3!

xr ¡xl

(x ¡ xl ) +

Derivamos el polinomio:

Despejamos la segunda derivada y obtenemos: ¡

f (x )¡f (x ) f (xr )¡f (xi ) ¡ xi ¡x l xr ¡xi i l

00

(xr ¡ xl ) + ³ ´ 000 + f 3!(xi ) (xr ¡ xi )2 ¡ (xl ¡ xi )2 f (xr )¡f (xi )

f (xi )¡f (xl ) xi ¡xl

P (x) ' f(xl ) +

(xl ¡ xi ) +

Restamos las expresiones anteriores: f (xr )¡f (xi ) (xr ¡xi )

+ O(h2 )

Polinomio de Lagrange:

f 00 (xi ) 2!

¼ f 0 (xi ) +

(h ¡ h) +

Solución: Por las diferencias divididas de Newton obtenemos lo siguiente: À xl ! f (xl ) f (xi )¡f (xl ) + f (xr )¡f (xi ) (xl ) ¡ f (xxi )¡f xi ! f (xi ) À xi ¡xl xr ¡xi i ¡xl xr ¡ xl f (xr )¡f (xi ) xr ¡xi xr ! f (xr )

f (xl ) ¼ f (xi ) + f 0 (xi ) (xl ¡ xi ) + +f

f 000 (xi ) 3!

Problema 52 Dados 3 puntos xl < xi < xr calcular el polinomio de Lagrange que interpola a f(x) en esos 3 puntos, 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: Desarrollo de Taylor de la función en el punto xi y evaluación en xr y xl :

00

f (xr )¡2f (xi )¡f (xl ) h2

¡

La aproximación de la segunda derivada queda de la forma, f (xr ) ¡ 2f (xi ) ¡ f (xl ) f 00 (xi ) ¼ h2

xr ¡ xl

aproxima la derivada segunda de f (x) en xi con un orden de aproximación de 1:

+f

=

f (xi )¡f (xl )

¡ x ¡x +h (xi +h¡xi ) ( i i ) xi +h¡xi +h

= 16

(xi ¡xl )

f (x )¡f (x ) f (xr )¡f (xi ) +(xr ¡xi ) xi ¡x l xr ¡xi i l

xr ¡xl

=

(0) (1) (0¡1) f (4)¡f +(4¡0) f (0)¡f 4¡0 0¡1 4¡1

=

1¡0 ¡ 9¡1 4 +4 ¡1 3

=

¡6 3

= ¡2

=

¡2¡4 3

=

=

f (xi )¡f (xl ) ¡ (xr ¡xi ) (xi ¡xl ) xr ¡xl

f (xr )¡f (xi )

f 00 (xi ) ¼ 2 =2

1¡0 9¡1 (4¡0) ¡ (0¡1)

4¡1

Fxx + Fyy =

=

F (x+h;y)+F (x¡h;y)+F (x;y+h)+F (x;y¡h)¡4F (x;y) , h2

= 2 2+1 3 =2

discretizando

Problema 54 Demostrar, utilizando el desarrollo de Taylor, que las siguientes expresiones son discretizaciones del laplaciano: ¢F =

¢F =

Fi+1;j + Fi¡1;j + Fi;j+1 + Fi;j¡1 ¡ 4Fi;j h2

Problema 55 Calcular una aproximación del laplaciano 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) = 41 ; F (¡ 12 ; 0) = 14 ; F (0; 21 ) = 41 ; F (0; ¡ 12 ) = 14 ; F ( 12 ; 12 ) = 12 ; F (¡ 12 ; ¡ 21 ) = 12 ; F (¡ 21 ; 12 ) = 21 ; F ( 12 ; ¡ 12 ) =

Fi+1;j+1 + Fi¡1;j+1 + Fi¡1;j¡1 + Fi+1;j¡1 ¡ 4Fi;j 2h2 Fi+1;j + Fi¡1;j + Fi;j+1 + Fi;j¡1 ¡ 4Fi;j ¢F = h2

1 2

Solución: A partir del desarrollo de Taylor de la función F , se obtiene lo siguiente:

Solución: Si representamos estos valores en una tabla, obtenemos lo siguiente:

F (x + h; y + h) = F (x; y) + hFx + hFy + + 12 h2 (Fxx + 2Fxy + Fyy )

1 2 1 4 1 2

F (x ¡ h; y ¡ h) = F (x; y) ¡ hFx ¡ hFy + + 21 h2 (Fxx + 2Fxy + Fyy )

+ 21 h2 (Fxx ¡ 2Fxy + Fyy ) F (x + h; y ¡ h) = F (x; y) + hFx ¡ hFy + + 21 h2 (Fxx ¡ 2Fxy + Fyy )

+Fi;j¡1 ¡4Fi;j + (1 ¡ °) Fi+1;j +Fi¡1;j +Fhi;j+1 ; 2

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 ) 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 (x; y + h) = F (x; y) + hFy +

h2 2 Fyy 2

h 2

8 3

+

4 3

1 1 1 1 2+2+2+2 2 14

=4

¡ ¢ + 1 ¡ 23

p ¢ ¡ ¡ ¡p 2¡ 2¢ 1 Fx = ¡2¡ 2 p ¡ 1¢ 4h ¡ 2¡ 2

Para demostrar la segunda igualdad, tomamos las siguientes ecuaciones: 2 F (x + h; y) = F (x; y) + hFx + h2 Fxx h2 2 Fxx

=

2 3

1 1 1 1 4+4+4+4 1 4

=

Problema 56 Demostrar que las máscaras

Fi+1;j+1 + Fi¡1;j+1 + Fi¡1;j¡1 + Fi+1;j¡1 ¡ 4Fi;j 2h2

F (x ¡ h; y) = F (x; y) ¡ hFx +

2 3

¢F (0; 0) =

Fxx + Fyy =

F (x; y ¡ h) = F (x; y) ¡ hFy +

1 4

i¡1;j¡1 +Fi+1;j¡1 ¡4Fi;j ¢F = ° Fi+1;j+1 +Fi¡1;j+1 +F2h + 2

Sumamos estas cuatro ecuaciones,

¢F =

0

1 2 1 4 1 2

El valor de h es 12 : Aproximamos el laplaciano promediando las dos expresiones del ejercicio anterior. Si no realizáramos este promediado, no se tendrían en cuenta todos los valores de la función.

F (x ¡ h; y + h) = F (x; y) ¡ hFx + hFy +

=

1 4

p ¢ ¡ ¡ 2¡ 2 1 0p ¢ Fy = 4h ¡ 2¡ 2

0 0 0

p ¢ ¡ 2 ¡ ¡p 2 ¢ 2¡ 2 p ¡ 1¢ 2¡ 2

¡p ¢ 2¡1 ¡p 0 ¢ 2 2¡1

¡2

p ¢ ¡ ¡ 2¡ 2 ¡ 0p ¢ 2¡ 2

dan lugar a una discretización del gradiente tal que su norma euclídea es invariante por rotaciones de 45 grados.

Fyy

Sumamos estas expresiones y obtenemos:

Solución: Procedemos de la misma forma que al calcular el valor de ° en el caso del laplaciano.

F (x + h; y) + F (x ¡ h; y) + F (x; y + h) + +F (x; y ¡ h) = 4F (x; y) + h2 Fxx + h2 Fyy 17

Consideramos una función que tiene los siguientes valores en un entorno de un punto (hi0 ; hj0 ) : 1 0 0

1 0 0

cuyas máscaras son las que se muestran en el enunciado del problema.

1 0 0

Problema 57 Calcular una aproximación del gradiente 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) = 21 ; F (¡ 12 ; 0) = ¡ 21 ; F (0; 12 ) = ¡ 21 ; F (0; ¡ 12 ) = 1 1 1 1 1 1 1 2 ; F ( 2 ; 2 ) = 0; F (¡ 2 ; ¡ 2 ) = 0; F (¡ 2 ; 2 ) = ¡1; 1 1 F(2; ¡2) = 1

Calculamos el valor del gradiente en el punto central de la siguiente manera: 0 0 Fx = ° 2h + (1 ¡ °) 4h =0 ¡2 1° 1 1¡° 1 Fy = (1 ¡ °) ¡1 2h + ° 4h = ¡ 2 h ¡ 2 h = ¡ 2h ¡ 1 ¢ r1 F (hi0 ; hj0 ) = (Fx ; Fy ) = ¡ 2h ;0

Solución: Los valores de la función en una tabla quedan de la siguiente manera: 0

Rotamos la función anterior 45o : 1 1 0

1 0 0

¡1 2

¡1

0 0 0

1° 4h

=

1 °¡2 4 h

¡1 1 1¡° Fy = (1 ¡ °) ¡1 2h + ° 4h = ¡ 2 h ¡

1° 4h

=

1 °¡2 4 h

r2 F (hi0 ; hj0 ) = (Fx ; Fy ) =

1 °¡2 4 h

=2

(1; 1)

=

kr1 F (hi0 ; hj0 )k = kr2 F (hi0 ; hj0 )k q ¡ ¢2 1 2 14 °¡2 2h = h p 1 1 2h = 4h 2 j° ¡ 2j ( p p2 = ¡ (° ¡ 2) ! ° = 2 ¡ 2 2 p p2 = (° ¡ 2) ! ° = 2 + 2 2

¡p

¢ 2¡1

¡p

¢ 2¡1

p 2¡1 h

= ¡ 12

0

1 2

¡

1 1 2+2

4h + 2 ¡ p 1 2¡ 2 1 = 2h 2 h

p ¢ 1+1 2 4h =

p 2¡1 h

¡

¡1 1 2 ¡2

4h

p ¢ ¡ + 2 ¡ 2 ¡1¡1 4h =

p 1 2¡ 2 2 h

1 = ¡ 2h

y obtenemos el valor del gradiente: µ ¶ µ 1 ¶ µ ¶ Fx 1 2h rF = = = 1 Fy ¡1 ¡ 2h Este vector nos da la dirección de máximo ascenso, que en este caso será en diagonal hacia arriba a la derecha.

Sustituyendo este valor en las expresiones de Fx ; Fy tenemos:

Problema 58 Calcular analíticamente y numéricamente la matriz gradiente en el punto (1; 1) (utilizar h = 0:1) de la función: ½ 2 x + y2 ¡ 1 f (x; y) = x¡y

¡Fi¡1;j Fx = (1 ¡ °) Fi+1;j2h + Fi+1;j+1 ¡Fi¡1;j+1 +Fi+1;j¡1 ¡Fi¡1;j¡1 4h

= ¡p ¢ Fi+1;j ¡Fi¡1;j =2 2¡1 + 4h p ¢ ¡ +Fi+1;j¡1 ¡Fi¡1;j¡1 + 2 ¡ 2 Fi+1;j+1 ¡Fi¡1;j+14h

Solución: Analíticamente µ ¶ µ ¶ 2x 2y 2 2 rf(x; y) = ! rf (1; 1) = 1 ¡1 1 ¡1

¡Fi;j¡1 Fy = (1 ¡ °) Fi;j+12h +

Fi+1;j+1 ¡Fi+1;j¡1 +Fi¡1;j+1 ¡Fi¡1;j¡1 4h

1 2

=2

p La solución válida es ° = 2 ¡ 2, ya que el gradiente r2 F debe ser negativo en sus dos derivadas.

+°

¡1 2

+ ¡p ¢ ¡Fi;j¡1 Fy = 2 2 ¡ 1 Fi;j+14h + p ¢ ¡ +Fi¡1;j+1 ¡Fi¡1;j¡1 + 2 ¡ 2 Fi+1;j+1 ¡Fi+1;j¡14h =

Calculamos las normas de los gradientes e igualamos:

+°

1

0

Sustituimos estos valores en las derivadas de la función: ¡p ¢ ¡Fi¡1;j Fx = 2 2 ¡ 1 Fi+1;j4h + p ¢ ¡ +Fi+1;j¡1 ¡Fi¡1;j¡1 + 2 ¡ 2 Fi+1;j+1 ¡Fi¡1;j+14h =

y calculamos su gradiente: ¡1 1 1¡° Fx = (1 ¡ °) ¡1 2h + ° 4h = ¡ 2 h ¡

1 2

=

¡p ¢F ¡F = 2 2 ¡ 1 i;j+14h i;j¡1 + p ¢F ¡ ¡F +F ¡F + 2 ¡ 2 i+1;j+1 i+1;j¡14h i¡1;j+1 i¡1;j¡1 ;

Numéricamente, si llamamos f1 (x; y) = x2 + y 2 ¡ 1 y f2 (x; y) = x ¡ y; entonces aplicando las máscaras vistas en teoría tenemos 18

2. n = 3

f1 (x; y) = x ¡ y @f1 (1;1) @x ¡ = p ¢ 1 4(0:1) 2 ¡ p2 (f1 (1 + 0:1; 1 ¡ 0:1) ¡ f1 (1 ¡ 0:1; 1 ¡ ¢ 1 4(0:1) 2 ¡ 2 (f1 (1 + 0:1; 1 + 0:1) ¡ f1 (1 ¡ 0:1; 1 ¡p ¢ 1 2 ¡ 1 (f1 (1 + 0:1; 1) ¡ f1 (1 ¡ 0:1; 1)) = 4(0:1) 2

3 P

¡ 0:1)) + + 0:1)) +

= 0:55555555555 ¢ P (0:7745966692) +

0:585 79 + 0:585 79 + 0:828 43 = 2:0

@f1 (1;1) @y ¡ = p ¢ 1 4(0:1) 2 ¡ p2 (f1 (1 ¡ 0:1; 1 + 0:1) ¡ f1 (1 ¡ 0:1; 1 ¡ ¢ 1 4(0:1) 2 ¡ 2 (f1 (1 + 0:1; 1 + 0:1) ¡ f1 (1 + 0:1; 1 ¢ ¡p 1 2 ¡ 1 (f1 (1; 1 + 0:1) ¡ f1 (1; 1 ¡ 0:1)) = 4(0:1) 2

+0:88888888 ¢ P (0) + +0:55555555555 ¢ P (¡0:7745966692) =

¡ 0:1)) +

= ¡: 4

¡ 0:1)) +

¢ R1 ¡ El valor exacto de la integral es ¡1 x3 ¡ x4 dx = ¡ 52 = ¡: 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.

0:585 79 + 0:585 79 + 0:828 43 = 2:0

De la misma forma, para f2 (x; y) tenemos @f2 (1;1) @x ¡ = p ¢ 1 4(0:1) 2 ¡ 2 (f2 (1 + 0:1; 1 ¡ 0:1) ¡ f2 (1 ¡ 0:1; 1 p ¢ ¡ 1 4(0:1) 2 ¡ 2 (f2 (1 + 0:1; 1 + 0:1) ¡ f2 (1 ¡ 0:1; 1 ¡p ¢ 1 2 ¡ 1 (f2 (1 + 0:1; 1) ¡ f2 (1 ¡ 0:1; 1)) = 4(0:1) 2

¡ 0:1)) + + 0:1)) +

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 cos(x)dx

0:292 89 + 0:292 89 + 0:414 21 = 1

@f2 (1;1) @y ¡ = p ¢ 1 4(0:1) 2 ¡ 2 (f2 (1 ¡ 0:1; 1 + 0:1) ¡ f2 (1 ¡ 0:1; 1 p ¢ ¡ 1 4(0:1) 2 ¡ 2 (f2 (1 + 0:1; 1 + 0:1) ¡ f2 (1 + 0:1; 1 ¡p ¢ 1 2 ¡ 1 (f2 (1; 1 + 0:1) ¡ f2 (1; 1 ¡ 0:1)) = 4(0:1) 2

¡ 0:1)) + ¡ 0:1)) +

¡¼ 2

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. ! µ ¶ 2 2 rf (1; 1) = 1 ¡1

Solución: ¼ R¼ 1. ¡2¼ cos(x)dx = sin(x)]¡2 ¼ = 1 ¡ (¡1) = 2 2

2

R

¼ 2

cos(x)dx =

R1

cos( ¼2 t) ¼2 dt =

¡1 ¡¼ 2 ¡¼¢ 2¼ 2¼ cos (0) + cos 32 32 4 p = 13 2¼ + 31 ¼ = 2: 528 2

Problema 59 Aproximar el valor de la siguiente integral, utilizando las fórmulas de Legendre para n = 2 y n = 3 Z 1 ¡ 3 ¢ x ¡ x4 dx

2¼ 32

¡ ¢ cos ¡ ¼4 +

Problema 61 Demostrar, utilizando los ceros y pesos asociados a los polinomios de Legendre, cual sería la fórmula de integración numérica de Legendre utilizando un sólo punto de interpolación. Cual sería su exactitud?

¡1

Cual es el valor exacto de la integral? Solución: N ¢ R1 ¡ 3 P 4 x ¡ x dx ' wk P (xk ) ¡1

Solución: El polinomio de Legendre para un solo punto:

k=0

L1 (x) = x ! x0 = 0

¡ ¢ P (x) = x3 ¡ x4

Calculamos el peso asociado: Z 1 1 w0 = dx = 2 ¡1 1

1. n = 2 2 P

wk P (xk ) =

k=1

Por lo tanto, la fórmula de integración numérica de Legendre es: Z

wk P (xk ) =

k=1

1

= 1 ¢ P (0:5773502692 + 1 ¢ P (¡0:5773502692) =

¡1

f (x) dx ' 2 ¢ f (0) ;

y su exactitud sería igual a 1(2N ¡ 1 = 1).

= ¡: 222 22 19

Problema 62 A partir de los ceros y de los pesos asociados 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]

Problema 64 Calcular de forma exacta la integral Z 1 ¢ ¡ 3 2 x ¡ x2 e¡x dx ¡1

utilizando los polinomios de Hermite.

Solución: Para encontrar los puntos x ^k ; y los pesos w ^k , hay que hacer un cambio de variable en la integral: Z b N X f (x) dx ' w ^k f (^ xk ) a

Solución: De forma analítica la integral da como resultado: ¢ R1 ¡ 3 p 2 x ¡ x2 e¡x dx = ¡ 21 ¼ = ¡: 886 23 ¡1

k=1

Hacemos el siguiente cambio de variable:

Utilizando el método de integración numérica:

x = (b¡a)t+(b+a) 2 dx = b¡a 2 dt

¡ ¢ f (x) = x3 ¡ x2

este cambio representa la recta que pasa por los puntos -1, 1 para t = a; b; respectivamente.

2 ¢ R1 ¡ 3 P 2 wk f (xk ) x ¡ x2 e¡x dx = ¡1 k=1

Z

b

Z

1

µ

¶

= w1 f (x1 ) + w2 f (x2) =

(b ¡ a) t + b + a b ¡ a dt 2 2 a ¡1 µ ¶ Z b N X b¡a (b ¡ a) x ~k + b + a f (x) dx ' w ~k f 2 2 a f (x) dx =

f

= 0:8862269255 ¢ f (¡0:707106781) + +0:8862269255 ¢ f (0:707106781) =

k=1

= ¡: 886 23

de donde se deduce que los cambios a realizar son de la forma x ^k = (b¡a)~xk2+(b+a) ; w ^k = (b¡a) ~k 2 w

Problema 65 Aproximar, utilizando dos puntos de aproximación, el valor de la integral: Z 1 1 dx 2 ¡1 1 + x

Problema 63 Utilizar el resultado del problema anterior para calcular de forma exacta la siguiente integral Z 1 ¡ 2 ¢ x ¡ x3 dx

Solución:

0

Solución: El resultado de la integral calculada de forma analítica, da el siguiente resultado: R1¡ 2 ¢ 1 3 ¡2 0 x ¡ x dx = 12 = 8: 333 3 £ 10

R1

1 dx ¡1 1+x2 2

f (x) =

Aplicando el método de integración numérica: ¡ ¢ f (x) = x2 ¡ x3

R1

ex 1+x2

1 ¡1 1+x2 dx

3 ¢ R1¡ 2 P 3 x ¡ x dx = wk f (xk ) = 0

=

R1

2

2 ex e¡x dx ¡1 1+x2

=¼

' w1 f(x1 ) + w2 f(x2 ) =

= 0:8862269255 ¢ f (¡0:707106781) + +0:8862269255 ¢ f (0:707106781) =

k=1

= w1 f (x1 ) + w2 f (x2 ) + w3 f (x3 ) = ¡ 1¡0 ¢ ¡ ¢ ¡ ¢ = (w ~0 f x~02+1 + w ~1 f x~12+1 + + 2 ¡ x~2 +1 ¢ w ~2 f )= 2 ¡ ¢ = 21 (0:55555556 ¢ f 0:7745966692+1 + 2 ¡ ¢ +0:8888888889 ¢ f 0+1 + 2 ¡ ¡0:7745966692+1 ¢ +0:55555556 ¢ f )= 2

= 1: 948 2 Problema 66 Calcular de forma exacta la integral Z 1 ¡ 3 ¢ x ¡ x2 e¡x dx 0

utilizando los polinomios de Laguerre.

= 8: 333 3 £ 10¡2

20

¢ R1 ¡ 3 4 ¡1 x ¡ x dx =

Solución: ¢ R1¡ 3 x ¡ x2 e¡x dx = 4 0

=

1 ¢ ¡x R1¡ 3 P 2 wk f (xk ) = e dx = x ¡ x 0

¢ ¢ R1¡ R0 ¡ 3 x ¡ x4 dx + 0 x3 ¡ x4 dx ' ¡1

'

k=0

= 0:8535533903 ¢ f (0:585786438) +

'

= 4:0

+ Problema 67 Calcular una fórmula de aproximación numérica de la integral siguiente Z 1 f (x)e¡x dx

Solución: Para calcular esta integral realizamos un cambio de variable ½ ¾ R1 t=x¡a ¡t = 0 f(t)e dx = dt = dx R1 R1 = a f (x ¡ a)e¡x+a dx = ea a f (x ¡ a)e¡x

ea

N P

R1 R1

¡1

=

a

k=1

6

k+1 +xk 2

f (0)+f (¡1)+4f ( ¡1+0 ) 2 6 f (1)+f (0)+4f ( 1+0 2 ) 6

F (x; y) dydx = N P

w ~k

k=1

= N P

f (xk+1 )+f (xk )+4f

¡1 ¡1

w ~k f (~ xk )

f (x ¡ a)e¡x = ea

+

#1

(xk+1 + xk )

=

(0 + 1) +

(1 ¡ 0) =

Solución:

k=0

R1

(xk+1 + xk )

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].

donde a es un número real cualquiera

f(t)e¡t dx =

´

#0

5 = ¡ 12 = ¡: 416 67

a

0

³x

k+1 +xk 2

0

+0:1464466093 ¢ f (3:414213562) =

R1

´

6

= w0 f (x0 ) + w1 f (x1 ) =

+

³x

f (xk+1 )+f(xk )+4f

N P

w ~k

k=1

wk f (xk ¡ a)

=

R1

¡1 F

Ã

N P

N R1 P ¡1

w ~k F (~ xk ; y) dy =

k=1

(~ xk ; y) dy

w ~j F (~ xk ; y~k )

j=1

!

N P ~ k F (~ W xk ; y~k ) ;

k;j=1

Para que estas dos igualdades sean equivalentes, basta hacer: xk = x ~k + a wk = e¡a w ~k

donde

~k = w W ~k w ~j w ~k = w ~j =

Problema 68 Aproximar, por el método de Simpson, la integral Z 1 ¡ 3 ¢ x ¡ x4 dx

Q (x¡~ xi ) Q i6=k ¡1 xk ¡~ xi ) i6=k (~

R ¡1

R ¡1 ¡1

Q (y¡~ yi ) Q i6=k yk ¡~ yi ) i6=k (~

y los x ~k e y~k son los ceros del polinomio de Legendre.

¡1

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].

utilizando únicamente el valor de la función en los puntos: ¡1; ¡ 12 ; 0; 12 y 1: Solución: ¢ R1 ¡ 3 4 dx = ¡ 25 = ¡: 4 ¡1 x ¡ x

Solución: RbRd

Aplicamos el método de Simpson: ¡ ¢ f (x) = x3 ¡ x4

a

21

c

F (x; y) dydx =

N Rd P c

k=1

w ~k F (~ xk ; y) dy =

N P

=

w ~k

k=1 N P

=

w ~k

=

c

Ã

k;j=1

Solución: Si calculamos el resultado de la integral de forma analíta, nos queda,

F (~ xk ; y) dy

N P

j=1

k=1 N P

Rd

w ~j0 F

!

R1 R2

x ¡1 0 1+ey2

(~ xk ; y~k )

=

w ~k w ~j0 F (~ xk ; y~k ) ;

R2 0

ahora bien, teniendo en cuenta los resultados obtenidos al integrar en una variable tenemos que :

R2 0

¡1 ¡1

=

=

1 ¡1 1+ey2

1. P (y) =

2 2 ¡1 x dx ¡1 y dy

k=1

=w ~k ;

dy =

por Hermite,

R1

w ~k P (~ xk )

=x ~k + 1;

= w1 P (x1 ) + w2 P (x2 ) =

R1

x2 y 2 dxdy =

2 P

(b¡a) ~k 2 w

(b+a) 2

= 2:0

2

R1

+

= : 444 44

P (x) = x2

R1 R1

(b¡a)~ xk 2

= (0:5773502692 + 1) + (¡0:5773502692 + 1) =

Utilizando la fórmula de integración numérica:

P (y) = y

wk P (xk ),

xdx =

Solución: El resultado de la integral es: y dxdy =

dy = 2: 144 3

P (x) = x

utilizando integración numérica.

¡1 ¡1 x

dy =

tenemos:

¡1

4 9

1 P

2 ¡1 1+ey2

k=0

1. wk =

Problema 71 Calcular de forma exacta la integral Z 1Z 1 x2 y 2 dxdy

2 2

2 ¡1 1+ey2

xdx =

xk =

donde wk son los pesos al integrar en una variable en el intervalo [¡1; 1]:

R1 R1

R1

R1

realizando un cambio de variables, y utilizando el polinomio de Legendre de segundo orden,

x ~k = (b¡a)xk2+(b+a) w ~k = (b¡a) 2 wk (d¡c)yk +(d+c) y~k = 2 w ~j0 = (d¡c) 2 wk

¡1

dxdy =

2 P

= 1 P

w ~j P (~ yk ) =

R1

2 1 e¡y dy ¡1 e¡y2 +1

=

2 P

w ~j P (~ yk ),

k=1

1 e¡y2 +1

wj P (yk ) =

k=0

k=1

= (w ~1 P (~ x1 ) + w ~2 P (~ x2 )) ¢

= w1 P (y1 ) + w2 P (y2 ) = = 0:8862269255 ¢ P (¡0:707106781) +

¢ (w ~1 P (~ y1 ) + w ~2 P (~ y2 )) = = (P (0:5773502692) + P (¡0:5773502692)) ¢

+0:8862269255 ¢ P (0:707106781) =

¢ (P (0:5773502692) + P (¡0:5773502692)) =

= 1: 103 3

= : 666 67 ¢ : 666 67 = : 444 45

El resultado de la aproximación numérica es,

Problema 72 Calcular una aproximación numérica de la integral Z 1Z 2 x dxdy 1 + ey2 ¡1 0

R1 R2

x ¡1 0 1+ey2

utilizando la evaluación de F (x; y) en 4 puntos.

22

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 F (x; y)dxdy ¼ F (x0 ; y0 )w0

1. Para 1 punto:

sea exacta para polinomios de grado 1 en x e y: Es decir P (x; y) = ax + by + c

=

R

T

-

x2 ydxdy = = F ( 23 ; 23 )Area(T ) = ¡ 2 ¢2 3

2 32

=

16 27

=

= : 592 59 Solución: 2. Para 3 puntos: R

Calculamos la integral de forma analítica: R 1 R 1¡x 0

0

(ax + by + c) dydx = 16 a + 16 b + 12 c

Igualamos el valor de la integral con la fórmula de integración numérica: 1 1 6a + 6b

a = c = 0; b = 1;

R

1 3

b = c = 0; a = 1; 16 a = w0 x0 a ! x0 = 13 ; luego para los valores w0 = de integración es exacta.

=

1 3 ; y0

=

1 3

la fórmula

8 15

8 15

=

Las propiedades que debe veri…car, para cumplir con la de…ción de norma, son:

= : 533 33

¯ ¯ 1 1 1 ¯ 1¯ El área del triángulo! Area(T ) = 2 ¯ 2 0 0 ¯ 0 2 0

¢ ¡ 4 4 12 4 4 12 = Area(T )[ 25 48 F ( 10 ; 10 ) + F ( 10 ; 10 ) + F ( 10 ; 10 ) ¡

Solución: En esta demostración vamos a generalizar para cualquier p. Al …nal particularizamos para p = 2 con el …n de hacer que la demostración sea más sencilla.

Utilizando las fórmulas de integración numérica:

F (x; y) = x2 y

=

Problema 75 Tomar N = 2 y demostrar que la norma k x k2 veri…ca las propiedades de la de…nición de norma q kxkp = p jx1 jp + jx2 jp

Solución: El cálculo de la integral de forma analítica nos da: R 2 R 2¡x ¡ 2 ¢ x y dydx = 0 0

2 - x ydxdy

ANALISIS NUMERICO MATRICIAL II

donde - es el triángulo de vértices (0; 0), (2; 0) y (0; 2) utilizando 1 punto, 3 puntos, y 4 puntos

x2 ydxdy =

=

= : 533 33

-

-

2 3

2 2 ¡ 27 48 F ( 3 ; 3 )] =

Problema 74 Calcular una aproximación numérica de la integral Z x2 ydxdy

R

¢1=

3. Para 4 puntos:

1 2

= w0 y0 b ! y0 =

1 2 ; x0

2 3

= : 666 67

Calculamos w0 ; x0 e y0 dando valores a a; b; c

1 6b

=

¡ ¢ = 13 Area(T ) F ( 22 ; 0) + F ( 22 ; 22 ) + F (0; 22 ) = =

+ 12 c = w0 (x0 a + y0 b + c)

a = b = 0; c = 1; 12 c = w0 c ! w0 =

2 - x ydxdy

1. kxkp = 0 () x = 0;

¯ ¯ ¯ ¯=2 ¯ ¯

p p jx1 jp + jx2 jp = 0 =) jx1 jp + jx2 jp = 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..

23

Solución:

2. k¸xkp = j¸j kxkp ; 8¸ 2 K y x 2 E; p k¸xkp = p j¸x1 jp + j¸x2 jp k¸xkp = k¸xkp =

Limp!1 kxkp = Limp!1

p p j¸jp jx1 jp + j¸jp jx2 jp

µq PN p

p

i=1 jxi j

¶

Extraemos el máximo componente de x, xmax : µq ¶ PN p p = Limp!1 i=1 jxi j

p p j¸jp (jx1 jp + jx2 jp )

µr ¶ P ³ jxi j ´p = Limp!1 p jxmax jp N = i=1 jxmax j

p k¸xkp = j¸j p jx1 jp + jx2 jp

r µ ¶ PN ³ i j ´p = Limp!1 jxmax j p i=1 jxjxmax = j

k¸xkp = j¸j kxkp ; c.q.d.

3. kx + ykp · kxkp + kykp ; 8x; y 2 E;

= jxmax j Limp!1

p p jx1 + y1 jp + jx2 + y2 jp · kxkp + kykp =)

= jxmax j Limp!1

=) jx1 + y1 jp + jx2 + y2 jp ·

Todos los elementos 1, con lo que

³p ´p p · p jx1 jp + jx2 jp + p jy1 jp + jy2 jp

Para p = 2 tenemos:

Limp!1

jx1 + y1 j2 + jx2 + y2 j2 ·

entonces

µq ¶2 q · jx1 j2 + jx2 j2 + jy1 j2 + jy2 j2 =)

³

jxi j jxmax j

jxmax j Limp!1

´p

µr PN ³ p i=1

³P ³ N i=1

jxi j jxmax j

=

½

³P ³ N i=1

jxi j jxmax j

jxi j jxmax j

´p ¶

=

´p ´1=p

son menores o iguales que

1 si xi = xmax ; 0 si xi 6= xmax

jxi j jxmax j

´p ´1=p

=

=) x21 + 2x1 y1 + y12 + x22 + 2x2 y2 + y22 ·

= jxmax j Limp!1 (0 + : : : + 0 + 1 + : : : + 1)1=p =

p p · x21 +x22 +2 (x21 + x22 ) (y12 + y22 )+y12 +y22 =)

= jxmax j, c.q.d.

=) x1 y1 + x2 y2 ·

p p (x21 + x22 ) (y12 + y22 ) =)

Problema 77 Tomar N = 2; y dibujar el lugar geométrico de los vectores x = (x1 ; x2 ) que veri…can

=) x21 y12 + 2x1 y1 x2 y2 + x22 y22 ·

1. kxk1 < 1

· x21 y12 + x21 y22 + x22 y12 + x22 y22 =)

2. kxk2 < 1 3. kxk1 < 1

=) 2x1 y1 x2 y2 · x21 y22 + x22 y12 =)

Solución: En las grá…cas 1, 2 y 3 se muestran los lugares geométricos de las normas 1, 2 e in…nito, respectivamente.

=) 0 · x21 y22 + 2x1 y1 x2 y2 + x22 y12 =)

1. kxk1 < 1 =) jxj + jyj < 1 =) y < 1 ¡ x

=) 0 · (x1 y2 + x2 y1 )2 ,

Esta ecuación representa, como borde, una recta de pendiente negativa. Tal y como se ve en la …gura 1, el lugar geométrico está contenido en un rombo. p ¡ ¢ 2. kxk2 < 1 =) (x2 + y2 ) < 1 =) x2 + y 2 < 1

que siempre es cierto, con lo que queda demostrado. Problema 76 Demostrar que Limp!1 kxkp = max jxi j

Esta es la ecuación de un círculo de radio menor que 1 y centro el origen. En la …gura 2 se muestra el lugar geométrico.

i

24

y

y

1

1

1

1

1

1

x

x

1

1

Figure 1: Lugar geométrico de kxk1

Figure 2: Lugar geométrico de kxk2

3. kxk1 < 1 =) max(x; y) < 1

Esto siempre es cierto ya que el producto de valores positivos siempre es positivo (o igual a cero si algún xi es cero).

Esto representa una recta de valor constante (x; y) menor que 1. En la …gura 3 se puede ver el lugar geométrico.

3. max(jx1 j ; jx2 j) · jx1 j + jx2 j : Es trivial (propiedad transitiva).

Problema 78 Tomar N = 2 y demostrar la siguiente desigualdad k x k1 ·k x k2 ·k x k1

De estas demostraciones se deduce que las distintas normas coinciden cuando el vector x descansa sobre uno de los ejes de coordenadas.

Solución: Esta desigualdad es equivalente a lo siguiente: p max(jx1 j ; jx2 j) · (x21 + x22 ) · jx1 j + jx2 j 1. max(jx1 j ; jx2 j) ·

p (x21 + x22 ) ()

() jxmax j · ()

x2max

·

Problema 79 Demostrar que si A; B son dos matrices de dimensión N xN; entonces para cualquier norma de matrices subordinada a una norma vectorial se veri…ca k AB k·k A k ¢ k B k Solución:

p (x21 + x22 ) ()

x21

+ x22

Esta desigualdad siempre es cierta ya que xmax es o bien x1 o bien x2 . p 2. (x21 + x22 ) · jx1 j + jx2 j ()

supx

kABxk kxk

supx

kABxk kBxk kBxk kxk

· supx

supx

kABxk kBxk kBxk kxk

· kBk ¢ kAk ;

= supx

kABxk kBxk kxk kBxk ; kBxk kxk

¢ supx

kABxk kBxk

entonces

¡ ¢ 2 () x21 + x22 · (jx1 j + jx2 j) ()

supx

() x21 + x22 · jx1 j2 + 2 jx1 j jx2 j + jx2 j2 ()

kABxk kxk

· kBk ¢ kAk ; c.q.d.

Problema 80 Demostrar que los autovalores de A son los ceros del polinomio característico P (¸):

() x21 + x22 · x21 + 2 jx1 j jx2 j + x22 () () 0 · 2 jx1 j jx2 j 25

8 8 < x1 + x2 = 0 < x1 = ¡x2 x1 + x2 = 0 : : 2x3 = 0 x3 = 0

y

0

1

p1 2 ¡1 p 2

x ¹1 = @

0

1 A

2. ¸2 ; ¸3 = 2 1

1

x

0

10 1 0 1 ¡1 1 0 x1 0 @ 1 ¡1 0 A @ x2 A = @ 0 A 0 0 0 x3 0

1

8 8 < ¡x1 + x2 = 0 < x1 = x2 x1 ¡ x2 = 0 : : x3 libre x3 libre

1 0 A;x ¹3 = @ 0 A x ¹2 = @ 1 0

Figure 3: Lugar geométrico de kxk1 Solución: De…nición de autovalor de una matriz A:

Axi = ¸i xi =) (A ¡ ¸i Id)xi = 0

0

p1 2 p1 2

0

0

1 0 0 A 1

Problema 82 Calcular las normas 2; 1 e in…nito de la matriz µ ¶ 1 0 A= 1 1

Problema 81 Calcular los autovectores de la matriz 0 1 1 1 0 @ 1 1 0 A 0 0 2

Solución:

y determinar una base ortonormal de R3 de autovectores de A:

1. kAk2 =

p ½ (AT A)

s µ ¶ q p 2 1 kAk2 = ½ = 3+12 5 = 1: 618 1 1

Solución: Calculamos los autovalores de A : jA ¡ ¸i Idj = 0

¯ ¯ ¯ 1¡¸ 1 0 ¯¯ ¯ ¯ 1 1¡¸ 0 ¯¯ = ¡4¸ + 4¸2 = 0 ¯ ¯ 0 0 2¡¸ ¯ ¸1 = 0 ¸2 = 2 ¸3 = 2

¯ ¯ 2¡¸ 1 ¯ ¯ 1 1¡¸

2. kAk1 = maxj

Calculamos los autovectores de A :

¯ ¯ ¯ = 1 ¡ 3¸ + ¸2 = 0, ¸ = ¯

³P 2

´ ja j ij i=1

kAk1 = max(1; 2) = 2

1. ¸1 = 0 0

p1 2 ¡1 p 2

1

contiene los autovectores de A que forman una base ortogonal en R3 :

jA ¡ ¸i Idj = 0 =) P (¸) = 0; c.q.d.

1

0

B=@

como xi 6= 0, entonces

10

p1 2 p1 2

La matriz,

xi 6= 0 2 E; ¸i 2 CÁAxi = ¸i xi

0

0

3. kAk1 = maxi

1

1 1 0 x1 0 @ 1 1 0 A @ x2 A = @ 0 A 0 0 2 x3 0

³P 2

kAk1 = max(2; 1) = 2 26

´

j=1 jaij j

=

3 2

§

1 2

p 5

0

1 2 ¡1 0 2. A = @ ¡1 2 ¡1 A 0 ¡1 2

Problema 83 Demostrar la siguiente igualdad: ½(At A) = ½(AAt ) Solución:

Solución: ° °El condicionamiento de una matriz Â(A) = kAk2 ¢ °A¡1 °2 : Calculamos los autovalores de ambas matrices:

Veamos que los polinomios característicos coinciden : ¡1

jAt A ¡ ¸i Idj = jAt j jAt A ¡ ¸i Idj jAt j = ¯ ¯ ¯ ¡1 ¡1 ¯ = ¯(At ) At AAt ¡ ¸i (At ) IdAt ¯ =

¯ ¯ ¯ 2¡¸ 2 ¡2 ¯¯ ¯ £ ¤ 1¡¸ 1 ¯¯ = (2 ¡ ¸) (1 ¡ ¸)2 ¡ 1 1. ¯¯ 2 ¯ ¡2 1 1¡¸ ¯ ¡2 [2(1 ¡ ¸) + 2] ¡2 [2 + 2(1 ¡ ¸)] = £ ¤ £ ¤ (2 ¡ ¸) ¸2 ¡ 2¸ ¡ 8 (2 ¡ ¸) = (2 ¡ ¸) ¸2 ¡ 2¸ + 8

¯ ¯ ¯ ¡1 ¯ = ¯AAt ¡ ¸i (At ) At ¯ =

de donde obtenemos

= jAAt ¡ ¸i Idj

¸1 = 2 ¸2 = ¡2 ¸3 = 4

Problema 84 Demostrar que si los autovectores de una matriz A de dimensión N xN forman una base ortonormal de RN ; entonces para la norma 2 se cumple:

Esta matriz es simétrica, luego posee una base ortonormal 3 de autovectores, con lo que el condicionamiento de A se puede calcular como:

° ° maxi fj¸i jg Â(A) = kAk2 ¢ °A¡1 °2 = mini fj¸i jg

Solución: Al ser una base de autovectores ortonormal, la norma kAk2 = ½(A) = maxi fj¸i jg

° ° 4 maxi fj ¸i jg = =2 Â(A) = kAk2 ¢ °A¡1 °2 = mini fj ¸i jg 2

Los autovalores de A¡1 vienen dados por:

¯ ¯ 2 ¡ ¸ ¡1 0 ¯ 2. ¯¯ ¡1 2 ¡ ¸ ¡1 ¯ 0 ¡1 2 ¡ ¸

Axi = ¸i xi =) =) A¡1 Axi = A¡1 ¸i xi =) =)

1 ¸i xi

¸1 = 2 p ¸1 = 2 + p2 ¸1 = 2 ¡ 2

= A¡1 xi =)

=) A¡1 xi = ¸0i xi ,

También es una matriz simétrica, con lo que sus autovectores forman una base ortonormal y su condicionamiento es:

donde ¸0i = ¸1i , 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 ) ½¯ ¯¾ ° ¡1 ° ¯ ¯ ©¯ ¯ª 1 °A ° = max ¯¸0i ¯ = max ¯ 1 ¯ = ; ¯ ¯ 2 i i ¸i mini fj¸i jg

p p ° ¡1 ° 2+ 2 maxi fj ¸i jg ° ° p Â(A) = kAk2 ¢ A = = 3+2 = 2 mini fj ¸i jg 2¡ 2

entonces

Problema 86 Sean las matrices A; R: Demostrar que la matriz A, y la matriz B = R¡1 AR poseen los mismos autovalores

° ° Â(A) = kAk2 ¢ °A¡1 °2 Â(A) = maxi fj¸i jg ¢ Â(A) =

maxi fj¸i jg mini fj¸i jg ;

¯ ¯ ¯ ¯ = 4 ¡ 10¸ + 6¸2 ¡ ¸3 = 0 ¯ ¯

1 mini fj¸i jg

Solución:

c.q.d.

Bxi = ¸i xi =) ¡ ¢ =) R¡1 AR xi = ¸i xi =)

Problema 85 Calcular el condicionamiento para la norma 2, de las siguientes matrices: 0 1 2 2 ¡2 1. A = @ 2 1 1 A ¡2 1 1

=) RR¡1 ARxi = R¸i xi =) 3 Vectores

ortonormales: dos vectores son ortonormales si cumplen ½ 0 si i 6= j T lo siguiente, xi xj = 1 si i = j 27

Problema 88 Demostrar trigonométricas

=) ARxi = ¸i Rxi =) =) Ayi = ¸i yi ,

siguientes

igualdades

q tan(®) = ¡ cot(2®) + sign(cot(2®)) 1 + cot2 (2®)

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.

¡ ¢ donde ® 2 ¡ ¼4 ; ¼4 , sign(x) = 1 si x ¸ 0 y sign(x) = ¡1 si x < 0;

Problema 87 Se considera la matriz µ ¶ 1 1 A= 1 1

1 cos ® = p 1 + tan2 (®) sin ® = tan(®) cos ® ¡ tan(®) + sin(2®) cot(2®) = 2 sin2 (®)

calcular el ángulo ® tal que la matriz µ ¶ cos ® sin ® R= ¡ sin ® cos ®

Solución:

veri…que que la matriz B = R¡1 AR sea diagonal.

1. cot(2®) =

Solución: Realizamos el cálculo de la matriz B : B = R¡1 AR = µ ¶µ ¶µ ¶ cos ® ¡ sin ® 1 1 cos ® sin ® = = sin ® cos ® 1 1 ¡ 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 r 1 cos ® = § 2 De esta igualdad se obtiene el valor del ángulo ® : ®=

las

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(®) =

p

¡2 cot(2®)§

4 cot2 (2®)+4

2

=

p 1 + cot2 (2®) p ½ ¡ cot(2®) + p1 + cot2 (2®) si ® ¸ 0 tan(®) = ¡ cot(2®) ¡ 1 + cot2 (2®) si ® < 0 = ¡ cot(2®) §

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 expresar la anterior igualdad de la siguiente forma: q tan(®) = ¡ cot(2®) + sign(cot(2®)) 1 + cot2 (2®)

¼ 3¼ ;® = 4 4

La matriz de rotación queda como sigue: p ¶ µ ¶ µ 1p sin ¼4 2 12 p2 cos ¼4 2 p = R1 = ¡ sin ¼4 cos ¼4 ¡ 21 2 12 2 p ¶ µ ¶ µ 1p cos 3¼ sin 3¼ ¡ 2 p2 12 p2 4 4 = R2 = ¡ sin 3¼ cos 3¼ ¡ 12 2 ¡ 12 2 4 4

2. cos ® = p

=

Calculamos los elementos de la diagonal: b1 = ¡2 cos ® sin ® + 1

r

1 1+tan2 (®)

=

1 cos2 (®)+sin2 (®) cos2 (®)

r

1 sin2 (®) 1+ cos 2 (®)

p cos2 (®) = cos ®

3. sin ® = tan(®) cos ® =

b1 = 0; b1 = 2 b2 = 2 cos ® sin ® + 1

=

b2 = 2; b2 = 0; luego las soluciones posibles son: µ ¶ µ ¶ 0 0 2 0 B1 = ; B2 = 0 2 0 0 28

sin(®) cos(®)

=

cos ® = sin ®

4. cot(2®) =

¡ tan(®)+sin(2®) 2 sin2 (®)

0

=

=

sin(®) ¡ cos(®) +2 sin(®) cos(®) 2 sin2 (®)

=

¡ sin(®)+2 sin(®) cos(®) cos(®) cos(®) 2 sin2 (®)

=

sin(®)(¡1+2 cos(®) cos(®)) 2 sin2 (®) cos(®)

=

cos(2®) sin(2®)

B B ¢B B @

=

0

= =

(2 cos2 (®)¡1) 2 sin(®) cos(®)

= 0

= cot(2®)

a0pq = 0 a0pp = app ¡ tan(®)apq a0qq = aqq + tan(®)apq = apj cos ® ¡ aqj sin ® j 6= p; q = apj sin ® + aqj cos ® j 6= p; q

0 0 1 0 0

De esta matriz se deducen las siguientes igualdades: qq ) a0pq = (app ¡a sin 2® + apq cos 2® 2 0 app = app cos2 ® + aqq sin2 ® ¡ apq sin 2® a0qq = app sin2 ® + aqq cos2 ® + apq sin 2® a0pj = apj cos ® ¡ aqj sin ® j 6= p; q aqj = apj sin ® + aqj cos ® j= 6 p; q

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.

En donde se iguala a0pq a cero para calcular el ángulo de rotación:

Vamos a expresar las matrices de la siguiente manera: 0 1 a11 ap1 ai1 aq1 an1 B ap1 app apj apq apn C B C C A=B B ai1 apj aij aqj ani C @ aq1 apq aqj aqq aqn A an1 apn ani aqn ann 0 1 1 0 0 0 0 B 0 cos ® : sin ® 0 C B C 0 : 1 : 0 C Rpq (®) = B B C @ 0 ¡ sin ® : cos ® 0 A 0 0 0 0 1 0 0 0 0 1

1 0 B 0 cos ® B 0 ¢B B 0 @ 0 ¡ sin ® 0 0

1 an1 apn C C ani C C¢ aqn A ann 1 0 0 sin ® 0 C C 0 0 C C= cos ® 0 A 0 1 aq1 apq aqj aqq aqn

ap1 sin ® + aq1 cos ® an1 (app ¡aqq ) sin 2® + apq cos 2® apn cos ® ¡ aqn sin ® 2 apj sin ® + aqj cos ® ani app sin2 ® + aqq cos2 ® + apq sin 2® apn sin ® + aqn cos ® apn sin ® + aqn cos ® ann

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:

0 0 0 ¡ sin ® 1 0 0 cos ® 0 0

ai1 apj aij aqj ani

ai1 apj cos ® ¡ aqj sin ® aij apj sin ® + aqj cos ® ani

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¡1 AR:

B = R¡1 AR = 0 1 0 B 0 cos ® B =B 0 B 0 @ 0 sin ® 0 0

ap1 app apj apq apn

a11 ap1 cos ® ¡ aq1 sin ® B ap1 cos ® ¡ aq1 sin ® app cos2 ® + aqq sin2 ® ¡ apq sin 2® B ai1 apj cos ® ¡ aqj sin ® =B B (app ¡aqq ) @ ap1 sin ® + aq1 cos ® sin 2® + apq cos 2® 2 an1 apn cos ® ¡ aqn sin ®

Problema 89 Dentro del método de Jacobi para el cálculo de autovalores demostrar las igualdades

a0pj a0qj

a11 ap1 ai1 aq1 an1

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 …nal. Para obtener a0pp y a0qq ; se opera de la siguiente manera: 1. a0pp = app cos2 ® + aqq sin2 ® ¡ apq sin 2® =

1

³ = app cos2 ® + app +

C C C¢ C A

2apq tan(2®)

¡apq sin 2® = app cos2 ®+ 29

´

sin2 ®¡

1 C C C C A

+

³

app sin(2®)+2apq cos(2®) 2 sin ® cos ®

´

sin2 ® ¡ apq sin 2® = a11

= app cos2 ®+ +

³

app 2 sin ® cos ®+2apq cos2 ®¡2apq sin2 ® 2 cos ®

´

sin ®¡

1 ® = ¡ arctan 2 = ¡: 553 57 2 = 2 ¡ tan(®) µ ¶ 1 a11 = 2 + tan arctan 2 = 2: 618 2

a33 = 1 + tan(®) 2

2

¡2apq sin ® cos ® = app cos ® + app sin ®+

a33 = 1 ¡ tan

+apq cos ® sin ® ¡ apq tan ® + apq sin ® cos ®¡

0

1 2: 618 0 0 A 1 0 B = R¡1 AR = @ 0 0 0 0: 381 97

2. a0qq = app sin2 ® + aqq cos2 ® + apq sin 2® = 2apq tan(2®)

´

2

Los autovalores son los elementos de la diagonal (2: 618; 1; 0: 381 97) : Como ® = ¡: 553 57; la matriz 0 1 0 cos ® 0 sin ® 0:850 65 0 ¡0:525 7 0 1 0 0 1 0 A=@ R (®) = @ 0:525 73 0 0:850 65 ¡ sin ® 0 cos ®

2

sin ® + aqq cos ®+

+apq sin 2® = =

³

aqq 2 sin ® cos ®¡2apq cos2 ®+2apq sin2 ® 2 sin ® cos ®

´

sin2 ®+

por tanto, en este caso, como con una única matriz de rotación conseguimos transformar A en una matriz diagonal, tendremos que los autovectores de A son simplemente los vectores columnas de R(®): Es decir el autovector del 0:850 65 0 autovalor 2: 618 es ( ), el autovector del auto0:525 73 0 valor 1 es ( 1 ); el autovector del autovalor 0: 381 97 es 0 ¡0:525 73 0 ( ): 0:850 65

+aqq cos2 ® + apq sin 2® = (aqq sin ® ¡ apq cos ®+ apq 2 + cos ® ¡apq cos ®) sin ®+aqq cos ®+apq sin 2® =

= aqq sin2 ® + aqq cos2 ® ¡ apq cos ® sin ®+ +apq tan ® ¡ apq cos ® sin ® + 2apq cos ® sin ® = = aqq + tan(®)apq Problema 90 Utilizar el método de Jacobi para aproximar los autovalores y autovectores de la matriz: 0 1 2 0 1 A=@ 0 1 0 A 1 0 1

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) µ ¶ 2 1 A= 0 1 µ ¶ ¡3 0 A= 1 1

Solución:

0

1 cos ® 0 sin ® 0 1 0 A R (®) = @ ¡ sin ® 0 cos ® (app ¡aqq ) 2

Solución: En este problema vamos a utilizar la norma euclídea aunque cualquier otra norma también sería válida. La norma in…nito, por ejemplo, simpli…caría los cálculos ya que es inmediato obtener el máximo de un vector.

sin 2® + apq cos 2® = 0 tan(2®) =

2apq (aqq ¡ app )

2 tan(2®) = (1 ¡ 2) ®=

¶ 1 arctan 2 = : 381 97 2

a21 = a32 = 0

¡2apq sin ® cos ® = app ¡ apq tan ®

³ = aqq ¡

µ

1. A =

arctan (¡2) 2

µ

2 1 0 1 1

¶

u2 = A kuu1 k = 30

=

µ

2 1 0 1

=

µ

3 2 p2 1 2 2

3

u =

4

¶Ã

p ¶

2 A kuu2 k

!

=

=

¶Ã

=

2 1 0 1

=

µ

7 10 p5p2 1 10 5 2

3

p !

3 p 2 2 5p 1 p 2 2 5

=

=

µ

2 1 0 1

¶µ

=

µ

3 2 p2 1 10 2

p

p ¶

7 10 p2 1 10 2

=

x¸ =

2. A =

µ

=

µ

p p p ¶ 27 ¡ 2132 1066 p p 26 p 2 , 2 533 1066 26 2

p ¶

¡3 0 1 1

=

µ

: 997 79 6: 651 9 £ 10¡2

¶

¶

¶Ã

=

µ

¡3 0 1 1

=

µ

p ¶ 3 ¡p 2 2 , 2

=

µ

p1 2 p1 2

!

¡ ¡®¢¢n sign u4 ; u3

1. ku1 k = 1 ! u2 =

µ

2 ¡1 ¡1 1

¶µ

1 1

¶

ku2 k = 1 ! u3 =

µ

2 ¡1 ¡1 1

¶µ

1 0

¶

=

=

¶µ

=

µ

1 0

¶

=

µ

2 ¡1

¶

Producto escalar (u2 ; u3 ) = 2 > 0. ! autovalor máximo = ku3 k = 2

° 2° 1 p °u ° = 26 = 2: 549 5 2

¡3 0 1 1

=

Solución:

1

2

¶

Problema 92 Calcular el autovalor µ mayor y¶el autovec2 ¡1 tor correspondiente de la matriz utilizando ¡1 1 el método de la potencia, dando 2 iteraciones del método a partir de u1 = (1; 1) y tomando como norma kuk = maxi jui j

u2 = A kuu1 k =

u3 = A kuu2 k =

9 1066 p26 p2 2132 p 1 ¡ 2132 1066 26 2

con signo positivo ya que (¡1)4 = 1.

El autovalor máximo aproximado es ¸ = 2: 126 y su autovector asociado es: p

p p

¡3 0 1 1

¸ = ¡3: 109 8; ¡®¢ con signo negativo ya que sign u4 ; u3 = ¡1 y su autovector asociado es à p p p ! µ ¶ p ¡ 213227¢82 1066 26 2 ¡: 958 8 65 026 p p p x¸ = = ; 2¢82 p : 284 09 1066 26 2 533 65 026

¶ ,

15 226 p113p2 1 226 113 2

¶µ

p

µ

El autovalor máximo aproximado es

° 4° 1 p °u ° = 113 = 2: 126 5

µ

p ° 3° °u ° = 1 1066 = 2: 511 5 13

p ° 4° °u ° = 1 65 026 = 3: 109 8 82

° 3° p °u ° = 5 = 2: 236 1

=

¶ ,

=

p p ¶ ,

3 A kuu3 k

p p

9 26 p2 26 p 1 ¡ 26 26 2

u4 = A kuu3 k =

° 2° p °u ° = 5 = 2: 236 1

=

µ

,

µ

u =

p1 2 p1 2

Autovector asociado normalizado

p p ¶ 3 ¡ 26 p 26p 2 = 1 13 26 2 31

u3 ku3 k

=

µ

1 ¡1=2

¶

0

1 0 5 1 ¡2 ¡1 0 ¡6 @ 0 1 ¡1 A u3 = 65 @ 32 A 0 0 ¡3 ¡ 13 0 1 1 30 ° ° A ; °u2 ° = 14 u3 = @ 14 15 15 2

Problema 93 Utilizar el método de la potencia inversa para aproximar el autovalor más pequeño de la matriz µ ¶ ¡2 1 A= 0 3 llegar hasta u3 partiendo de u = (1; 1):

15

El autovalor máximo (A®¢¡ 2Id)¡1 es ¸max = ¡- de 3 2 con signo positivo (sign u ; u = 1)

Solución: Aun = µ

un¡1 kun¡1 k

¡2 1 0 3

¶

u2 =

Ã

p1 2 p1 2

(A ¡ 2Id)¡1 x ¹ = ¸max x ¹

!

Para calcular el autovalor más cercano a 2, realizamos las siguientes operaciones:

p ¶ ° ° 1 ¡ 61p 2 ; °u2 ° = = : 333 33 1 3 6 2 µ ¶ µ 1p ¶ ¡2 1 ¡ 2p 2 u3 = 1 0 3 2 2 u2 =

u3 =

µ

µ

p ¶

1 3 p2 1 6 2

1 x ¹ = (A ¡ 2Id) x ¹ ¸max µ ¶ 1 A ¡ 2Id ¡ Id x ¹=0 ¸max ¶ ¶ µ µ 1 Id x ¹ = (A ¡ ¸prox Id) x ¹=0 A¡ 2+ ¸max µ ¶ µ ¶ 1 1 43 ¸prox = 2 + = 2 + 14 = ¸max 14 15

° ° 1p 10 = : 527 05 ; °u3 ° = 6

° 3° °u ° es el autovalor máximo de A¡1 , con lo que el p ¡1 6 autovalor mínimo de A es ¸min = ku 10 = ¡1: 3 k = ¡ 10 ¡- 3 2 ®¢ 897 4; con signo negativo ya que sign u ; u = ¡1.

¸prox = 3: 071 4

Problema 95 Escribir en Fortran las funciones siguientes: SIGNO_PRODUCTO_ESCALAR(uf,vf,Nf ) que devuelve el signo del producto escalar de los vectores uf y vf de dimensión Nf (12 líneas de código como máximo), y la función AUTOVALOR_MAXIMO(Af,uf,Nf,Nfmax,N…ter,Tolf) que devuelve el autovalor máximo de una matriz y su autovector por el método de la potencia. Los parámetros son la matriz Af, el vector candidato inicial uf, Nf la dimensión real, Nfmax, la dimensión para coger memoria, N…ter número máximo de iteraciones, y Tolf la tolerancia. Esta función devuelve el valor 2.**120 si P no termina correctamente. Tomar como norma kuk = i ABS(ui ) (28 líneas de código como máximo)

Problema 94 Calcular el autovalor y autovector más cercano a 2 de la matriz 0 1 0 ¡1 0 @ 0 3 ¡1 A 0 0 ¡1

para ello calcular dos iteraciones del método de la potencia inversa partiendo de u1 = (1; 1; 1): Solución: 0

1 ¡2 ¡1 0 1 ¡1 A A0 = A ¡ 2Id = @ 0 0 0 ¡3

Solución:

Vamos a utilizar la norma in…nito con el …n de simpli…car los cálculos.

01 02 03 04 05 06 07 08 09 10 11 12

n¡1

A0 un = kuun¡1 k 0 1 0 ¡2 ¡1 0 1 @ 0 1 ¡1 A u2 = @ 1 0 0 ¡3 1 0 5 1 ¡6 ° ° u2 = @ 23 A ; °u2 ° = ¡ 13

14 15

1 A 5 = : 833 33 6

32

SIGNO_PRODUCTO_ESCALAR(uf,vf,Nf) DIMENSION uf(*),vf(*) PE=0 DO 06 I=1,Nf PE=PE+uf(I)*vf(I) CONTINUE IF (PE.GT.0) THEN SIGNO_PRODUCTO_ESCALAR=1 ELSE SIGNO_PRODUCTO_ESCALAR=-1 ENDIF END

01 AUTOVALOR_MAXIMO(Af,uf,Nf,Nfmax,N…ter,Tolf) 02 DIMENSION Af(Nfmax,*),uf(*),vf(Nfmax) 03 DO 26 I=1,N…ter 04 NORMA=0 05 DO 11 J=1,Nf 06 vf(J)=0 07 DO 09 K=1,Nf 08 vf(J)=vf(J)+Af(J,K)*uf(K) 09 CONTINUE 10 NORMA=NORMA+ABS(vf(J)) 11 CONTINUE 12 IF (NORMA.EQ.0) THEN 13 AUTOVALOR_MAXIMO=2.**120 14 RETURN 15 ENDIF 16 AUTOVALOR_MAXIMO=SIGNO_PRODUCTO_ ESCALAR(uf,vf,Nf)*NORMA 17 ERROR=0 18 DO 22 J=1,Nf 19 vf(J)=vf(J)/AUTOVALOR_MAXIMO 20 ERROR=ERROR+ABS(uf(J)vf(J))/(ABS(vf(J)+1.) 21 uf(J)=vf(J) 1. 22 CONTINUE 23 IF(ERROR
un = Mun¡1 + c = 1 0 1 0 ¡1 0 1 0 = @ 12 0 0 A un¡1 + @ 23 A 1 0 13 0 3

Iteraciones: 1 0 1 0 1 0 ¡1 ¡1 0 u2 = M @ 0 A + @ 32 A = @ 32 A 1 1 0 3 3 0

u3 = M @ 0

¡1 3 2 1 3

1 2

1

0

A+@

1

0

5 6

¡1 3 2 1 3

¡1 3 2 1 3

1

1

un = M un¡1 + c

donde las matrices L; D; U son de la forma: 0 1 0 0 0 L = @ ¡1 0 0 A 0 ¡1 0 0 1 1 0 0 D=@ 0 2 0 A 0 0 3

0

A=@

1 2 5 6

0 7 4 2 3

1 A

80 19 < ¡1 = @ 0 A : ; 1

Problema 98 Calcular 3 iteraciones Gauss-Seidel para resolver el sistema 0 10 1 0 1 ¡1 0 x @ ¡1 2 0 A @ y A = @ 0 ¡1 3 z partiendo de u1 = (0; 0; 0)

33

1

una base 1 ortogonal de autovec0 1 2 0 A, 0 1

80 1 0 19 0 = < 1 0; @ 0 A ; @ 1 A $ 2 : ; 1 0

b

0

A=@ 1 A

1. Autovectores y autovalores:

Du = (¡L ¡ U)u + b (¡L ¡ U )u + D

La ecuación iterativa queda:

Solución:

Solución: La expresión matricial del método de Jacobi es como sigue: (L + D + U )u = b

u=D

Calculamos la matriz M y el vector c : 0 1 0 1 0 M = D¡1 (¡L ¡ U) = @ 12 0 0 A 0 13 0 0 1 0 1 ¡1 ¡1 c = D¡1 b = D¡1 @ 3 A = @ 32 A 1 1 3

Problema 97 Calcular 0 1 tores de la matriz @ 0 1

partiendo de u1 = (0; 0; 0)

¡1

1 0 ¡1 0 U =@ 0 0 0 A 0 0 0

3. u4 = M @ 1 A + @

Problema 96 Calcular 3 iteraciones del método de Jacobi para resolver el sistema 0 10 1 0 1 1 ¡1 0 x ¡1 @ ¡1 2 0 A @ y A = @ 3 A 0 ¡1 3 z 1

¡1

0

$

del método de 1 ¡1 3 A 1

Solución: La expresión matricial del método de GaussSeidel es de la forma:

Al ser A tridiagonal, el wopt se puede calcular como

(L + D)u = ¡U u + b

u = (L + D)¡1 (¡U )u + (L + D)¡1 b

de:

un = M un¡1 + c

MJ es la matriz del método de Jacobi que se obtiene MJ = D¡1 (¡L ¡ U ) = 10 1 0 1 0 0 0 1 0 = @ 0 12 0 A @ 1 0 0 A = 0 0 13 0 1 0 1 0 0 1 0 = @ 12 0 0 A 0 13 0

Si construimos el sistema de ecuaciones y despejamos las incógnitas: 8 8 < x ¡ y = ¡1 < x = ¡1 + y y = 3+x ¡x + 2y = 3 2 : : ¡y + 3z = 1 z = 1+y 3 Iteraciones:

Los autovalores de MJ : ½ (MJ )2 = 12

1. x = ¡1 y = 3¡1 2 =1 2 z = 1+1 3 = 3

wopt =

2. x = 0 y = 32 z = 65 3. x = ¡1 + 32 = 3+ 1 y = 2 2 = 47 1+ 7 z = 3 4 = 11 12

2 q 1 + 1 ¡ ½ (MJ )2

wopt =

(L + D + U )u = b

wopt =

=

p2

1+

1¡ 12

4p 2+ 2

Iteraciones del sistema: xn = w (yn¡1 ¡ 1) + (1 ¡ w) xn¡1 n yn = w 3+x + (1 ¡ w) yn¡1 2

Problema 99 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.

n zn = w 1+y + (1 ¡ w) zn¡1 3 1 0 1 ¡w ¡1: 171 6 A = @ 1: 071 1 A u2 = @ w 3¡w 2 : 808 83 w 1+1:3071 1 0 w (1: 071 1 ¡ 1) ¡ (1 ¡ w) 1: 171 6 35 + (1 ¡ w) 1: 071 1 u3 = @ w 3+: 284 2 1+1: 740 2 w + (1 ¡ w) : 808 83 3 0 1 : 284 35 = @ 1: 740 2 A : 931 34 0 w (1: 740 2 ¡ 1) + (1 ¡ w) : 284 35 42 + (1 ¡ w) 1: 740 2 u4 = @ w 3+: 818 2 1+1: 938 2 w + (1 ¡ w) : 931 34 3 0 1 : 818 42 = @ 1: 938 2 A : 987 65

0

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 100 Calcular 3 iteraciones del método de relajación para resolver el sistema 0 10 1 0 1 1 ¡1 0 x ¡1 @ ¡1 2 0 A @ y A = @ 3 A ; 0 ¡1 3 z 1

partiendo de u1 = (0; 0; 0): Calcular previamente el parámetro de relajación óptimo

Cálculo del wopt :

2p 1+ 12 2

=

wopt = 1: 171 6

1 2

Solución:

2 1¡½(MJ )2

p

1+

p ¢ ¡ 1p 0; 2 2; ¡ 12 2 , luego

1

A=

1

A=

Problema 101 Escribir en Fortran la función siguiente: CONDICIONAMIENTO(Af,Nf,Nfmax,TOLf,N…ter) que devuelve el condicionamiento de una matriz utilizando el método de Jacobi para calcular los autovalores, se supondrá implementada la función JACOBI(A,N,Nmax,TOL,Niter)

x ¡ y = ¡1 ¡x + 2y = 3 ¡y + 3z = 1

34

2. Para demostrar que jAj = 0 cuando la por ¯ suma ¯ columnas es cero, basta saber que jAj = ¯AT ¯ : Pn i aij = 0, esto es equivalente a lo siguiente: 0 1 0 1 0 1 1 0 1 B 1 C B 0 C B 1 C B C B C B C AT B . C = B . C = 0 ¢ B . C @ .. A @ .. A @ .. A

que devuelve 0 si termina bien y 1 si termina mal. La función CONDICIONAMIENTO devuelve 2.*120 si termina mal porque Jacobi da un error o se produce una división por cero. Los parámetros son la matriz Af, Nf la dimensión real, Nfmax, la dimensión para coger memoria, N…ter número máximo de iteraciones, y Tolf la tolerancia (máximo .21 líneas de instrucciones) Solución:

1

01 CONDICIONAMIENTO(Af,Nf,Nfmax,N…ter,Tolf) 02 DIMENSION Af(Nfmax,*) 03 CONDICIONAMIENTO=2.*120 04 IF(JACOBI(Af,Nf,Nfmax,TOLf,N…ter).NE.0) THEN 05 RETURN 06 END IF 07 AMAX=ABS(A(1,1)) 08 AMIN=ABS(A(1,1)) 09 DO I=2,Nf 10 IF(AMAX.LT.ABS(A(I,I))) THEN 11 AMAX=ABS(A(I,I)) 12 END IF 13 IF(AMIN.GT.ABS(A(I,I))) THEN 14 AMIN=ABS(A(I,I)) 15 ENDIF 16 END DO 17 IF(AMIN.EQ.0) THEN 18 RETURN 19 ENDIF 20 CONDICIONAMIENTO=AMAX/AMIN 21 END

Problema 103 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 de M correspondientes a autovalores de módulo menor que 1; entonces el método converge. Solución: Sean xi los autovectores de M correspondientes a autovalores menores que 1: P u1 = ni=1 ai xi P c = ni=1 ci xi

Realizando iteraciones obtenemos las siguientes expresiones: u2 = Mu1 + c

¡ ¢ u3 = Mu2 + c = M M u1 + c + c = M 2 u1 + M c + c .. . un = M n¡1 u1 + M n¡2 c + : : : Mc + c = ¡ ¢ = M n¡1 u1 + M n¡2 + : : : M + 1 c

Solución: Si jAj = 0, entonces la matriz A no es invertible y el sistema no tiene solución. 1. Vamos a demostrar que si la suma por …las de A es igual a cero, entonces jAj = 0 : Pn j aij = 0, esto es equivalente a lo siguiente: 0 1 0 1 0 1 1 0 1 B 1 C B 0 C B 1 C B C B C B C AB . C = B . C = 0 ¢ B . C . . @ . A @ . A @ .. A 1 0 1

Tomando el primer sumando: Pn Pn M n¡1 u1 = M n¡1 i=1 ai xi = M n¡2 i=1 ai M xi = P P = M n¡2 ni=1 ai ¸i xi = : : : ni=1 ai ¸n¡1 xi i

Como u1 depende linealmente de los xi (autovectores) cuyos autovalores ¸i son menores que uno, entonces ¸n¡1 i tiende a 0 cuando n tiende a in…nito, luego este término converge.

Esto signi…ca que la matriz A posee un autovalor igual a cero (¸ = 0).

Para el segundo sumando: ¡ n¡2 ¢ M + :::M + 1 c = ¡ ¢ Pn = M n¡2 + : : : M + 1 i=1 ci xi = P Pn n = M n¡2 i=1 ci xi + : : : M i=1 ci xi +

El determinante de una matriz es igual al producto de sus autovalores: n Y

i=1

1

T La matriz ¯ T ¯A posee un autovalor igual a cero (¸ = 0), ¯ ¯ luego A = 0. ¯ ¯ jAj = ¯AT ¯ = 0

Problema 102 Demostrar que si una matriz A veri…ca que por …las o columnas su suma es siempre igual a 0, entonces el determinante de A es cero, y por tanto el sistema asociado a A no tiene solución.

jAj =

0

¸i = ¸1 ¢ : : : 0 : : : ¢ ¸n = 0 35

+

Pn

i=1 ci xi

=

Pn

Pn

n¡2 xi + i=1 ci ¸i

Problema 105 Plantear el algoritmo necesario para calcular, utilizando el método de Newton-Raphson, las raíces complejas o reales de un polinomio de grado 3:

Pn

+ : : : i=1 ci ¸i xi + i=1 ci xi = ¡ ¢ Pn + : : : + ¸i + 1 · = i=1 ci xi ¸n¡2 i | {z }

Solución:

Serie geométrica convergente

·

Pn

1 i=1 ci xi 1¡¸i ,

P (z) = az 3 + bz 2 + cz + d = 0

con lo que este término también converge.

Un polinomio de grado 3 posee al menos una raíz real. Las otras dos raíces pueden ser también reales o imaginarias conjugadas. Sea z un número complejo: z = x + yi, sustituyendo en la anterior ecuación,

Problema 104 Calcular 2 iteraciones del método de Newton-Raphson no-lineal para aproximar una raíz del sistema de ecuaciones x2 + y 2 ¡ 1 = 0 y¡x =0

P (x + yi) = a (x + yi)3 + b (x + yi)2 + c (x + yi) + d

partiendo de (x; y) = (1; 1):

P (x + yi) = ax3 + 3iax2 y ¡ 3axy 2 ¡ iay 3 + +bx2 + 2ibxy ¡ by2 + cx + icy + d = 0

Solución: ½ 2 ½ µ ¶ x + y2 ¡ 1 = 0 2x 2y rf(x; y) = y¡x =0 ¡1 1

Separamos la parte real de la parte imaginaria: ½ ax3 ¡ 3axy 2 + bx2 ¡ by 2 + cx + d = 0 f= 3ax2 y ¡ ay3 + 2bxy + cy = 0 µ 3ax2 ¡ 3ay 2 + 2bx + c ¡6axy ¡ 2by rf = 6axy + 2by 3ax2 ¡ 3ay2 + 2bx + c

un = (xn ; y n ) u0 = (1; 1) ½ rf (un )z = ¡f (un ) un+1 = un + z µ ¶µ ¶ µ 2 ¶ 2x 2y z1 x + y2 ¡ 1 =¡ ¡1 1 z2 y¡x

El proceso iterativo es de la forma: un = (xn ; yn ) ½ rf (un ) ¢ z = ¡f (un ) un+1 = un + z

Iteraciones: µ ¶µ ¶ µ ¶ 2 2 z1 1 1. =¡ ¡1 1 z2 0 µ

¶

z1 z2

1

=

µ

¡ 14 ¡ 14

0

u =u +z =

u1 =

2.

µ

µ

µ

3 2

¡1 z1 z2

3 4 3 4 3 2

1 ¶

=

¶ ¶µ µ

u =

µ

17 24 17 24

¶

1 1

µ

=

¶

¶

z1 z2

+

3 4 3 4

µ

µ

=¡

¶

1 ¡ 24 1 ¡ 24

u2 = u1 + z =

2

µ

¶

¶

+

µ

: 708 33 : 708 33

1 8

0

Funcion F (u) f (1) = a ¢ u(1)3 ¡ 3a ¢ u(1) ¢ u(2)2 + +b ¢ u(1)2 ¡ b ¢ u(2)2 + c ¢ u(1) + d f (2) = 3a ¢ u(1)2 ¢ u(2) ¡ a ¢ u(2)3 + +2b ¢ u(1) ¢ u(2) + c ¢ u(2) devolver f Fin funcion

¶

1 ¡ 24 1 ¡ 24

¶

Las funciones F (u) y rF (u) se utilizan para evaluar la función y el gradiente de la función en un punto, respectivamente.

¶

¡ 41 ¡ 14

µ

Algoritmo: Este algoritmo utiliza una función, "Sistema(A; u)"; para resolver un sistema de ecuaciones.

Funcion rF (u) 2 rf (1; 1) = 3a ¢ u (1) ¡ 3a ¢ u(2)2 + 2b ¢ u (1) + c rf (1; 2) = ¡6a ¢ u (1) ¢ u (2) ¡ 2b ¢ u (2) rf (2; 1) = ¡rf (1; 2) rf (2; 2) = rf (1; 1) devolver rf Fin funcion

¶

36

Solución:

Algoritmo un¡1 = (x0 ; y0 ) /* calculamos¡ la primera ¢¢ */ ¢ aproximación ¡ ¡ z = Sistema rF un¡1 ; ¡F un¡1 un (1) = un¡1 (1) + z(1) un (2) = un¡1 (2) + z(2) n=0 ¯ ¡¯ ¢ Mientras ¯un ¡ un¡1 ¯ ¸ T OL y (n < T OP ) un¡1 = un /* calculamos¡ la siguiente ¢¢ */ ¢ aproximación ¡ ¡ z = Sistema rF un¡1 ; ¡F un¡1 un (1) = un¡1 (1) + z(1) un (2) = un¡1 (2) + z(2) n=n+1 Fin Mientras Si (n = T OP ) Entonces ERROR: No se ha encontrado solución Fin Si Fin Algoritmo

0

yzexyz @ 0 rf (x; y; z) = 4 (z ¡ 1) x3 ½

0

un+1 = un + z = 0 1 0 1 1 ¡ e (e ¡ 1) ¡ 11 =@ 1 A+@ 2 1 3 0 15 1 1 ¡ 2 ¡ e (e ¡ 1) 13 A =@ 2 4

Solución: µ

¶ y x¡1 1. rf (x; y) = ! rf(1; 1) = y¡2 x µ ¶ µ ¶¡1 µ ¶ 1 0 1 0 0 ! u2 = + ¡1 1 ¡1 1 1 µ ¶ µ ¶ 1 1 = 1 2 2 0 rf (1; 2) = µ ¶ µ 0¶ 1 1 1 = 2 2

! u3 =

2 0 0 1

¶¡1 µ

0 0

¶

10 1 z1 xyexyz ¡3z 2 A @ z2 A = z3 x4 1

Iteración: 0 10 1 0 1 e e e z1 e¡1 @ 0 2 ¡3 A @ z2 A = ¡ @ ¡2 A 0 0 1 ¡3 z3 0 1 0 1 1 z1 ¡ e (e ¡ 1) ¡ 17 2 11 @ z2 A = @ A 2 z3 3

A partir de u1 = (1; 1), calcular u2 y u3 utilizando el método de Newton-Raphson para aproximar un cero del sistema no-lineal.

µ

xzexyz 2y 0

exyz ¡ 1 2 @ y ¡ z3 ¡ 2 A =¡ (z ¡ 1)x4 ¡ 3

(x ¡ 1)y = 0 (y ¡ 2)x = 0

¶

rf (x; y; z) z = ¡f (x; y; z) un+1 = un + z

yzexyz @ 0 4 (z ¡ 1) x3 0

Problema 106 Se considera el sistema no-lineal

µ

1 xyexyz ¡3z 2 A x4

xzexyz 2y 0

17 2

1

A=

INTERPOLACION DE FUNCIONES II Problema 108 Calcular los polinomios base de Hermite que corresponden a tomar como puntos de interpolación x0 = ¡1; x1 = 1; y el orden de derivación M = 1:

+

Solución: Los polinomios de Hermite que corresponden a esos puntos de interpolación vienen dados por las grá…cas 4, 5, 6 y 7.

Problema 107 Calcular 1 iteración del método de Newton-Raphson no-lineal para aproximar una raíz del sistema de ecuaciones

1. La grá…ca 4 se hace cero en 1 y sus derivadas, tanto en ese punto como en ¡1, valen cero. Este polinomio tiene dos raíces en 1 (la segunda debido al valor de su derivada en 1), con lo que la forma de este polinomio es como sigue:

exyz ¡ 1 = 0 y2 ¡ z 3 ¡ 2 = 0 (z ¡ 1)x4 ¡ 3 = 0 partiendo de (x; y; z) = (1; 1; 1):

2

0 H¡1 (x) = (x ¡ 1) (a (x + 1) + b)

El valor de este polinomio en -1 es 1: 37

y

y

1

1

0.75

0.75

0.5

0.5

0.25

0.25

0 -1

-0.5

0 0

0.5

1

-1

-0.5

0

0.5

x

0 Figure 4: Polinomio de Hermite H¡1

Figure 6: Polinomio de Hermite H10 -1

-0.5

0

y

0.5

x 1

0 0.25

-0.05

0.2

-0.1

0.15

-0.15

0.1

-0.2

0.05

-0.25

0 -1

1 x

-0.5

y 0

0.5

1 x

Figure 7: Polinomio de Hermite H11

1 Figure 5: Polinomio de Hermite H¡1

Por la misma razón que en el caso anterior, sabemos que la función posee dos raíces en 1, con lo que el polinomio tiene la forma,

0 H¡1 (¡1) = 1

(¡1 ¡ 1)2 (a (¡1 + 1) + b) = 4b = 1 b=

1 H¡1 (x) = (x ¡ 1)2 (a (x + 1) + b)

1 4

1 H¡1 (¡1) = (¡1 ¡ 1)2 (a (¡1 + 1) + b) = 4b = 0

Al ser la derivada en -1 igual a cero tenemos: 2

00 H¡1 (x) = 2 (x ¡ 1) (a (x + 1) + b) + (x ¡ 1) a = 0

b=0

2

para calcular el valor de a, derivamos el polinomio y evaluamos en 1,

00 H¡1 (x) = 2 (¡2) (a (0) + b) + (¡2) a = 0

¡4b + 4a = 0

10 H¡1 (x) = 2 (x ¡ 1) (a (x + 1) + b) + (x ¡ 1)2 a

a = b = 14 ;

10 H¡1 (¡1) = 2 (¡1 ¡ 1) (a (¡1 + 1) + b) +

luego el polinomio queda, 0 H¡1 (x) =

1 4

+ (¡1 ¡ 1)2 a = 1

(x ¡ 1)2 (x + 2)

4a = 1

2. Para calcular el segundo polinomio partimos de la grá…ca 5. En ésta, La función se anula en -1 y 1, la derivada en -1 es igual a 1 y su derivada en 1 es cero.

a = 14 ;

38

µ

luego el polinomio nos queda: 1 H¡1 (x) = (x ¡ 1)2

¡1

4

¢ (x + 1)

di =

3. Para calcular los otros dos polinomios, basta considerar que son funciones simétricas a las dos anteriores. En la grá…ca 6 se puede ver que esta función 0 es simétrica a H¡1 (x) (ver grá…ca 4) con respecto al eje de las y:

=

bi =

¶

2 5 ¡ 85

ci+1 ¡ci 3hi

i = 0; : : : N ¡ 1

(ai+1 ¡ai ) hi

2 5

3 ¡8

2 5 ¡5

8 5

3

3

¡

1

0 2 1 C @ 152 A ¡3 A= 8 15

hi (2ci +ci+1 ) 3

i = 0; : : : N ¡ 1

1 0 1 0 2 1 ¡ 15 b0 @ b1 A = @ 1 ¡ 45 ¡ 85 A = @ 3 b2 ¡1 + 16 15 0

2

(¡x ¡ 1) (¡x + 2)

13 15 19 15 1 15

1 A

Los splines cúbicos nos quedan de la siguiente manera:

2 H10 (x) = ¡ 14 (x + 1) (x ¡ 2)

P1 (x) =

4. Por último, la función representada en la grá…ca 7, es 1 simétrica al polinomio H¡1 (grá…ca 5) con respecto al origen, con lo que,

2 15

(x + 1)3 +

13 15

(x + 1) ¡ 1

x 2 [¡1; 0] P2 (x) = ¡ 23 x3 + 52 x2 +

1 H11 (x) = ¡H¡1 (¡x)

H11 (x) = ¡ (¡x ¡ 1)2

µ

1 0 d0 @ d1 A = B @ d2

0 H10 (x) = H¡1 (¡x) 1 4

=

0

El polinomio es por tanto,

H10 (x)

¶

c1 c2

19 15 x

x 2 [0; 1] ¡1 4

H11 (x) =

1 4

¢ (¡x + 1)

P3 (x) =

8 15

(x ¡ 1)3 ¡

8 5

(x ¡ 1)2 +

1 15

(x ¡ 1) + 1

x 2 [1; 2]

(x + 1)2 (x ¡ 1)

Problema 109 Calcular los polinomios que determinan la interpolación por splines cúbicos de la función f(x) = ¡ ¢ sin ¼2 x para los puntos x = ¡1; 0; 1; 2

y

1

0.5

0

Solución: Los polinomios son de la forma:

-1

-0.5

0

0.5

1

1.5

2 x

P (x) = dx3 + cx2 + bx + a

-0.5

Vamos a calcular los coe…cientes para cada intervalo:

-1

8 (xi ¡ xi¡1 )

hi = 1 ai = f (xi )

Figure 8: Comparación entre la función sin aproximación por splines cúbicos.

i = 0; : : : N

0

1 0 a0 f (x0 ) B a1 C B f (x1 ) B C B @ a2 A = @ f (x2 ) a3 f (x3 )

1

0

1 ¡1 C B 0 C C=B C A @ 1 A 0

Problema 110 Calcular la función que interpola, utilizando la función sinc(x) a la función f(x) = sin(x) en los puntos x = ¡¼; ¡ ¼2 ; 0; ¼2 ; ¼:

hi¡1 ci¡1 + 2 (hi¡1 + hi ) ci + hi ci+1 = =

3(ai+1 ¡ai ) hi

¡

¡¼ ¢ 2 x y su

3(ai ¡ai¡1 ) hi¡1

Solución:

c0 = c3 = 0 µ ¶µ ¶ µ ¶ 4 1 c1 0 = 1 4 c2 ¡6

sinc (x) =

sin(x) x

La interpolación a través de la función sinc(x) : 39

P sin(¼( a ¡i)) f~(x) ¼ N i=M f(xi ) ¼( x ¡i) a x

¼ 2

xi = a ¢ i =

[¡2; ¡1; 0; 1; 2]

sin(¼( +2)) f~(x) ¼ f (¡¼) ¼ 2x¼+2 +f (¼ ) 2x

+f (0)

sin(¼( 2x ¼ )) ¼( 2x ¼ )

+ sin =

+f

sin(¼( 2x ¼ ¡2))

+f (¼)

¡ ¼ ¢ sin(2x¡¼) 2

2x¡¼

¡

sin 2x 2x¡¼

=

¡ ¡¼ ¢ sin(¼( 2x ¼ +1)) + 2 ¼( 2x ¼ +1)

¡ ¡¼ ¢ sin(2x+¼) 2

sin(2x¡¼) 2x¡¼

2x+¼

La interpolación por polinomios trigonométricos tiene la forma:

+

sin(2x+¼) 2x+¼

=

sin 2x(2x¡¼)¡sin 2x(2x+¼) 4x2 ¡¼2

=

=

Solución: ½ ¡x si ¡¼ · x · 0 jxj = x si 0 · x · ¼

¡ ¼ ¢ sin(¼( 2x ¼ ¡1)) + 2 ¼( 2x ¼ ¡1)

= sin

¼( 2x ¼ ¡2)

sin 2x 2x+¼

Problema 111 Calcular el polinomio trigonométrico tomando N = 2; que interpola a la función f (x) = jxj en el intervalo [¡¼; ¼]:

¡

~ ¼ f(x)

En la …gura 9 se muestran el sin(x) y su aproximación por el seno cardinal

ck =

1

R¼ f (x)eikx dx jxjeikx dx ¡¼ = 2¼ 2¼ R0 R ¼ ikx ikx ¡ ¡¼ xe dx xe dx + 0 2¼ 2¼

¡¼

c2 = c¡2 = 0 c1 = c¡1 = ¡2 ¼ c0 = ¼2

0 0

5

10 x

Sustituimos en el sumatorio que aproxima a la función y obtenemos:

-0.5

f~(x) ¼

-1

¡2 ¡ix ¼ e

= 12 ¼ ¡

4 ¼

+

¼ 2

¡ ¼2 eix =

cos x

El resultado de la aproximación es, por tanto,

Figure 9: Comparación del sin x con su aproximación numérica utilizando sin c(x), tomando como puntos de in¼ terpolación x=¡¼; ¡¼ 2 ; 0; 2 ; ¼:

y

=

Los valores de estos coe…cientes son:

0.5

-5

k=¡2

R¼

=

-10

ck eikx ;

donde los coe…cientes se calculan a partir de la siguiente expresión:

2x = ¡2¼ 4xsin 2 ¡¼ 2

y

2 X

1 4 f~(x) ¼ ¼ ¡ cos x 2 ¼ La siguiente grá…ca compara f(x) = jxj con su aproximación f~(x) para N = 2 en el intervalo [¡¼; ¼]: jxj

1

y

5

0.5

3.75

0 -10

-5

0

5

10 x 2.5

-0.5

1.25

-1

0 -2.5

-1.25

0

1.25

2.5 x

Figure 10: Comparación del sin x con su aproximación numérica utilizando sin c(x), tomando como puntos de in¡¼ ¼ 3¼ terpolación x=¡2¼; ¡3¼ 2 ; ¡¼; 2 ; 0; 2 ; ¼; 2 ; 2¼

Polinomio trigonométrico (N = 2; [¡¼; ¼]) En la siguiente grá…ca se realiza la misma comparación tomando 20 muestras en el intervalo [¡¼; ¼]:

40

jxj y

5

3.75

2.5

1.25

0 -2.5

-1.25

0

1.25

2.5 x

Polinomio trigonométrico (N = 10; [¡¼; ¼])

Problema 112 Calcular la cuadrática lineal de la tabla xi 0 1 2 3

aproximación

mínimo

yi 0 1 0 2

Solución: Aplicando las fórmulas para calcular los coe…cientes de la recta que más se aproxima a estos puntos, obtenemos: a=

PN

N

P PN xi yi ¡ N i=1 xi i=1 yi PN P 2 N i=1 x2i ¡( N i=1 xi )

4(1+6)¡(1+2+3)(1+2) 4(1+22 +32 )¡(1+2+3)2

= b=

=

i=1

PN

i=1

x2i N

=

1 2

PN P yi ¡ N i=1 xi i=1 xi yi P 2 N 2 x ¡ x ( ) i i=1 i i=1

PN

Pi=1 N

(1+22 +32 )(1+2)¡(1+6)(1+2+3)

=

4(1+22 +32 )¡(1+2+3)2

= =0

P (x) = ax + b = 21 x y

4

3

2

1

0 0

1

2

3

4 x

Figure 11: Aproximación mínimo cuadrática

41