Diferencia entre revisiones de «Ecuación de ondas. Grupo C2»

De MateWiki
Saltar a: navegación, buscar
(→‎Método de Heun)
(→‎Método del trapecio)
Línea 11: Línea 11:
  
 
=== Método del trapecio ===
 
=== Método del trapecio ===
 +
[[Archivo:T3Apartado2_graf.jpg|400px|miniaturadeimagen|derecha||Gráfica de la posición de cada punto del cable]]
 
{{matlab|codigo=
 
{{matlab|codigo=
 
clear all
 
clear all

Revisión del 21:50 7 may 2015

Trabajo realizado por estudiantes
Título Ecuación de ondas. Grupo C2
Asignatura Ecuaciones Diferenciales
Curso Curso 2014-15
Autores Nuestros nombres
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Consideramos un cable de una estructura civil de longitud L = 10 m sujeto por ambos extremos. Supondremos que el cable tiene una sección pequeña respecto a su longitud y que las vibraciones pueden modelizarse mediante la ecuación de ondas. Si denotamos su desplazamiento vertical por u(x,t), podemos plantear el problema de su movimiento según el siguiente sistema de ecuaciones:

2 Vibración sin amortiguamiento. Condiciones Dirichlet. Resolución por el método de líneas

2.1 Método del trapecio

Gráfica de la posición de cada punto del cable
clear all
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(2:N);
xx=xx';
ua=0;ub=0; %Condiciones de contorno.
U0=zeros(size(xx)); %Preasignación de U0.
%Recorremos mediante un bucle U0, y añadimos los valores que correspondan.
for j=1:length(xx);
    if xx(j)<b/3
        U0(j)=3*xx(j)/b;
    else
        U0(j)=1.5-1.5*xx(j)/b;
    end
end
V0=zeros(size(xx)); %Preasignación de V0.
%Matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
%Término F y valor inicial
F=0*xx;
F(1)=F(1)+ua/h^2;
F(end)=F(end)+ub/h^2;
%Resolución del sistema de ecuaciones de EDO de orden 1.
t0=0;tM=40;
k=h; %Paso en t.
t=t0:k:tM; %Discretización del vector de tiempos.
M=length(t)-1; %Número de subintervalos.
%Añadimos en la primera columna las condiciones iniciales.
U(:,1)=U0;
V(:,1)=V0;
for i=1:M
    %Resolución del sistema de ecuaciones por el método del trapecio
    U(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(U(:,i)+0.5*k*(2*V(:,i)+0.5*k*(-K*U(:,i)+2*F)));
    V(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(V(:,i)+0.5*k*(-K*U(:,i)+2*F)-0.5*k*K*(U(:,i)+0.5*k*V(:,i)));
end
%Incluimos condiciones Dirichlet.
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
U=[UA;U;UB];
%Dibujamos el gráfico.
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');


2.2 Método de Euler explícito

clear all
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(2:N);
xx=xx';
ua=0;ub=0; %Condiciones de contorno.
U0=zeros(size(xx)); %Preasignación de U0.
%Recorremos mediante un bucle U0, y añadimos los valores que correspondan.
for j=1:length(xx);
    if xx(j)<b/3
        U0(j)=3*xx(j)/b;
    else
        U0(j)=1.5-1.5*xx(j)/b;
    end
end
V0=zeros(size(xx)); %Preasignación de V0.
%Matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
%Término F y valor inicial
F=0*xx;
F(1)=F(1)+ua/h^2;
F(end)=F(end)+ub/h^2;
%Resolución del sistema de ecuaciones de EDO de orden 1.
t0=0;tM=40;
k=h; %Paso en t.
t=t0:k:tM; %Discretización del vector de tiempos.
M=length(t)-1; %Número de subintervalos.
%Añadimos en la primera columna las condiciones iniciales.
U(:,1)=U0;
V(:,1)=V0;
for i=1:M
    %Resolución del sistema de ecuaciones por el método de Euler (explícito).
    U(:,i+1)=U(:,i)+k*V(:,i);
    V(:,i+1)=V(:,i)+k*(-K*U(:,i)+F);
end
%Incluimos condiciones Dirichlet.
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
U=[UA;U;UB];
%Dibujamos el gráfico.
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');


2.3 Método de Heun

clear all
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(2:N);
xx=xx';
ua=0;ub=0; %Condiciones de contorno.
U0=zeros(size(xx)); %Preasignación de U0.
%Recorremos mediante un bucle U0, y añadimos los valores que correspondan.
for j=1:length(xx);
    if xx(j)<b/3
        U0(j)=3*xx(j)/b;
    else
        U0(j)=1.5-1.5*xx(j)/b;
    end
end
V0=zeros(size(xx)); %Preasignación de V0.
%Matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
%Término F y valor inicial
F=0*xx;
F(1)=F(1)+ua/h^2;
F(end)=F(end)+ub/h^2;
%Resolución del sistema de ecuaciones de EDO de orden 1.
t0=0;tM=40;
k=h; %Paso en t.
t=t0:k:tM; %Discretización del vector de tiempos.
M=length(t)-1; %Número de subintervalos.
%Añadimos en la primera columna las condiciones iniciales.
U(:,1)=U0;
V(:,1)=V0;
for i=1:M
    %Resolución del sistema de ecuaciones por el método de Heun
    K1U=V(:,i);
    K2U=V(:,i)+K1U*k;
    K1V=-K*U(:,i)+F;
    K2V=-K*(U(:,i)+K1V*k)+F;
    U(:,i+1)=U(:,i)+(k/2)*(K1U+K2U);
    V(:,i+1)=V(:,i)+(k/2)*(K1V+K2V);
end
%Incluimos condiciones Dirichlet.
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
U=[UA;U;UB];
%Dibujamos el gráfico.
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');


2.4 Cálculo de la energía

3 Vibración con amortiguamiento. Cálculo de la energía

Gráfica de la energía para distintos amortiguamientos
clear all, clf
%Coeficiente de amortiguamiento
am=[0,1,4,10,100];
%Hacemos un bucle donde calcular la energía para cada coeficiente.
for n=am
    %Datos en x
    a=0; b=10; %Longitud del cable L=10.
    h=0.1; %Paso.
    x=a:h:b; %Discretización espacial del cable.
    N=round((b-a)/h);
    %Definimos las matrices de la ecuación
    xx=x(2:N);
    xx=xx';
    ua=0;ub=0; %Condiciones de contorno.
    U0=zeros(size(xx)); %Preasignación de U0.
    %Recorremos mediante un bucle U0, y añadimos los valores que correspondan.
    for j=1:length(xx);
        if xx(j)<b/3
            U0(j)=3*xx(j)/b;
        else
            U0(j)=1.5-1.5*xx(j)/b;
        end
    end
    V0=zeros(size(xx)); %Preasignación de V0.
    %Matriz K
    K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
    %Término F y valor inicial
    F=0*xx;
    F(1)=F(1)+ua/h^2;
    F(end)=F(end)+ub/h^2;
    %Resolución del sistema de ecuaciones de EDO de orden 1
    t0=0;tM=40;
    k=h; %Paso en t.
    t=t0:k:tM; %Discretización del vector de tiempos.
    M=length(t)-1; %Número de subintervalos.
    %Añadimos en la primera columna las condiciones iniciales.
    U(:,1)=U0;
    V(:,1)=V0;
    for i=1:M
        %Sistema de ecuaciones por el método del trapecio
        V(:,i+1)=((1+0.5*k*n)*eye(size(K))+0.25*(k^2)*K)\(V(:,i)+0.5*k*(-n*V(:,i)-K*U(:,i)+2*F)-0.5*k*K*(U(:,i)+0.5*k*V(:,i)));
        U(:,i+1)=U(:,i)+0.5*h*(V(:,i)+V(:,i+1));
    end
    %Incluimos condiciones Dirichlet.
    UA=ua*ones(1,length(t));
    UB=ub*ones(1,length(t));
    U=[UA;U;UB];
    %Como las condiciones Dirichlet son nulas, las velocidades de estos puntos
    %también lo serán
    V=[UA;V;UB];
    %Energía
    E=zeros(size(t)); %Preasignación.
    Ux=zeros(size(x));
    for l=1:M+1
        for m=2:N
            Ux(m)=(U(m+1,l)-U(m-1,l))/(2*k); %Cálculo de la derivada Ux mediante la aproximación por diferencias finitas.
        end
        Ux=Ux';
        E(l)=trapz(x,V(:,l).^2)+trapz(x,Ux.^2); %Cálculo de la energía.
    end
    %Dibujamos la gráfica de la energía.
    hold on
    plot(t,E)
    xlabel('Tiempo (s)'); ylabel('Energía (J)');
    %Borramos todos los datos para realizar el bucle de nuevo.
    clear all
end
legend('a=0','a=1','a=4','a=10','a=100','Location','Best');
hold off


4 Cambio de condiciones en los extremos. Cálculo de la energía

Consideramos que nuestro cable está sujeto a una estructura que sufre vibraciones periódicas con frecuencia F0 Herzios. Vamos a tomar la función [math]f(t)=\sin(2\pi*F0*t)[/math] que será la que defina la posición del extremo izquierdo, que está sujeto a la estructura, en función del tiempo. 

clear all
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(2:N);
xx=xx';
ub=0; %condición de contorno en el extremo derecho.
%Preasignación de la posición y la velocidad incial.
U0=zeros(size(xx));
V0=U0;
%Matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
%Término F y valor inicial
F=0*xx;
%Resolución del sistema de ecuaciones de EDO de orden 1
t0=0;tM=60; %Discretización del vector de tiempos.
k=h; %Paso en t.
t=t0:k:tM;
M=length(t)-1; %Número de subintervalos.
%Añadimos en la primera columna las condiciones iniciales.
U(:,1)=U0;
V(:,1)=V0;
%Pedimos por teclado al usuario los distintos valores de la frecuencia que
%le transmite la estructura al cable. Estas serán F0=1/L+0.01 Hz, 
%F0=1/L-0.01 Hz y F0=1/L Hz, siendo L=b=10.
F0=input('Introduzca la frecuencia (Hz) transmitida al cable: ');
for i=1:M
    F(1)=F(1)+sin(2*pi*F0*t(i))/h^2;
    F(end)=F(end)+ub/h^2;
    %Resolución del sistema de ecuaciones por el método del trapecio.
    U(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(U(:,i)+0.5*k*(2*V(:,i)+0.5*k*(-K*U(:,i)+2*F)));
    V(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(V(:,i)+0.5*k*(-K*U(:,i)+2*F)-0.5*k*K*(U(:,i)+0.5*k*V(:,i)));
end
%Incluimos las condiciones Dirichlet en nuestra solución. 
UA=sin(2*pi*F0*t);
UB=ub*ones(1,length(t));
U=[UA;U;UB];
VB=zeros(1,length(t));
VA=VB;
V=[VA;V;VB];
%Dibujamos el gráfico.
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');
%Cadena de texto con la frecuencia introducida, para el título de la gráfica.
Frec=sprintf('Frecuencia = %.2f Hz',F0); 
title(Frec);
%Energía
E=zeros(size(t)); %Preasignación.
Ux=zeros(size(x));
for l=1:M+1
    for m=2:N
        Ux(m)=(U(m+1,l)-U(m-1,l))/(2*k); %Cálculo de la derivada Ux mediante la aproximación por diferencias finitas.
    end
    Ux=Ux';
    E(l)=trapz(x,V(:,l).^2)+trapz(x,Ux.^2); %Cálculo de la energía.
end
%Dibujamos la gráfica de la energía.
figure(2)
plot(t,E)
xlabel('Tiempo (s)'); ylabel('Energía (J)');
Gráfica de la posición del cable en cada punto y en cada instante.
Gráfica de la posición del cable en cada punto y en cada instante.
Gráfica de la posición del cable en cada punto y en cada instante.

5 Cambio de condiciones en los extremos. Condición Neumann. Cálculo de la energía

6 Vibración sin amortiguamiento. Método de Fourier

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_de_ondas._Grupo_C2&oldid=29609»