Source: https://es.scribd.com/doc/26435699/Simulacion-de-ondas-con-diferencias-finitas-Matlab-Tr
Timestamp: 2016-05-30 13:09:56+00:00

Document:
SubirSign inJoinBooksAudiobooksComicsSheet MusicEditors' Picks BooksHand-picked favorites from our editorsEditors' Picks AudiobooksHand-picked favorites from our editorsEditors' Picks ComicsHand-picked favorites from our editorsEditors' Picks Sheet MusicHand-picked favorites from our editorsTop BooksWhat's trending, bestsellers, award-winners & moreTop AudiobooksWhat's trending, bestsellers, award-winners & moreTop ComicsWhat's trending, bestsellers, award-winners & moreTop Sheet MusicWhat's trending, bestsellers, award-winners & moreCategoriesArts & IdeasBiography & MemoirBusiness & LeadershipChildren'sComputers & TechnologyCooking & FoodCrafts & HobbiesFantasyFiction & LiteratureHappiness & Self-HelpHealth & WellnessHistoryHome & GardenHumorLGBTMystery, Thriller & CrimePolitics & EconomyReferenceReligionRomanceScience & NatureScience FictionSociety & CultureSports & AdventureTravelYoung AdultCategoriesArts & IdeasBiography & MemoirBusiness & LeadershipChildren'sComputers & TechnologyCooking & FoodFantasyFiction & LiteratureHappiness & Self-HelpHealth & WellnessHistoryHome & GardenHumorLGBTMystery, Thriller & CrimePolitics & EconomyReferenceReligionRomanceScience & NatureScience FictionSociety & CultureSports & AdventureTravelYoung AdultCategoriesAdaptationsChildren’sCrime & MysteryFictionHumorMangaNonfictionRomanceSciFi, Fantasy & HorrorSuperheroesYoung AdultPublishersArcanaArchie ComicsBOOM! StudiosDynamiteIDW PublishingKingstone ComicsMarvel ComicsSpace Goat ProductionsTop Cow ComicsTop Shelf ProductionsValiant Comics ZenescopeDifficultyBeginnerIntermediateAdvancedMixedInstrumentBrassDrums & PercussionGuitar, Bass, and FrettedPianoStringsVocalWoodwindsGenreClassicalCountryFolkJazz & BluesMovies & MusicalsPop & RockReligious & HolidayStandardsWelcome to Scribd! Start your free trial and access books, documents and more.Find out moreSIMULACIÓN DE ONDAS EN DOMINIOS INFINITOS TRUNCADOS APROXIMACIÓN MEDIANTE DIFERENCIAS FINITAS EN MATLAB ®Asignatura: Métodos Numéricos, Modelización y Sistemas Dinámicos Profesora: Ana Carpio Rodríguez Autor: Alfonso de la Fuente Ruiz Fecha: Septiembre de 2009
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. 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.
Los resultados espaciotemporales se aproximan por un esquema explícito de segundo orden basado en diferencias finitas (FDTD) para dos juegos de condiciones de contorno diferentes. En tercer lugar se 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 exportaron en un Pentium IV sobre una máquina virtual (VirtualBox) que ejecutaba Windows XP SP2 con 384 MB RAM. Se incluye una bibliografía de referencias consultadas. MIM-UCM Página 2
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 𝑢𝑖,𝑙 .
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.
A continuación se presentan algunos cálculos intermedios y diagramas realizados:
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: u0=reshape(U0_red,(M-1)^2,1); u1=reshape(U1_red,(M-1)^2,1); %Utilizamos la matriz U para trabajar cómodamente con el primer índice %para la coordenada ’x’ y el segundo para la coordenada ’y’ U=zeros(M+1,M+1); U_save=zeros(M+1,M+1,Nframes); t_save=zeros(Nframes,1); % Vector de cuadros de animación cuentaframe=1; %% marca=floor(tsteps/(Nframes-1)); %% Para saber qué iteración almacenar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Condiciones de frontera: U_save(1,:,:)=zeros(1,M+1,Nframes); % Lado izquierdo U_save(M+1,:,:)=zeros(1,M+1,Nframes); % Lado derecho U_save(:,1,:)=zeros(M+1,1,Nframes); % Lado inferior U_save(:,M+1,:)=zeros(M+1,1,Nframes); % Lado superior 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 = 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% '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); end %pause(.1); end 'Animación completa. Exportando vídeo...' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% save onda2D_P1.mat P1; % Guardar la animación MATLAB a un archivo mpgwrite(P1,jet,'onda2D_P1.mpg'); % Convertir la animación en vídeo MPEG 'Programa completo. Que tenga un buen día :)' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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: 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: u0=reshape(U0_red,(M-1)^2,1); u1=reshape(U1_red,(M-1)^2,1); %Utilizamos la matriz U para trabajar cómodamente con el primer índice %para la coordenada ’x’ y el segundo para la coordenada ’y’ U=zeros(M+1,M+1); U_save=zeros(M+1,M+1,Nframes); t_save=zeros(Nframes,1); % Vector de cuadros de animación cuentaframe=1; %% marca=floor(tsteps/(Nframes-1)); %% Para saber qué iteración almacenar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 = 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)); % Inicialización de nueva frontera ++ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% '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; %calcular nuevos contornos ++ 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); %asignar nuevos contornos ++
U_save(:,:,indice)=U; t_save(indice)=tiempos(n); 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 %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 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: u0=reshape(U0_red,(M-1)^2,1); u1=reshape(U1_red,(M-1)^2,1); %Utilizamos la matriz U para trabajar cómodamente con el primer índice %para la coordenada ’x’ y el segundo para la coordenada ’y’ U=zeros(M+1,M+1); U_save=zeros(M+1,M+1,Nframes); t_save=zeros(Nframes,1); % Vector de cuadros de animación cuentaframe=1; %% marca=floor(tsteps/(Nframes-1)); %% Para saber qué iteración almacenar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Condiciones de frontera: U_save(1,:,:)=zeros(1,M+1,Nframes); % Lado izquierdo U_save(M+1,:,:)=zeros(1,M+1,Nframes); % Lado derecho U_save(:,1,:)=zeros(M+1,1,Nframes); % Lado inferior U_save(:,M+1,:)=zeros(M+1,1,Nframes); % Lado superior 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 = 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% '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); 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: P3((i-1)*NFramesPorIt+nf)=getframe(gcf); end %pause(.1); end 'Animación completa. Exportando vídeo...' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% save onda2D_P3.mat P3; % Guardar la animación MATLAB a un archivo mpgwrite(P3,jet,'onda2D_P3.mpg'); % Convertir la animación en vídeo MPEG 'Programa completo. Que tenga un buen día :)' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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, 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, 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, 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
Dependiendo de las condiciones iniciales y del tiempo de ejecución, podría ser conveniente elegir un recinto aún mayor, con multiplicador superior a tres (se sugiere ocho). De no hacerlo así podría suceder que el retorno debido al rebote de la ondas en los contornos del recinto distorsionara las condiciones de contorno para el recinto pequeño, anulando así la hipótesis de infinitud para el dominio. MIM-UCM Página 15
El código fuente para esta aproximación se ha parametrizado para observar el comienzo del rebote de 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.
Autor: Alfonso de la Fuente Ruiz http://www.linkedin.com/in/alfonsofr
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) *5+: “Computational Physics 430 (course)” by Ross Spencer and Michael Ware, Department of Physics and Astronomy, Brigham Young University. http://www.physics.byu.edu/Courses/Computational/phys430/phys430.pdf [6]: “Tiny FDTD v1.0” by Nick Clark. MATLAB Central (7 Aug. 2008) http://www.mathworks.com/matlabcentral/fileexchange/21000 *7+: “Web personal Ana Carpio: Laboratorio de Realidad Virtual – Ecuaciones de ondas (Ejemplos VRML en Cortona3D)” por Ana Carpio Rodríguez y Guillermo Garrido Yuste, UCMUNED. http://www.mat.ucm.es/~acarpio/ondas.htm *8+: “Notes on perfectly matched layers (PMLs)” by Steven G. Johnson. MIT (19 July 2008) http://math.mit.edu/~stevenj/18.369/pml.pdf *9+: “Electromagnetic simulation using the FDTD method” by Dennis M. Sullivan, University of Idaho. IEEE press (2000) [10]: “Finite Difference Methods for Differential Equations” by Randall J. LeVeque. 585 course from the University of Washington (January 2006) http://www.amath.washington.edu/~rjl/pubs/am58X/am585winter06.pdf *11+: “UNIVERSIDAD AUTÓNOMA DE MADRID: Prácticas de Ecuaciones en Derivadas Parciales y Análisis Numérico” de Pablo Angulo Ardoy y María López Fernández (UAM) http://www.uam.es/personal_pdi/ciencias/pangulo/
Más de este usuarioAcreditación de Interés - Premios Doctor t-cue 2015El error de KantResumen y conclusiones de la tesis doctoral de Alfonso de la Fuente RuizAtentados Del 11-MSummarized Curriculum Vitae in English 2015MULTIOBJECTIVE OPTIMIZATION OF ALTERNATIVE ROUTESProgramación Matemática con SAS®, FORTRAN-NAG y MPLTécnicas de Análisis Multivariante en SASResumen ponencia I Jornada Doctorandos UBUEntrevista Alfonso de La Fuente Transcripción-Audio Artículos selectos de D. Ernesto Ruiz y González de LinaresRecurso de reposiciónIntroducción al Chino Mandarín - Eibur CefyntAhorro energético en alumbrado públicoQuantum annealing algorithm descriptionQuantum annealing, descripción del método algorítmicoSesión demostrativa con Mathematica 8Temple paralelo - Estado del arteSituación laboral de ingeniería en informática vF2VotarCongreso RSME 2011 - Sesion informática matematicasSiempre Se Van Los MejoresFuga de cerebros ¿Siempre se van los mejores?Entrevista a Ángel Gil Cuevas en InnovatiaPanorama de la industria del ocio digital - INNOVATIA
Simulación de ondas con diferencias finitas (Matlab)Tr by Alfonso4,6K viewsEmbedDownloadDescriptionTrabajo de simulación de propagación de ondas por un dominio plano infinito truncado. Aproximación mediante diferencias finitas en Matlab. Realizado por Alfonso de la Fuente Ruiz para el Master Ofi...Trabajo de simulación de propagación de ondas por un dominio plano infinito truncado. Aproximación mediante diferencias finitas en Matlab. Realizado por Alfonso de la Fuente Ruiz para el Master Oficial en Ingeniería Matemática de la Universidad Complutense de Madrid.Categories: Types, School Work, Essays & ThesesRead on Scribd mobile: iPhone, iPad and Android.Copyright: Attribution Non-Commercial (BY-NC)Download as PDF, TXT or read online from ScribdFlag for inappropriate contentShow moreShow less
RelatedDiferencias finitasby Saul Perez Perez3º MetMatII - Solución de la ecuación de ondas bidimensional mediante métodos numéricosby Eloy Mora VargasOnPlanas en Matlabby Shelene Guzman GarzaEcuacion de Laplace en Matlabby LeninlmvMÉTODOS NUMÉRICOS EN MATLABby eLeCienciaMetodo de Runge Kuttaby José Manuel Limachi Chávezmatlab con metodos numericosby Sol MirandaIMPLEMENTACI´ON DE LOS M´ETODOS DE diferencias finitasby sam19740Ecuación de la membrana vibratoria by Klever GuamanMétodo de Gauss-y-de-Gauss-Jordanby Vilma FajardoMétodos numéricos básicos para ingenieríaby Carlos Armando De CastroDiferencias Finitas en Transferencia de Calorby Francisco TipanMetodos Numericos Para Ingenieros Quimicos Con Matlabby mdmorelloOptica Hechtby analopez0014724Ecuaciones Usando Matlabby Daniel Zuluaga Meza3.12 Funcion Delta Dirac.by -Misael Canela Román-Programa en matlab del método de Eulerby Marcela MorenoAcreditación de Interés - Premios Doctor t-cue 2015by AlfonsoEl error de Kantby AlfonsoAtentados Del 11-Mby AlfonsoMULTIOBJECTIVE OPTIMIZATION OF ALTERNATIVE ROUTES: A PAN-HEURISTIC PRIMER WITH HAZMAT APPLICATIONSby AlfonsoMetodos Numericos Con Matlabby JoseMendozaResumen y conclusiones de la tesis doctoral de Alfonso de la Fuente Ruizby AlfonsoMetodos Numericos Para Ingenieros 5eby Nomar Esahid Duno SanchezMetodos Numericos Para Ecuaciones Diferencialesby Roman DominguezGUÍA 5 CURVAS PLANASby Dany JimpiECUACIONES DIFERENCIALESby Dev Saez AyalaEjercicios Euler Mejoradoby Ivonne CasasSimilar to Simulación de ondas con diferencias finitas (Matlab)TrDiferencias finitas3º MetMatII - Solución de la ecuación de ondas bidimensional mediante métodos numéricosOnPlanas en MatlabEcuacion de Laplace en MatlabMÉTODOS NUMÉRICOS EN MATLABMetodo de Runge Kuttamatlab con metodos numericosIMPLEMENTACI´ON DE LOS M´ETODOS DE diferencias finitasEcuación de la membrana vibratoria Método de Gauss-y-de-Gauss-JordanMétodos numéricos básicos para ingenieríaDiferencias Finitas en Transferencia de CalorMetodos Numericos Para Ingenieros Quimicos Con MatlabOptica HechtEcuaciones Usando Matlab3.12 Funcion Delta Dirac.Programa en matlab del método de EulerAcreditación de Interés - Premios Doctor t-cue 2015El error de KantAtentados Del 11-MMULTIOBJECTIVE OPTIMIZATION OF ALTERNATIVE ROUTESMetodos Numericos Con MatlabResumen y conclusiones de la tesis doctoral de Alfonso de la Fuente RuizMetodos Numericos Para Ingenieros 5eMetodos Numericos Para Ecuaciones DiferencialesGUÍA 5 CURVAS PLANASECUACIONES DIFERENCIALESEjercicios Euler MejoradoSismica de ReflexionDiferencias Finitas-Ecuaciones Elipticas

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución