Cosa

5
Tecnológico de Monterrey Física Computacional 2 Dr. Servando López Aguayo Tarea # 1 Número del equipo 3 Nombres Matrícula José Jame Paredes Peralta A01099631 Adriana E. Amaro A01139083 María Martínez Casales A01175585 Fecha de entrega: 19/08/2015 Calificación: _____ de ______ puntos

description

cosa

Transcript of Cosa

Page 1: Cosa

!Tecnológico!de!Monterrey!

!Física!Computacional!2!

!Dr.!Servando!López!Aguayo!

!Tarea!#!1!

!Número!del!equipo!3!

!!Nombres! Matrícula!José!Jame!Paredes!Peralta! A01099631!Adriana!E.!Amaro! A01139083!María!Martínez!Casales! A01175585!!!!Fecha!de!entrega:!19/08/2015!

Calificación:**_____**de*______*puntos*

Page 2: Cosa

1. Errores Computacionales

El programa tiene como proposito calcular los termi-nos de la serie de Taylor para la funcion exponencialy evaluarlos en dos valores de x. Se analizara el errorque existe al utilizar la serie respecto al que utiliza Matlab.

Codigo

%n-avo termino de la serien = 99;%Valores a evaluar de xx = [50 -50];l = length(x);%Valor real para calcular el errorvReal = exp(x);

%Calcular n-avos terminos de la serie%Matriz para almacenar n-avo termino para cada xterminos = zeros(l,n+1);%Matriz para almacenar estimacion utilizando el%n-avo terminovExp = zeros(l,n+1);for k = 1:lfor i = 0:n%n-avo termino con x(k)terminos(k,i+1) = x(k)^i/factorial(i);%Calcular n-ava estimacion para x(k)vExp(k,i+1) = sum(terminos(k,1:i+1));endend

%Calcular erroreserrorAbsoluto = zeros(l,n+1);errorRelativo = zeros(l,n+1);for k = 1:lfor i = 0:n\begin{itemize}\item\end{itemize}errorAbsoluto(k,i+1) = abs(vReal(k)-vExp(k,i+1));errorRelativo(k,i+1) = abs(errorAbsoluto(k,i+1)...

/vReal(k));endend

Graficas �! Ver figura 1

a ¿Que diferencia encuentras para la aproximacion de

x=+50 y x=-50? ¿Por que ocurre este comporta-

miento?

La aproximacion para x = 50 tiene un comporta-miento esperado (decreciente) conforme se agreganterminos a la serie, Al principio no parece mejorar,sin embargo alrededor del 40o termino comienza adecrecer notablemente. Por otro lado, para x = 50 elerror pasa de ser menor a mayor en la primera mitadde los terminos y comienza a decrecer despues.

Se puede notar que la aproximacion para este valorno es nada buena debido a que el logaritmo del errorse encuentra entre los ordenes de 1020 y 1045. Estainexactitud en el calculo se puede deber al proceso de

n0 20 40 60 80 100

Log(ErrorRelativo)

10-10

10-9

10-8

10-7

10-6

10-5

10-4

10-3

10-2

10-1

100Error Relativo con x = 50

n0 20 40 60 80 100

Log(ErrorRelativo)

1020

1025

1030

1035

1040

1045Error Relativo con x = -50

Figura 1: Errores relativos en escala logarıtmica

resta que ocurre debido al signo negativo del numeroevaluado. Con x = 50 no existe tal comportamientodebido a que todos los terminos que se suman sonpositivos, por otro lado la serie con x = �50 esalternante. Esto, computacionalmente hablando,puede afectar el resultado.

b ¿Si no pudieras usar la funcion exp, o revisar tablas

de valores numericos. ¿Como pudieras determinar el

posible error en dicha funcion? ¿Como encontrarıas

el rango de valores de n que da resultados optimos?

En primera instancia se podrıa generar una mejoraproximacion para cada valor si se expandiera la se-rie alrededor de ese mismo punto. Esta serıa nuestrareferencia, confiable hasta cierto error de truncamien-to. El rango de valores de “n” que darıan resultadosoptimos serıan aquellos que hagan que la suma de elresultado mas cercano a la referencia propuesta.

2. Derivadas Numericas Se utilizan 3 diferentes esquemaspara hacer una aproximacion numerica de la derivada de lafuncion seno cardinal (sin(x)/x). Se definio el seno cardinalcomo sin(x)/x para x 6= 0 y 1 para x = 0.A pesar deque se conoce la derivada analıtica exacta, tener metodosnumericos para aproximar derivadas es util si solo se tienenpuntos discretos sobre la funcion o si es mas sencillo evaluarnumeros en la funcion original que en la misma derivada.

Codigo

D = [5:-1:1 1:-1.1:.1 .1:-.01:.001];%IntervaloI = 30;x = -10*pi:D:10*pi;y = sinc(x);DyRe = dsinc(x);DyExp = zeros(3,length(x),length(D));%Derivada numericafor j = 1:length(D)for i = 1:length(x)%a)Derivada hacia adelanteDyExp(1,i,j) = (sinc(x(i)+D(j))-sinc(x(i)))/D(j);

