Simulación de ondas con diferencias finitas (Matlab)Tr

SIMULACIÓN DE ONDAS EN DOMINIOS INFINITOS TRUNCADOS APROXIMACIÓN MEDIANTE DIFERENCIAS FINITAS EN MATLAB ®

Asignatura: Asignatura: Métodos Numéricos, Modelización y Sistemas Dinámicos Profesora: Profesora: Ana Carpio Rodríguez Autor: Autor: Alfonso de la Fuente Ruiz Fecha: Fecha: Septiembre de 2009

Trusted by over 1 million members

Try Scribd FREE for 30 days to access over 125 million titles without ads or interruptions! Start Free Trial Cancel Anytime.









INTRODUCCIÓN Este trabajo se presenta para optar a la superación parcial de la asignatura “ Métodos

Numéricos, Modelización y Sistemas Dinámicos ” (en adelante: “MNMSD”) del Master Oficial de Postgrado en Ingeniería Matemática de la Universidad Complutense de Madrid. Madrid. Su contenido constituye una extensión a dos dimensiones más el tiempo, discretizando por el “Método numérico de las Líneas” (Numeric MOL: ‘Method Of Lines’ *6 pág. 184+) de una simulación de ondas en dominios infinitos truncados con forma cuadrada. Puede interpretarse físicamente como la oscilación de una membrana. membrana.

Los resultados espaciotemporales se aproximan por un esquema explícito de segundo orden basado en diferencias finitas (FDTD ( FDTD)) para dos juegos de condiciones de contorno diferentes. En tercer lugar se aproxima aproxima la solución en un recinto mayor también cuadrado, con el triple de lado. Los programas se han implementado sobre MATLAB R2006a (v7.2.0.232) y los vídeos se









ENUNCIADO Y NOTAS









NOTAS DE LA VERSIÓN: Se utiliza como base el algoritmo ya implementado para la ecuación de ondas en una dimensión (ming7c.pdf )  +1 Se desarrolla la ecuación por diferencias finitas hasta despejar el término t érmino  , .

Para simular dos juegos de condiciones de contorno, seguidamente se despejan también las  +1  +1  +1  +1 condiciones de contorno no nulas:  , , 0, ,  , ,  ,0 . A continuación se inicializan las condiciones de partida. Finalmente se introduce una rutina que permite exportar el vídeo generado a un fichero en formato MPEG. Para ello se emplea el ‘plug-in’ “MPGwrite 6.0 ” de MATLAB, que debe descomprimirse en la misma carpeta que contenga el fichero principal. Dicho plugin puede descargarse gratuitamente de la siguiente dirección: http://www.mathworks.com/matlabcentral/fileexchange/309 A partir del código fuente desarrollado para la primera aproximación, se obtiene el de las simulaciones segunda y tercera, mediante una serie de extensiones y modificaciones.



























CÓDIGO FUENTE PRIMERA APROXIMACIÓN, CONDICIONES DE CONTORNO NULAS %UCM-MIM 2009 MNMSD:ANA CARPIO %AUTOR:ALFONSO DE LA FUENTE RUIZ %function [U_save,x,y,t_save]= ondas2D(c,t_f,M,tsteps); %Resuelve utt - c2 (uxx +uyy) = 0 mediante un esquema explicito de segundo orden %Primera aproximación. Contorno nulo. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Parámetros y datos de entrada: Nframes=180; % % Nº de cuadros de vídeo totales Nframes=input('Introduzca el Nº de fotogramas (Def.=180): '); NFramesPorIt=3; % Nº cuadros de vídeo por cada iteración dibujada c=1; % coeficiente c de la ecuacion (Velocidad de la onda) c=input('Introduzca el valor de c, velocidad de onda (Def.=1): '); t_f=6; % Tiempo final hasta el que calculamos la solución t_f=input('Introduzca el instante final de tiempo (Def.=6): '); M=12; %% N=length(x)=lenght(y) % Nº de partes para discretizar cada intervalo [-L,L] M=input('Introduzca el Nº de divisiones de la rejilla (Def.=12): '); ResX=640; ResY=240; ResX=input('Introduzca la resolución de vídeo en el eje X (Def.=640): '); ResY=input('Introduzca la resolución de vídeo en el eje Y (Def.=240): '); 'Inicializando primera aproximación' L=1; % Semieje de intervalo cuadrado [-1,1] h=2*L/M; % dx=dy Longitud de paso espacial tsteps=360; %% Nº de partes iguales en que se descompone [0,t_f] k=c/(4*h^2); %% Paso espacial dt=t_f/tsteps; % Paso temporal tiempos=0:dt:t_f; % Mallado espacio-temporal % Recinto pequeño: x=-L:h:L; y=-L:h:L; [X,Y]=meshgrid(x,y); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%Condiciones iniciales: U0=exp(-10*(X.^2+Y.^2)); %% u0 % Posición inicial U1=zeros(size(X)); % Velocidad inicial %Condiciones iniciales, sin contornos: U0_red=U0(2:M,2:M); % Posición sin contorno U1_red=U1(2:M,2:M); % Velocidad sin contorno %Versiones aplanadas:









