Diferencia entre revisiones de «Ecuación de ondas. Grupo C2»
De MateWiki
(→Cambio de condiciones en los extremos. Cálculo de la energía) |
(→Método del trapecio) |
||
| Línea 11: | Línea 11: | ||
=== Método del trapecio === | === Método del trapecio === | ||
| + | {{matlab|codigo= | ||
| + | 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)'); | ||
| + | }} | ||
=== Método de Euler explícito === | === Método de Euler explícito === | ||
Revisión del 21:44 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 | |
Contenido
- 1 Introducción
- 2 Vibración sin amortiguamiento. Condiciones Dirichlet. Resolución por el método de líneas
- 3 Vibración con amortiguamiento. Cálculo de la energía
- 4 Cambio de condiciones en los extremos. Cálculo de la energía
- 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
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
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
2.3 Método de Heun
2.4 Cálculo de la energía
3 Vibración con amortiguamiento. Cálculo de la energía
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)');