Page 3: Cosa

%b)Derivada centralDyExp(2,i,j) = (sinc(x(i)+D(j))-sinc(x(i)-D(j)))/(2*D(j));%c)Derivada de sexto ordenDyExp(3,i,j) = (256*(sinc(x(i)+D(j))-sinc(x(i)...-D(j)))+40*(sinc(x(i)-2*D(j))-sinc(x(i)+2*D(j)))+...sinc(x(i)+4*D(j))-sinc(x(i)-4*D(j)))/(360*D(j));endend%Errork = zeros(1,length(D));eAbsAd=k;eRelAd=k;eAbsCe=k;eRelCe=k;eAbsSe=k;eRelSe=k;for i = 1:length(D)%Hacia adelanteeAbsAd(i) = mean(abs(DyExp(1,:,i)-DyRe));eRelAd(i) = abs(eAbsAd(i)./mean(DyRe));%CentraleAbsCe(i) = mean(abs(DyExp(2,:,i)-DyRe));eRelCe(i) = abs(eAbsCe(i)./mean(DyRe));%Sexto OrdeneAbsSe(i) = mean(abs(DyExp(3,:,i)-DyRe));eRelSe(i) = abs(eAbsSe(i)./mean(DyRe));end

Graficas �! Ver figura 2.

1/D. D = Desplazamiento entre puntos.0 10 20 30 40 50 60 70 80 90 100

log(E

rror

Rela

tivo

)

10-12

10-10

10-8

10-6

10-4

10-2

100

102Errores Relativos Para diferentes esquemas de derivadas

AdelanteCentralSexto Orden

Figura 2: Error en el calculo de la derivada en funcion del des-plazamiento entre puntos muestra.

El programa varıa la cantidad de puntos utilizados en elintervalo en donde se calcula la derivada con el parametro”D”, el cual es la diferencia entre puntos. Notar que secomienza con un intervalo espaciado en 1, despues .1 yfinalmente .01. El criterio de error que se utilizo, paracada valor de D, fue obtener el valor absoluto del pro-medio de los errores relativos a lo largo de todo el intervalo.

a Utilizando diferencias finitas hacia adelante, obten el

error al calcular la derivada de dicha expresion con-

tra su version analıtica. Grafica dicho error en escala

logarıtmica, para diferentes numeros de puntos utili-

zados en la discretizacion numerica.

El error disminuye, sin embargo este comportamientono es drastico a partir de mas de 10 puntos utilizadospara el calculo.

b Repite el inciso anterior, pero utilizando diferencias

centrales.

La derivada central muestra ser una mejor aproxima-cion de la derivada en cada punto, debido a que esun esquema de segundo orden. Sin embargo, el errormuestra un comportamiento similar al de la derivadahacia adelante a partir del intervalo con alrededor de10 puntos.

c Nuevamente repite el inciso anterior, pero utilizando

la expresion de la primera derivada de sexto orden

(investigar o deducir tal expresion).

Para deducir la aproximacion de sexto orden se uti-lizo un metodo conocido como la “Extrapolacion deRichardson” [1]. El metodo consiste en comenzar conla aproximacion por derivadas centrales y manipularalgebraicamente su expansion en serie de Taylor pa-ra incrementar el orden de aproximacion. La ecuacionfinal de sexto orden obtenida fue:

f 0(x) = [256(f(x+ h)� f(x� h)) + 40(f(x� 2h)

�f(x+ 2h)) + f(x+ 4h)� f(x� 4h)] /24h+O(h6)

En cuanto al calculo numerico, es notoria la diferenciaen la precision de este esquema. En contraste con loslogaritmos de los errores anteriores, el error para laderivada de sexto orden es aproximadamente 9 orde-nes de magnitud mas pequeno. El metodo es superiorsi el objetivo es mayor exactitud en un intervalo conla menor cantidad de puntos.

3. Integrales Numericas

a ¿Que problemas presenta esta integral desde el punto

de vista numerico?

Como se muestra en la siguiente figura, la funcionoscila mucho conforme x ! 0, razon por la cual elcalcular la integral cerca de cero se vuelve compli-cado, pues dependiendo de la discretizacion, puedentomarse valores muy grandes (positivos o negativos)que arrojan resultados muy diferentes para el mismointervalo de integracion.

Figura 3: Grafica de la funcion.

Page 4: Cosa

b ¿Cuales son los mejores resultados que se pueden

obtener utilizando el metodo trapezoidal?

Metodo trapezoidal

function Ints=ej3(eps,n)%eps: limite inferior%n: numero de pasos

%discretizacion de xx=linspace(eps,1,n);

Ints=0;%integracion trapezoidalfor i=1:n-1Int=(x(i+1)-x(i))*((x(i)^(-1)*cos(x(i)^(-1)*log(x(i))))+(x(i+1)^(-1)*cos(x(i+1)^(-1)*log(x(i+1)))))/2;Ints=Ints+Int;endend