t_save(1)=0; % Cuadro inicial de animación % OPERADOR D: Usamos el comando diag(vector,k) para crear la matriz tridiagonal % El comando kron se usa para poner unos en los pares de índices: Tridiag = diag(-4*ones((M-1)^2,1))+... kron(eye(M-1),diag(ones(M-2,1),1))+... kron(eye(M-1),diag(ones(M-2,1),-1))+... kron(diag(ones(M-2,1),1),eye(M-1))+... kron(diag(ones(M-2,1),-1),eye(M-1)); D=(1/h^2)*Tridiag; % Crear el operador discretizado % u_ant=u0; % Inicializar paso anterior %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u = u0 + dt*u1 + 0.5*(dt^2)*(c^2)*D*u0; % Cálculo del paso 2 %u(i,l,j+1)=2*u(i,l,j)-u(i,l,j-1)+(c^2*k^2/h^2)*(u(i+1,l,j)+u(i-1,l,j)+u(i,l+1,j)+u(i,l1,j)-4*u(i,l,j)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 'Calculando. Espere, por favor...' % ITERAR para el resto de pasos temporales: for n=1:tsteps temp=u; u = 2*u - u_ant + (dt^2)*(c^2)*D*u; u_ant=temp; if mod(n,marca)==0 % Guardamos los valores de u sólo para algunos tiempos: indice = 1 + n/marca; U(2:M,2:M)=reshape(u,M-1,M-1); U_save(:,:,indice)=U; t_save(indice)=tiempos(n); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 'Cálculo completo. Creando visualización...' % Representación: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Parámetros de visualización: camlight(camlight('right'),'right'); %view(0,90); % Angulo de visualización shading interp; lighting phong; colormap copper; axis off; zlim([-1 1]); colorbar; set(gca,'NextPlot','replacechildren'); set(gcf, 'Renderer', 'ZBuffer'); set(gcf,'MenuBar','none'); % Quitar la barra de menú set(gcf,'Position',[20 40 ResX ResY]);% Situar la figura a 20*40 píxeles con tamaño ResX*ResY: P1=moviein(NFramesPorIt*indice); % Crear película 1 [X,Y]=meshgrid(x,y); for i=1:indice surf(U_save(:,:,i)) set(gcf,'Color', [0 0 0], 'Number', 'off', 'Name', sprintf('Onda2D: N=%i, T=%.2f msec',i-1,t_save(i)*1e3)); for nf=1:NFramesPorIt % Guardar 'NFramesPorIt' cuadros en el vídeo: P1((i-1)*NFramesPorIt+nf)=getframe(gcf);









SEGUNDA APROXIMACIÓN, CONDICIONES DE CONTORNO ESPECIALES

Nota: para evitar colisiones entre asignaciones a contornos horizontales y verticales, las esquinas (2*2) del contorno se igualarán a 0 en cada iteración, en las tres celdillas accesibles (marcadas en rojo sobre la figura).

%UCM-MIM 2009 MNMSD:ANA CARPIO %AUTOR:ALFONSO DE LA FUENTE RUIZ %function [U_save,x,y,t_save]= ondas2D(c,t_f,M,tsteps); %Resuelve utt - c2 (uxx +uyy) = 0 mediante un esquema explicito de segundo orden %Segunda aproximación. Contorno especial. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Parámetros y datos de entrada: Nframes=120; % % Nº de cuadros de vídeo totales Nframes=input('Introduzca el Nº de fotogramas (Def.=120): '); NFramesPorIt=3; % Nº cuadros de vídeo por cada iteración dibujada c=1; % coeficiente c de la ecuacion (Velocidad de la onda) c=input('Introduzca el valor de c, velocidad de onda (Def.=1): '); t_f=6; % Tiempo final hasta el que calculamos la solución t_f=input('Introduzca el instante final de tiempo (Def.=6): '); M=11; %% N=length(x)=lenght(y) % Nº de partes para discretizar cada intervalo [-L,L] M=input('Introduzca el Nº de divisiones de la rejilla (Def.=11): '); ResX=640; ResY=240; ResX=input('Introduzca la resolución de vídeo en el eje X (Def.=640): '); ResY=input('Introduzca la resolución de vídeo en el eje Y (Def.=240): '); 'Inicializando segunda aproximación' L=1; % Semieje de intervalo cuadrado [-1,1] h=2*L/M; % dx=dy Longitud de paso espacial tsteps=360; %% Nº de partes iguales en que se descompone [0,t_f] k=c/(4*h^2); %% Paso espacial dt=t_f/tsteps; % Paso temporal tiempos=0:dt:t_f; % Mallado espacio-temporal % Recinto pequeño: x=-L:h:L; y=-L:h:L; [X,Y]=meshgrid(x,y); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%Condiciones iniciales:









% Condiciones de frontera: U0_izq=zeros(1,M+1,Nframes); % Lado izquierdo U0_dch=zeros(1,M+1,Nframes); % Lado derecho U0_inf=zeros(M+1,1,Nframes); % Lado inferior U0_sup=zeros(M+1,1,Nframes); % Lado superior % Inicialización de frontera: U0int_izq=U0(2,2:M);% Lado izquierdo interior U0int_dch=U0(M,2:M);% Lado derecho interior U0int_inf=U0(2:M,2);% Lado inferior interior U0int_sup=U0(2:M,M);% Lado superior interior for i=3:abs(M+1-2) U1_izq(i)=U0_izq(i)+c*k/h*(U0int_izq(i-1)-U0_izq(i)); end for i=3:abs(M+1-2) U1_dch(i)=U0_dch(i)-c*k/h*(U0_dch(i)-U0int_dch(i-1)); end for i=3:abs(M+1-2) U1_sup(i)=U0_sup(i)-c*k/h*(U0_sup(i)-U0int_sup(i-1)); end for i=3:abs(M+1-2) U1_inf(i)=U0_inf(i)+c*k/h*(U0int_inf(i-1)-U0_inf(i)); end %4 esquinas 2*2-1 nulas (decisión de diseño) U1_izq(2)=0; U1_izq(1)=0; U1_dch(2)=0; U1_dch(1)=0; U1_sup(1)=0; U1_sup(2)=0; U1_inf(1)=0; U1_inf(2)=0; U1_izq(M)=0; U1_izq(M+1)=0; U1_dch(M)=0; U1_dch(M+1)=0; U1_sup(M)=0; U1_sup(M+1)=0; U1_inf(M)=0; U1_inf(M+1)=0; % 4 redund. %Inicialización U_save(1,:,:)=U0_izq; U_save(M+1,:,:)=U0_dch; U_save(:,1,:)=U0_inf; U_save(:,M+1,:)=U0_sup; % U_save(:,:,1)=U0; % Posición inicial t_save(1)=0; % Cuadro inicial de animación % OPERADOR D: Usamos el comando diag(vector,k) para crear la matriz tridiagonal % El comando kron se usa para poner unos en los pares de índices: Tridiag = diag(-4*ones((M-1)^2,1))+... kron(eye(M-1),diag(ones(M-2,1),1))+... kron(eye(M-1),diag(ones(M-2,1),-1))+... kron(diag(ones(M-2,1),1),eye(M-1))+... kron(diag(ones(M-2,1),-1),eye(M-1)); D=(1/h^2)*Tridiag; % Crear el operador discretizado % u_ant=u0; % Inicializar paso anterior %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%









U_save(:,:,indice)=U; t_save(indice)=tiempos(n); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 'Cálculo completo. Creando visualización...' % Representación: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Parámetros de visualización: camlight(camlight('right'),'right'); %view(0,90); % Angulo de visualización shading interp; lighting phong; colormap copper; axis off; zlim([-1 1]); colorbar; set(gca,'NextPlot','replacechildren'); set(gcf, 'Renderer', 'ZBuffer'); set(gcf,'MenuBar','none'); % Quitar la barra de menú set(gcf,'Position',[20 40 ResX ResY]);% Situar la figura a 20*40 píxeles con tamaño ResX*ResY: P2=moviein(NFramesPorIt*indice); % Crear película 1 [X,Y]=meshgrid(x,y); for i=1:indice surf(U_save(:,:,i)) set(gcf,'Color', [0 0 0], 'Number', 'off', 'Name', sprintf('Onda2D: N=%i, T=%.2f msec',i-1,t_save(i)*1e3)); for nf=1:NFramesPorIt % Guardar 'NFramesPorIt' cuadros en el vídeo: P2((i-1)*NFramesPorIt+nf)=getframe(gcf); end end %pause(.1); end 'Animación completa. Exportando vídeo...' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% save onda2D_P2.mat P2; % Guardar la animación MATLAB a un archivo mpgwrite(P2,jet,'onda2D_P2.mpg'); % Convertir la animación en vídeo MPEG 'Programa terminado. Que tenga un buen día :)' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%











TERCERA APROXIMACIÓN, RECINTO MAYOR %UCM-MIM 2009 MNMSD:ANA CARPIO %AUTOR:ALFONSO DE LA FUENTE RUIZ %function [U_save,x,y,t_save]= ondas2D(c,t_f,M,tsteps); %Resuelve utt - c2 (uxx +uyy) = 0 mediante un esquema explicito de segundo orden %Tercera aproximación. Contorno nulo. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Parámetros y datos de entrada: Nframes=180; % % Nº de cuadros de vídeo totales Nframes=input('Introduzca el Nº de fotogramas (Def.=180): '); NFramesPorIt=3; % Nº cuadros de vídeo por cada iteración dibujada c=1; % coeficiente c de la ecuacion (Velocidad de la onda) c=input('Introduzca el valor de c, velocidad de onda (Def.=1): '); t_f=4; % Tiempo final hasta el que calculamos la solución t_f=input('Introduzca el instante final de tiempo (Def.=4): '); M=36; %% N=length(x)=lenght(y) % Nº de partes para discretizar cada intervalo [-3*L,3*L] M=input('Introduzca el Nº de divisiones de la rejilla (Def.=36): '); ResX=640; ResY=240; ResX=input('Introduzca la resolución de vídeo en el eje X (Def.=640): '); ResY=input('Introduzca la resolución de vídeo en el eje Y (Def.=240): '); 'Inicializando primera aproximación' L=3; % Semieje de intervalo cuadrado [-3,3] h=2*L/M; % dx=dy Longitud de paso espacial tsteps=360; %% Nº de partes iguales en que se descompone [0,t_f] k=c/(4*h^2); %% Paso espacial dt=t_f/tsteps; % Paso temporal tiempos=0:dt:t_f; % Mallado espacio-temporal % Recinto pequeño: x=-L:h:L; y=-L:h:L; [X,Y]=meshgrid(x,y); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%Condiciones iniciales: U0=exp(-10*(X.^2+Y.^2)); %% u0 % Posición inicial









kron(diag(ones(M-2,1),1),eye(M-1))+... kron(diag(ones(M-2,1),-1),eye(M-1)); D=(1/h^2)*Tridiag; % Crear el operador discretizado % u_ant=u0; % Inicializar paso anterior %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u = u0 + dt*u1 + 0.5*(dt^2)*(c^2)*D*u0; % Cálculo del paso 2 %u(i,l,j+1)=2*u(i,l,j)-u(i,l,j-1)+(c^2*k^2/h^2)*(u(i+1,l,j)+u(i-1,l,j)+u(i,l+1,j)+u(i,l1,j)-4*u(i,l,j)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 'Calculando. Espere, por favor...' % ITERAR para el resto de pasos temporales: for n=1:tsteps temp=u; u = 2*u - u_ant + (dt^2)*(c^2)*D*u; u_ant=temp; if mod(n,marca)==0 % Guardamos los valores de u sólo para algunos tiempos: indice = 1 + n/marca; U(2:M,2:M)=reshape(u,M-1,M-1); U_save(:,:,indice)=U; t_save(indice)=tiempos(n); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 'Cálculo completo. Creando visualización...' % Representación: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Parámetros de visualización: camlight(camlight('right'),'right'); %view(0,90); % Angulo de visualización shading interp; lighting phong; colormap copper; axis off; zlim([-1 1]); colorbar; set(gca,'NextPlot','replacechildren'); set(gcf, 'Renderer', 'ZBuffer'); set(gcf,'MenuBar','none'); % Quitar la barra de menú set(gcf,'Position',[20 40 ResX ResY]);% Situar la figura a 20*40 píxeles con tamaño ResX*ResY: P3=moviein(NFramesPorIt*indice); % Crear película 3 [X,Y]=meshgrid(x,y);









RESULTADOS Se presentan a continuación los resultados de ejecución de los programas. Se permite elegir los parámetros por interpelación interactiva vía consola. Las animaciones se exportan a ficheros en formato estándar MPEG y vídeo MATLAB en el mismo directorio. Para exportar los vídeos es necesario instalar el ‘plug-in’ “MPGwrite 6.0 ” ya mencionado. El vídeo resultante para la primera aproximación, aproximación, con condiciones de contorno nulas, puede reproducirse en la siguiente dirección del portal de vídeo en internet YouTube: http://www.youtube.com/watch?v=J6UkKmmAq0









Se muestra el primer fotograma de la ejecución de la segunda aproximación, aproximación, inicializada con condiciones de contorno especiales. Como ya se ha indicado, se tomó la decisión de diseño de anular las celdillas de las cuatro esquinas, a falta de indicaciones a este respecto en el enunciado. Junto con el código fuente se incluyen también algunos cálculos intermedios realizados a mano para deducir las condiciones de contorno y traducirlas a código MATLAB.

El vídeo resultante para la tercera aproximación aproximación,, con condiciones de contorno nulas sobre un recinto con el triple de lado, puede reproducirse en la siguiente dirección del portal de vídeo en internet YouTube: http://www.youtube.com/watch?v=o4s_thweckM









El código fuente para esta aproximación se ha parametrizado para observar el comienzo del rebote de retorno. retorno.

Como posible mejora se sugiere combinar la segunda aproximación con la tercera, ejecutando ambas en paralelo, de modo que las condiciones de contorno para el recinto menor se actualicen, en cada iteración, con los valores de función que en esas mismas celdillas almacena el recinto mayor.











BIBLIOGRAFÍA Esta bibliografía incluye una serie de libros, artículos y otras referencias que se han consultado para la realización del presente trabajo.

[1]: “Apuntes y guiones de prácticas de la asignatura MNMSD“ del Máster Oficial de Postgrado en Ingeniería Matemática adaptado al Espacio Europeo de Educación Superior (EEES/EHEA), Universidad Complutense de Madrid, curso 2006-2007. *2+: “Numerical Methods for Problems with Moving Fronts ” by Bruce A. Finlayson, University

of Washington. Ravenna Park Publishing Inc. (1992) *3+: “Numerical Analysis. Mathematics of Scientific Computing” by David Kincaid and Ward

Cheney, University of Texas at Austin. Brooks/Cole (1991) *4+: “A compendium of partial differential equation models: Method Of Lines analysis with

Matlab” by William E. Schiesser and Graham W. Griffiths. Cambridge University Press (2009)

Simulación de ondas con diferencias finitas (Matlab)Tr

Recommend Documents