Guiones de laboratorio de prácticas de computación
Full text
M´etodos Num´ericos y Computacionales I Ingenier´ıa F´ısica (UPC - ETSETB) Guiones de Laboratorio de Pr´acticas de Computaci´on ` A. Meseguer
Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) M´etodos Num´ericos y Computacionales (2014-2015) A. Meseguer Pr´actica 1: introducci´on al entorno Matlab – I. 0. Comandos b´asicos: Inicio (prompt, ejecutado en terminal): matlab -nodesktop omatlab -nojvm Salir (prompt, ejecutado en terminal): quit Documentacion - ayuda: help ... Variables escalares, memoria, borrado: clear, clear all, clear x y z, who, whos Constantes: i, 1i, pi 1. Representaci´on en coma flotante de R(p. 2 a p. 6 de Quarteroni & Saleri): Formatos: format long/short e/f/g Underflow y overflow (realmin realmax) y variable Inf Epsilon m´aquina: eps 2. Operaciones escalares b´asicas y funciones: Algebraicas: x±y, x*y, x/y, x^y Funciones b´asicas: sqrt, log, log10, exp, sin, cos, tan, asin, acos, atan, abs 3. Vectores, manipulaci´on y operadores vectorizados: Vector fila: x = [x1 x2 x3 ...]=[x1, x2, x3, ...], linspace(a,b,N), [a:inc:b] Vector columna y transposici´on real: x = [x1;x2;x3;...], x.’=[x1, x2, x3, ...] Suma componentes, media y diferencia comp. consecutivas: sum(x), mean(x), diff(x) Ordenaci´on (sorting), m´ax. m´ın.: [sx ii]=sort(x), max(x), min(x) Re-ordenaci´on: flipud([x1;x2;x3]), fliplr([x1 x2 x3]) Indexado: x(j), x(j:k), x(:), x(j:end),x(j:end-k) Extracci´on componentes condicionadas: j = find(c(x)) Ejemplo: x = [0.8234, 0.6948, 0.3171, 0.9502, 0.03444];j = find(abs(x)<.5); j = [3 5] x(j) = [0.3171 0.03444] Vectorizaci´on de funciones: f(x) = [f(x1) f(x2) f(x3) ...]. Ejemplos: x.^p = [x1^p x2^p x3^p ...],1./x = [1/x1 1/x2 1/x3 ...] Operaciones escalar-ay vector-x: a±x = [x1 ±a x2 ±a x3 ±a ...],a*x = [a*x1 a*x2 a*x3 ...] a./x = [a/x1 a/x2 a/x3 ...] Operaciones vector-vector (x=[x1 x2 ...] y = [y1 y2 ...] ): Suma-diferencia: x±y=[x1 ±y1 x2 ±y2 ...] Producto y cociente comp-to-comp: x.*y = y.*x = [x1*y1 x2*y2 ...];x./y=[x1/y1 x2/y2 ...] Producto escalar: x*y.’ = sum(x.*y) = sum(y.*x) Norma eucl´ıdea: norm(x,2)= norm(x) = sqrt(x*x’) Concatenaci´on-actualizaci´on como variable: v = [1 3 7],u = [8 pi],v = [v u] : [1 3 7 8 pi], es decir: v←(v1,...,vn, u1,...,um)
4. Matrices, manipulaci´on, indexado y operaciones: Formato fila-columna: A = [1 pi ; i 4] = [1 , pi ; i , 4] Inicializacion para optimizar velocidad: A = zeros(N) = zeros(N,N),B = zeros(M,N) Matriz identidad (Aij =δij) y matriz de unos (Bij = 1): A = eye(N),B = ones(M,N) Indexaci´on extracci´on columnas: A = [1 2 3 ; 4 5 6],v = A(:,2) −→ v = [2 ; 5] Indexaci´on extracci´on filas: A = [1 2 3 ; 4 5 6],u = A(2,:) −→ u = [4 5 6] Transposici´on real: A = [1 2 3 ; 4 5 6],−→ B = A.’ −→ B = [1 4 ; 2 5 ; 3 6] Suma columnas: A = [1 2 3 ; 4 5 6],sum(A) = [5 7 9],sum(A’) = [6 15] Indexado vectorial: A = [1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ; 13 14 15 16] ,idx = [2 ; 4]: A(idx,2) = [6 ; 14] ,A(2,idx) = [6 8] ,A(2,end) = 8 Diferencia filas consecutivas: A = [1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ; 13 14 15 16] −→ diff(A) = [4 4 4 4 ; 4 4 4 4 ; 4 4 4 4] M´aximo-m´ınimo por columna y absolutos: A = [1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ; 13 14 15 16] : max(A) = [13 14 15 16] ,min(A) = [1 2 3 4],min(min(A)) = 1,max(max(A)) = 16 Ordenaci´on (sorting) por columnas-filas: A = [3 2 1 ; 8 9 9 ; 0 0 5] : sort(A) = [0 0 1 ; 3 2 5 ; 8 9 9],(sort(A’))’ = [1 2 3 ; 8 9 9 ; 0 0 5] Re-ordenaci´on: A=[321;899;005]: flipud(A) =[0 0 5 ; 8 9 9; 3 2 1] , fliplr(A) = [1 2 3 ; 9 9 8 ; 5 0 0] Extracci´on componentes condicionadas: A = [1 2 3 ; 4 5 6 ; 7 8 9] : [i j] = find(A>5),[i j] = [ 3 1 ; 3 2 ; 2 3 ; 3 3] Evaluaci´on de funciones sobre elementos de matriz: A = [1 2 3 ; 4 5 6 ; 7 8 9] f(A) = [f(1) f(2) f(3) ; f(4) ...]. Ejemplos: A.^p = [1 2^p 3^p ; 4^p ...],1./A = [1 1/2 1/3 ; 1/4 ...] Operaciones matriz-vector-escalar: A=[123;456;789],b = [1 ; 0 ; -1], c = 10 : A*b = [-2 ; -2 ; -2], b’*A = [-6 -6 -6] A*c = c*A = [10 20 30 ; 40 50 60 ; 70 80 90] A/c = [0.1 0.2 0.3 ; 0.4 0.5 0.6 ; 0.7 0.8 0.9] Actualizaci´on-ampliaci´on como variable: A = [1 2 3 ; 4 5 6 ; 7 8 9] ,b = [1 ; 0 ; -1]: A = [A b] = [1 2 3 1; 4 5 6 0 ; 7 8 9 -1] , A = [A ; b] = [1 2 3 ; 4 5 6 ; 7 8 9 ; 1 0 -1] Operaciones matriz-matriz (elemento a elemento): A = [1 2 ; 4 5 ] ,B = [1 1 ; -1 -1 ] Suma-diferencia: A + B= [2 3 ; 3 4 ] ... Producto-cociente: A.*B= [1 2 ; -4 -5 ] ... 2
Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) M´etodos Num´ericos y Computacionales (2014-2015) A. Meseguer Pr´actica 2: introducci´on al entorno Matlab – II. 1. Archivos .m, funciones y par´ametros: Loops: for ii=1:N,for ii=1:L:N,end Condicionales, operadores flujo: if a == b, c = ..., end,==, <, >, <=, >=, ~=. Loops condicionados: while jj ~= N, for ii=1:N ... end, end Funciones externas: y = function(x,a,b,c,...),[y z] = function(x,a,b,c,...). Funciones con argumentos funcionales: y = function(x,a,b,...,@fun,param); Tiempo de c´alculo: tic ... toc,cputime Opcional: utilidades debugging: dbstop, etc. 2. Representaci´on gr´afica: plot(x,y),loglog(x,y),semilogy(x,y),semilogx(x,y). Ejemplos de representaci´on de funciones de crecimiento exponencial. Identificaci´on de leyes y∼1/xp,y∼log x, , y∼ex. Evaluar funciones externas y representarlas. Exportar gr´aficos: print -deps figurename,print -depsc figurename 3. Input-Output: Guardar-recuperar datos: load file.dat,save filename A B C ..., archivos .mat. Escritura/formato: fprintf,sprintf,int2str,num2str
Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) M´etodos Num´ericos y Computacionales (2011-2012) Pr´actica 3: resoluci´on num´erica de ecuaciones trascendentes (aplicaciones f´ısicas). 1. Problema: nos encontramos en el valle de un mont´ıculo de unos 50 metros de altura aproximadamente y queremos lanzar una pelota al otro lado de la cima del mismo. El perfil del mont´ıculo viene dado por la ecuaci´on y(x) = ax2e−bx. Tanto ycomo xvienen dados en metros, siendo a= 0.15 m−1yb= 0.04 m−1. Lanzamos la pelota desde (x0, y0) = (0,0) m con una velocidad inicial v0= 37 ms−1y un ´angulo inicial de α0= 67.5o(ver figura). x(m) y(m) 0 50 100 150 200 0 20 40 60 (a) Plantear la ecuaci´on F(x, v0, α0) = 0 que determina la intersecci´on de la par´abola con el perfil del mont´ıculo. Editar una funci´on .m de matlab para F(x, v0, α0) y representarla para los par´ametros dados para x∈[0,200]. (b) Desde un c´odigo principal .m de matlab, determinar el punto de impacto con el m´etodo de la secante para los valores dados. (c) Con el mismo ´angulo de tiro, determinar la velocidad m´ınima necesaria para que la pelota pase a la derecha de la cima del mont´ıculo. Comentar qu´e tipo de condici´on hay que imponer y qu´e consecuencias tiene en la convergencia del m´etodo de la secante.
Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) M´etodos Num´ericos y Computacionales (2011-2012) Pr´actica 4: resoluci´on num´erica de ecuaciones trascendentes (evaluaci´on). 1. Problema: un pavimento muy irregular tiene el perfil vertical que se muestra en la figura de la izquierda. El perfil representado viene dado por la expresi´on f(x) = asin2(bx) e−cx2, en donde xeyvienen dados en metros, a= 5 m, b= 1 m−1yc= 0.06 m−2. x(m) y(m) 0 2 4 6 8 10 0 1 2 3 4 5 θ ~g y x µ f(x) (a) Suponiendo que la superficie no tiene rozamiento, determinar los puntos de equilibrio en x∈[0,8] en los cuales un cuerpo en reposo se mantendr´ıa siempre en esa posici´on sin deslizar por la superficie. Para ello plantear la ecuaci´on F(x, a, b, c) = 0 que determina la condici´on de equilibrio. Representar Fen el intervalo x∈[0,8] y determinar los puntos de equilibrio con un m´etodo adecuado. (b) Supongamos que ahora el pavimento tiene un rozamiento est´atico uniforme µconocido. Determinar las regiones (subintervalos de x∈[0,8]) en las cuales podr´ıamos dejar un bloque en reposo sin que posteriormente este se deslizase por la acci´on de la gravedad (tomar µ= 0.3 y g = 10 ms−2).
Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) M´etodos Num´ericos y Computacionales (2015-2016) A. Meseguer Pr´actica 5: interpolaci´on de Lagrange equiespaciada. 1. Problema: dado un conjunto de puntos equiespaciados {x0=−1, x1, . . . , xn−1, xn= 1}, con: xj=−1 + 2j n,(j= 0,1, . . . , n). Se pide: (a) determinar la matriz de polinomios cardinales de Lagrange P∈Mm+1,n+1(R)1evaluada en un conjunto de puntos {z0, . . . , zm}denso (m≈600, por ejemplo): P= `0(z0)`1(z0). . . `n(z0) `0(z1)`1(z1). . . `n(z1) `0(z2)`1(z2). . . `n(z2) . . .. . ..... . . `0(zm)`1(zm). . . `n(zm) ,con `i(z) = n Y j=0 (j6=i) z−xj xi−xj , i = 0,1, . . . , n. Nota: procurar no hacer uso de bucles anidados para calcular los productos. Un s´olo bucle para los denominadores y otro para los numeradores deber´ıa ser suficiente. (b) para n= 3,6 y 9, representar los `i(zk) en el conjunto de puntos zkpara poder visualizar el comportamiento de dichos polinomios entre nodos y comprobar que cumplen `i(xj) = δij . (c) para n= 8,16,24,32, calcular la funci´on de Lebesgue: λn(x) = n X j=0 |`j(x)| sobre el mismo conjunto de puntos zky comprobar su crecimiento en los extremos del intervalo. Utilizar una escala adecuada en los ejes. (d) interpolar la funci´on f(x)=exen [−1,1] para n= 4,6,8,10,...,60 mediante un bucle. Para cada n, evaluar la interpolaci´on en el conjunto de nodos zkhaciendo el producto matriz-vector Πnf(zi) = PjPij f(xi), evaluar el error m´aximo εn= m´ax 0≤k≤m|Πnf(zk)−ezk| y almacenarlo en un vector. Representar en una gr´afica semilogar´ıtmica εnfrente a ny comprobar la ley εn∼2n/n log(n). Comprobar que la extrapolaci´on de εnpara n→0 intercepta el eje ycerca de εmach ∼10−16. 1Mmn(R) es el espacio de matrices reales de mfilas y ncolumnas
Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) M´etodos Num´ericos y Computacionales (2014-2015) Pr´actica 6: interpolaci´on de Chebychev. 1. Problema: interpolar la funci´on: f(x) = tanh {20 sin(12x)}+2 100e3xsin(300x), en el intervalo [0,1]. Utilizar el conjunto de nodos de Chebychev: zj= cos jπ n, j = 0,1, . . . , n. El intervalo es el [0,1], con lo que habr´a que realizar un cambio de variable adecuado, es decir, los nodos sobre los que se realiza la interpolaci´on son: xj=1 2(1 + zj). Se puede utilizar la forma de Lagrange u, opcionalmente, la baric´entrica (la f´ormula no cambia a pesar del cambio de variable): Πnf(x) = n X j=0 0(−1)jfj(x−xj)−1 n X j=0 0(−1)j(x−xj)−1 ,∀x6=xj. Recu´erdese que el s´ımbolo 0en los sumatorios indica que el primer y ´ultimo t´erminos de las sumas deben dividirse por 2. (a) Representar la funci´on en un conjunto de 2000 puntos para visualizarla primero. (b) Utilizar n= 102,103, . . ., etc., nodos de interpolaci´on e ir visualizando el error local: ε(x) = |Πnf(x)−f(x)|, cometido a medida que se aumenta el n´umero de puntos. Para ello utilizar una gr´afica semilogar´ıtmica. (c) Utilizar el n´umero de nodos necesario para que el ε(x)<10−12 en todo el intervalo [0,1].
Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) M´etodos Num´ericos y Computacionales (2016-2017) A. Meseguer Pr´actica 7: derivaci´on num´erica (I). Introducci´on: la construcci´on de matrices de diferenciaci´on con Matlab es relativamente sencilla si dichas matrices tienen una estructura simple. En particular, las matrices de diferenciaci´on resultantes de aplicar f´ormulas locales de orden bajo suelen ser de tipo Toeplitz, es decir con sus diagonales principales adoptando un mismo valor constante. 1. Comando toeplitz:este comando construye una matriz de Toeplitz partiendo de las componentes de dos vectores u= (u1, u2, . . . , un) y v= (v1, v2, . . . , vn) dados (con v1=u1), ubic´andolas diagonalmente de la siguiente forma: u1v2v3v4· · · vn−1vn u2u1v2v3· · · vn−2vn−1 u3u2u1v2· · · vn−3vn−2 u4u3u2u1· · · vn−4vn−3 . . .. . .. . .. . .. . .. . . un−1un−2un−3un−4· · · u1v2 unun−1un−2un−3· · · u2u1 Ejemplo: introducir en la l´ınea de comandos: toeplitz([1 3 5 7],[1 2 3 4]). 2. Utilizando el comando anterior, construir las matrices de diferenciaci´on fd ycd vistas en clase de teor´ıa para n= 5,20,40,80: −1 1 0 0 . . . 0 0 0−110. . . 0 0 0 0 −1 1 . . . 0 0 . . .. . .. . .. . .. . .. . . 0 0 0 0 . . . 1 0 0 0 0 0 . . . −1 1 y −1/2 0 1/2 0 . . . 0 0 0 0 0−1/201/2. . . 0 0 0 0 0 0 −1/2 0 . . . 0 0 0 0 . . .. . .. . .. . .. . .. . .. . .. . . 0 0 0 0 . . . −1/201/2 0 0 0 0 0 . . . 0−1/201/2 , asociadas a las f´ormulas uFD i=h−1(fi+1 −fi) y uCD i= (2h)−1(fi+1 −fi−1). Al aplicar el comando toeplitz, recordar que hay que eliminar algunas filas y que no podemos obtener las derivadas en todos los puntos originales. No olvidar los factores hen las matrices. Visualizar las matrices para n= 5 en la l´ınea de comandos para comprobar que tienen el n´umero de filas y columnas adecuado. 3. Considerar la funci´on f(x) = sin(πx) y su derivada f0(x) = πcos(πx) en el intervalo [0,2]. Evaluar f(x) en un conjunto de nodos equiespaciado xj= 2j/n, (j= 0,1, . . . , n), con n= 5,20,40,80. Aplicar las matrices anteriores sobre el vector de valores (f0, f1, . . . , fn) para obtener una aproximaci´on de f0(x). Representar la derivada exacta sobre los nodos f0(xj) = πcos(πxj) y, en la misma gr´afica, las aproximaciones obtenidas para ver la discrepancia. En otra gr´afica en escala semilogar´ıtmica representar |f0(xi)−uFD i|y|f0(xi)−uCD i|. 4. (opcional): Para n= 100,200,300,400,...,2000, calcular εn= m´ax |f0(xi)−ui|para fd ycd. En una gr´afica loglog, representar εnfrente al hutilizado en cada caso e interpretar la pendiente.