Figura 4: Salidas de integracion trapezoidal

El mejor resultado que se obtuvo con el metodotrapezoidal fue 0,32504. Utilizando el intervalo deintegracion de 0,01 a 1. El utilizar mas puntos parala integracion no mejora el resultado. Utilizando unlimite inferior mas pequeno, el comportamiento de laintegral hace que el resultado numerico se aleje muchomas del valor correcto, como se muestra en la figura 4.

c ¿Y cuales son los mejores resultados utilizando una

cuadratura Gaussiana?

Utilizando el metodo de cuadratura Gauss-Legendre, elmejor resultado se obtuvo utilizando el lımite inferior 0,005y 10000 puntos. De manera similar al metodo trapezoi-dal, el comportamiento del metodo se vuelve inestable paralımites inferiores mas pequenos que 0,005, como se observaen la figura 5 al introducir distintos valores de limite infe-rior y de puntos. Enseguida el codigo utilizado, adaptadode [2]Metodo cuadratura Gauss-Legendre

function int=ej3gauss(n,eps)

[x,w]=gauleg(eps,1,n);y=x.^(-1).*cos(x.^(-1).*log(x));int=w*y’;end

function [x,w]=gauleg(x1,x2,n)

m=ceil((n+1)/2);xm=.5*(x2+x1);xl=.5*(x2-x1);eps=1e-14;

%Se hacen las matrices vacıas que se%van a llenar con los pesos y raicesx=zeros(1,n);w=zeros(1,n);

for i=1:m

z=cos(pi*(i-.25)/(n+.5));z1=z+1;

while abs(z-z1)>epsp1=1;p2=0;% algoritmo de Newtonfor j=1:n

p3=p2;p2=p1;p1=((2*j-1)*z*p2-(j-1)*p3)/j;

end

pp=n*(z*p1-p2)/(z*z-1);z1=z;z=z1-p1/pp;

end%aqui se llenan los vetores de pesos y raicesx(i)=xm-xl*z;x(n+1-i)=xm+xl*z;w(i)=2*xl/((1-z*z)*pp*pp);w(n+1-i)=w(i);

endend

Page 5: Cosa

Figura 5: Salidas de integracion por cuadratura gaussiana

4. Sobre investigacion de algunos comandos en

MATLAB

El objetivo del ejercicio era decifrar un mensaje en cifradoCesar. La unica clave fue que el desplazamiento se encuen-tra entre 1 y 20. El codigo programado convierte el mensajecifrado a ASCII, almacenando cada caracter en un elemen-to, y corre un ciclo for restando un numero diferente entre20 y 1 a todo el vector.

%% a) Descifrarclear allclc%Cambiar texto a asciis = double(’mensaje cifrado’);%Prealocar matriz para almacenar datosJ = zeros(20,length(s));for i = 1:20

%Restar i correspondiente a sJ(i,:) = s - i;

end%Convertir codigo de J a caracteresU = char(J);%Despliegar el unico con sentido

U(5,:) %En este caso indica que el%el corrimiento fue de 5.

%% b) cifrar%Desplegar mensaje con cifrado +4char(4+double(U(5,:)));

a Si los valores de las letras son los correspondientes va-

lores en codigo ASCII, ¿que desplazamiento se utilizo

con el siguiente cifrado cesar?

Al imprimir el vector U, la unica entrada que repre-sentaba un mensaje legible era la U(5,:), lo cual indicaque el desplazamiento del cifrado fue de 5. El mensajeobtenido fue:

“No se como puedo servisto por el mundo, pero en mi

opinion, me he comportado como un ni:o que juega al

borde del mar, y que se divierte buscando de cuando

en cuando una piedra mas pulida y una concha mas

bonita de lo normal, mientras que el gran oceano de

la verdad se exponia ante mi completamente descono-

cido.”

b Utilizando un cifrado Cesar con desplazamiento de +4

(con valores en codigo ASCII), ¿como quedarıa cifra-

do el texto anterior?

Con la segunda parte del codigo, el mensaje cifradocon desplazamiento +4 es:

Rs$wi$gsqs$tyihs$wivzmwxs$tsv$ip$qyrhs0$tivs$ir$qm$stmrmsr0$qi$li$gsqtsvxehs$gsqs$yr$rm>s$uyi$nyike$ep$fsvhi$hip$qev0$}$uyi$wi$hmzmivxi$fywgerhs$hi$gyerhs$ir$gyerhs$yre$tmihve$qew$typmhe$}$yre$gsrgle$qew$fsrmxe$hi$ps$rsvqep0$qmirxvew$uyi$ip$kver$sgiers$hi$pe$zivheh$wi$i|tsrme$erxi$qm$gsqtpixeqirxi$hiwgsrsgmhs2

Referencias

[1] Levy D., Introduction to Numerical Analysis. Maryland:CSCAMM, 2010.

[2] Press, H. W. et al., Numerical Recipes in Fortran 77. Cam-bridge University Press, 1992.