Ecuación de ondas (grupo 2B)

De MateWiki
Revisión del 22:42 19 may 2014 de Araceli Martin (Discusión | contribuciones) (→‎Sujeción a un aparato que responde a una vibración externa.)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Ecuación de ondas. Grupo 2-B
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores

Ignacio Díaz-Caneja Camblor

Alberto Fernández Pérez

Adela González Barbado

Lucia López Sánchez

Araceli Martín Candilejo

Diego Solano López

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Modelización de los desplazamientos del Cable.

El problema que se tratará a continuación describe una ecuación de tipo hiperbólico, una ecuación de ondas. Consideramos un cable de 10 metros de longitud sujeto por sus extremos. Suponiendo que éste tiene una sección pequeña respecto a su longitud someteremos al cable a pequeñas vibraciones que estudiaremos con una modelización de la ecuación de ondas. Se caracteriza al cable de una masa constante por unidad de volumen, es decir, será homogéneo. Éste, flexible a la tracción, únicamente ofrece resistencias en su dirección longitudinal, tangenciales, pero no a esfuerzos de flexión o cortes. Se supone, además, que es inextensible por lo que será lo suficientemente rígido longitudinalmente como para poder despreciar su extensibilidad.


Con esta modelización se estudiarán las pequeñas vibraciones, desplazamientos transversales, a las que es sometido el cable. Para iniciar el movimiento del cable lo sujetaremos por el centro subiéndolo dos metros, lo que nos proporciona una primera condición inicial determinada por la función g(x), y soltándolo con una velocidad nula.


Obtenemos la ecuación [math]u_{tt}-u_{xx}=0[/math] , de la cual se puede observar que tanto el módulo de la tensión como la densidad son constantes ya que el módulo de la elasticidad es constante, [math]c^2=1[/math] . Al ser ésta primera ecuación homogénea deducimos que transversalmente no actúa ninguna fuerza densidad espacial y temporal. Se deduce del enunciado anterior que las condiciones de frontera son homogéneas, [math]u(0,t)=0[/math] y [math]u(10,t)=0[/math] que nos indican que los extremos del cable no describen ningún tipo de movimiento en todo instante t al no estar sometidos a ningun tipo de fuerza. Finalmente, añadimos al sistema las condiciones iniciales, la primera de ellas, [math]u(x,0)=g(x)= 2- |\frac{2}{5}x -2|[/math], define el movimiento descrito por la onda en el instante t=0 en cualquier punto del cable; y la segunda, [math]u_t(x,0)=0[/math], que la onda parte de una velocidad inicial nula.

Deslizamiento vertical del cable.


[math] \begin{cases} u_{tt}- u_{xx}=0 \ x∈[0,10] \ t\gt0\\ u(0,t)=0\\ u(10,t)=0\\ u(x,0)= 2- |\frac{2}{5}x -2|\\ u_t(x,0)=0 \end{cases} [/math]

En los siguientes apartados se realizaran aproximaciones del sistema expuesto con tres métodos distintos de los que más tarde se hará una comparativa. Para esto último nos ayudaremos del error por lo que necesitaremos la solución analítica. Como el problema consta de condiciones de frontera homogéneas, aplicamos directamente el método de separación de variables para su resolución analítica.

1.1 Aproximación por el método del trapecio.

El primer método que usaremos será el del Trapecio que consiste en RELLENAR 

%u_tt-u_xx=0
%u(0,t)=u(l,t)=0
%u(x,0)=2-abs(0.4x-2)
%ut(x,0)=0
clear all

%Datos del problema
L=10;
T=40;

%Discretización
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;
dt=h;
t=0:dt:T;

%Aprox. de u_xx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1);
K=K/h^2;
F=zeros(N-1,1);

%Posición inicial.
u0=(2-abs(xint*2/5-2))';
v0=zeros(N-1,1);

%Aprox. de u_tt
M=[zeros(N-1),eye(N-1);-K,zeros(N-1)];
G=[zeros(N-1,1);F];

%Aproximación por el método del Trapecio
W0=[u0;v0];
WW=W0;
sol=zeros(length(t),N+1);
sol(1,:)=[0,W0(1:N-1)',0];
for j=1:length(t)-1
    WW=(eye(2*N-2)-dt/2*M)\(eye(2*N-2)+dt/2*M)*WW;
    sol(j+1,:)=[0,WW(1:N-1)',0];
end

%Dibujo de la solución
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);
xlabel('Desplazamiento horizontal');
ylabel('Tiempo');
zlabel('Desplazamiento vertical');


Aproximación por el método del Trapecio.

1.2 Aproximación por el método de Euler.

Se exponen en este apartados dos tipos de aproximaciones de Euler.

Este método se basa en la posibilidad de aproximar una función original mediante rectas tangentes partiendo de un punto dado. Sabiendo que la recta tangente de una función en un punto es la derivada en dicho punto, este método consiste en calcular derivadas de la función en diferentes puntos pertenecientes a ella. Estos puntos serán tomados con la misma distancia a lo largo de la recta, a esta distancia la denominamos paso de discretización.

Aproximación de una función por el método de Euler


1.2.1 Euler Explícito.
%utt-uxx=0
%u(0,t)=u(l,t)=0
%u(x,0)=2-abs(0.4x-2)
%ut(x,0)=0
clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;
dt=h;
t=0:dt:T;
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1);
K=K/h^2;
F=zeros(N-1,1);
u0=2-abs((xint*2/5-2))';
v0=zeros(N-1,1);
M=[zeros(N-1),eye(N-1);-K,zeros(N-1)];
G=[zeros(N-1,1);F];
W0=[u0;v0];
WW=W0;
sol=zeros(length(t),N+1);
sol(1,:)=[0,W0(1:N-1)',0];
for j=1:length(t)-1
    WW=WW+dt*M*WW;
    sol(j+1,:)=[0,WW(1:N-1)',0];
end
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);
1.2.2 Euler Modificado.

ELEGIR UNA DE LAS DOS DEFINICIONES

Este método se basa en la misma idea del método anterior, pero hace un refinamiento en la aproximación, tomando un promedio entre las pendientes.

El método de Euler modificado toma como base que es más razonable pensar que obtendríamos un valor más preciso si aproximáramos tomando un promedio entre las pendientes, en vez de tomarla igual a su valor como en el método anterior de Euler explícito. 

%utt-uxx=0
%u(0,t)=u(l,t)=0
%u(x,0)=2-abs(0.4x-2)
%ut(x,0)=0
clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;
dt=h;
t=0:dt:T;
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1);
K=K/h^2;
F=zeros(N-1,1);
u0=2-abs((xint*2/5-2))';
v0=zeros(N-1,1);
M=[zeros(N-1),eye(N-1);-K,zeros(N-1)];
G=[zeros(N-1,1);F];
W0=[u0;v0];
WW=W0;
sol=zeros(length(t),N+1);
sol(1,:)=[0,W0(1:N-1)',0];
for j=1:length(t)-1
    k1=M*WW;
    WW=(eye(2*N-2)+dt*M/2)*WW+(dt/2)*k1+(dt^2)*M*k1/2;
    sol(j+1,:)=[0,WW(1:N-1)',0];
end
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);


1.3 Aproximación por Fourier con diferentes términos de series.

En este apartado se ofrece una aproximación alternativa a las anteriores. Se muestran cinco iteraciones distintas del método del Fourier que se compararan entre sí y con los métodos anteriores, Trapecio y Euler.
Puesto que la realización de este método es la misma sea cual sea el número de serie elegido, se han expuesto únicamente dos casos de Q para una mayor claridad del programa. 

%utt-uxx=0
%u(0,t)=u(l,t)=0
%u(x,0)=2-abs(0.4x-2)
%ut(x,0)=0

%1 iteracion
clear all
L=10;
T=40;
Q=1;
h=0.1;
N=L/h;
x=0:h:L;
dt=h;
t=0:dt:T;
sol=zeros(length(t),N+1);
for k=1:Q
    phi=sin(k/L*pi*x);
    a=trapz(x,(2-abs(2/5*x-2)).*phi)/(trapz(x,phi.^2));
    b=0;
    T=a.*cos(k*pi*t/L);
    sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
figure (1)
surf(xx,tt,sol)
xlabel('Desplazamiento horizontal');
ylabel('Tiempo');
zlabel('Desplazamiento vertical');

%3 iteracion
clear all
L=10;
T=40;
Q=3;
h=0.1;
N=L/h;
x=0:h:L;
dt=h;
t=0:dt:T;
sol=zeros(length(t),N+1);
for k=1:Q
    phi=sin(k/L*pi*x);
    a=trapz(x,(2-abs(2/5*x-2)).*phi)/(trapz(x,phi.^2));
    b=0;
    T=a.*cos(k*pi*t/L);
    sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
figure (2)
surf(xx,tt,sol)
xlabel('Desplazamiento horizontal');
ylabel('Tiempo');
zlabel('Desplazamiento vertical');



Aproximación por el método del Trapecio con un término de serie.
Aproximación por el método del Trapecio con tres términos de serie.


Aproximación por el método del Trapecio con cinco términos de serie.
Aproximación por el método del Trapecio con diez términos de serie.


Aproximación por el método del Trapecio con veinte términos de serie.



2 Energía del Cable.

La energía del cable que viene definida por la función: \begin{equation} E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx \end{equation}

El método utilizado en el desarrollo numérico de está expresión es el método de diferencias finitas (Funciona calculando de manera aproximada las soluciones a las ecuaciones diferenciales usando ecuaciones diferenciales finitas para aproximar derivadas.).

FUNCIÓN: 

function v = deriva(x,y,t)

n = length(x);
del=zeros(1,n);

if t==1
    
    i=1;
    del(1)=(y(i)-y(i+1))/(-x(i)+x(i+1));

    for i=2:n
    del(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
    end 
    
elseif t==-1

    for i=1:n-1
    del(i)=(y(i)-y(i+1))/(x(i)-x(i+1));
    end 
    
    i=n;
    del(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
    
elseif t==0
    
    i=1;
    del(1)=(y(i)-y(i+1))/(+x(i)-x(i+1));
    
    for i=2:n-1
    del(i)=(y(i+1)-y(i-1))/(2*(x(i)-x(i-1)));
    end 
    
    i=n;
    del(n)=(y(i)-y(i-1))/(x(i)-x(i-1));  
    
end

v=del;

end



PROGRAMA: 

clc;clear all;
L=10;dx=0.1;N=L/dx;
x=dx:dx:L-dx;
xtot=0:dx:L;
T=40;dt=0.1;t=0:dt:T;
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/dx^2)*K;
U=(2-abs((x*2/5-2)))';
V=(0*x)';
sol(1,:)=[0 U' 0];
for i=1:length(t)-1;
    U=U+dt*V;
    V=V-dt*K*U;
    sol(i+1,:)=[0 U' 0];
end
%DERIVAMOS LA SOLUCION EN x Y EN t:
dsol_dx=ones(size(sol));
for i=1:length(t)
    dsol_dx(i,:)=deriva(sol(i,:),xtot,0);
end

dsol_dt=ones(size(sol));
for i=1:length(xtot)
    dsol_dt(:,i)=deriva(sol(:,i),t,1); 
end
%SUMAMOS LAS DERIVADAS AL CUADRADO E INTEGRAMOS:
dE=zeros(1,size(t,1));
for i=1:length(xtot)
    dE=dE+(dsol_dt(i,:).^2+dsol_dx(i,:).^2)*dx;
end
xtot=0:dx:L;
[xx,tt]=meshgrid(xtot,t);
figure(1)
surf(xx,tt,dsol_dx);
figure(2)
surf(xx,tt,dsol_dt);
plot(dE,t)


El resultado habría de ser una gráfica que representase la evolución de la energía de la cuerda frente al transcurso del tiempo. No obstante, en nuestra resolución aparece una gráfica en blanco que ajudicamos o a un posible error en el código, o que podemos interpretar como que la energía tiende a infinito.

3 Aplicaciones.

3.1 Sumersión en un medio viscoso.

Suponemos ahora la inmersión del cable en una medio viscoso. Este medio viscoso produce un amortiguamiento en el movimiento del cable, con lo que la anterior ecuación diferencial se convertiría en \begin{equation} u_{tt}-u_{xx}+au_{t}=0 \end{equation} , siendo [math] a [/math] una constante de amortiguamiento propia del medio viscoso. Puesto que la función [math] u(x,t) [/math] se ve afectada ahora por el medio viscoso, y más directamente, por la constante [math] a [/math] que actúa en su nombre, la energía del cable ,que se expresaba como una integral en función de [math] u(x,t) [/math] también se va a ver afectada: \begin{equation} E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx \end{equation}

El obetivo del siguiente código será hallar una representación de la evolución de la energía en función del tiempo a raíz de la modificación de la función [math] u(x,t) [/math]. Para ello asignaremos distintos valores a la constante [math] a [/math]. Éstos van a ser: [math] a=0 [/math] (la energía para este valor resultará igual que la hallada en el apartado anterior) [math] a=1 [/math] [math] a=4 [/math] [math] a=10 [/math] [math] a=100 [/math]'

El procedimiento numérico volverá a seguir el método de diferencias finitas:

PROGRAMA: 

clc;clear all;close all;
%Datos
dx=0.1;L=10;
dt=0.001;NT=40;
A=[0,1,4,10,100];
x=0:dx:L;
t=0:dt:NT;
ENE=zeros(5,length(t));

%Condición de contorno: U(0,t)=U(L,t)=0;

for k=1:5

%Vector solución y derivada temporal de la solución
%t=0<>T(1) es donde se guardan las condiciones iniciales

U=zeros(length(t),length(x));
Dt_U=zeros(length(t),length(x));
Dx_U=zeros(length(t),length(x));    

%Condición inicial: Estiramos 2 m del centro
%Modelizado como dos rectas antisimétricas que se cortan en x=5 y
%sujetas en x=0 y x=L.

U(1,:)=2*sin(x*pi/L);
Dt_U(1,:)=zeros(1,length(x))';
D2t_U=Ecuacion5(U(1,:),x,-1);

%Inicio del esquema numérico
%empieza en 2 ya que el 1 es la cond. inicial

%ESQUEMA NEWMARK

 for i=2:length(t)
    iter=[k,i]
    Ua=U(i-1,:);
    Dt_Ua=Dt_U(i-1,:);
    D2t_Ua=D2t_U;
    %Derivada segunda respecto del tiempo en el instante t, se obtiene de
    %eq. diferencial u_tt=f(x,u,u_x,u_xx)
    D2t_U=Ecuacion5(U(i-1,:),x,t(i))-A(k)*Dt_Ua;
    %Derivada primera
    Dt_U(i,:)=Dt_Ua+0.5*dt*(D2t_Ua+D2t_U);
    %Ecuacion
    U(i,:) = Ua+Dt_Ua*dt+(1/6)*(D2t_Ua+D2t_U)*dt^2;
 end


 %ENERGIA
 for i=1:length(t)
    Dx_U(i,:)=deriva(x,U(i,:),0);
 end
dE=zeros(length(t),1);
for i=1:length(x)
    dE=dE+(Dt_U(:,i).^2+Dx_U(:,i).^2)*dx;
end
 ENE(k,:)=dE;
end

 figure(1)
 hold on
 grid
 aa=find(x==L/4);
 ab=find(x==3*L/4);
 legend('A=0','A=1','A=4','A=10','A=100');
 plot(t,ENE(1,:),t,ENE(2,:),t,ENE(3,:),t,ENE(4,:),t,ENE(5,:));


FUNCIONES: DERIVA: 

function v = deriva(x,y,t)
n = length(x);
del=zeros(1,n);
if t==1
    i=1;
    del(1)=(y(i)-y(i+1))/(-x(i)+x(i+1));
    for i=2:n
    del(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
    end 
elseif t==-1
    for i=1:n-1
    del(i)=(y(i)-y(i+1))/(x(i)-x(i+1));
    end 
    i=n;
    del(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
elseif t==0
    i=1;
    del(1)=(y(i)-y(i+1))/(+x(i)-x(i+1));
    for i=2:n-1
    del(i)=(y(i+1)-y(i-1))/(2*(x(i)-x(i-1)));
    end 
    i=n;
    del(n)=(y(i)-y(i-1))/(x(i)-x(i-1));  
end
v=del;
end


ECUACIÓN 5: 

function [ V ] = Ecuacion5( u,x,t )
%se guarda aquí la ecuación diferencial del tipo u_tt=f(x,u,u_x,u_xx)
%INCLUIDAS LAS CC AQUÍ
%PARA u_tt=u_xx
%CC: u(0,t)=u(L,t)=0;

Ax=x(2)-x(1);
n=length(u);
V=zeros(1,n);

%Las condiciones de contorno vienen en x(1)=0 y x(Nx)=L
u(1)=0;

if(t==-1)
    u(n)=0;
else
    u(n)=0;
    %u(n)=sin(2*pi*F0*t);
    %u(n)=(1/(b*Ax-1))*u(n-1);
end
for i=2:n-1
    
    V(i)=(1/Ax^2)*(u(i+1)-2*u(i)+u(i-1));
end

i=1;
V(i)=(1/Ax^2)*(u(i+2)-2*u(i+1));
i=n;
V(i)=(1/Ax^2)*(u(i)-2*u(i-1)+u(i-2));
end


La gráfica buscada es:

Gráfica de la energía(ordenadas) del cable respecto del tiempo (Abcisas).

En la gráfica se aprecian cuatro curvas de distintos colores correspondientes a [math] a=0 [/math]',[math] a=1 [/math],[math] a=4 [/math],[math] a=10 [/math],[math] a=100 [/math].

Como podemos apreciar el valor de la energía disminuye a medida que aumentamos el valor de "a". Esto es debido a que "a" define el coeficiente de amortiguamiento del líquido en el cual se encuentra la cuerda. Luego, dicha disminución de energía se produce por la disipación en el medio de gran parte de la energía inicial.

3.2 Sujeción a una estructura de vibración periódica.

En este caso supondremos que el extremo derecho que sujeta el cable está sometido a una vibración periódica de frecuencia [math] F_0 [/math] [Herzios]. Al haber una vibración externa que afecta directamente al cable la expresión de su movimiento y, por consecuencia, de su energía, se va a ver modificada. Nosotros incluiremos esta modificación en la función [math] f(t) [/math] de problema de valor inicial que venimos tratando. Esta función va a pasar a ser: \begin{equation} f(t)= sin(2 \pi F_0t) \end{equation}

El objetivo será el de ver cómo afecta esta alteración en el desarrollo de la energía del cable frente al tiempo. Para nuestro cálculo numérico especificaremos para las siguientes situaciones: En primer lugar, para: [math] F_0= \frac{1}{L} + 0,01 [/math] Después,para: [math] F_0= \frac{1}{L} - 0,01 [/math]

Y para finalizar, lo estudiaremos en un intervalo de tiempo de t=0 hasta t=60 : [math]t \in (0,60)[/math] [segundos].

El código matemático para matlab / octave:


PROGRAMA: 

clc;clear all;close all;

%Datos

dx=0.1;L=10;dt=0.0001;NT=60;
A=[0,1,4,10,100];
x=0:dx:L;t=0:dt:NT;
ENE=zeros(5,length(t));

%Condición de contorno: U(0,t)=U(L,t)=0;
for k=1
%Vector solución y derivada temporal de la solución
%t=0<>T(1) es donde se guardan las condiciones iniciales

U=zeros(length(t),length(x));
Dt_U=zeros(length(t),length(x));
Dx_U=zeros(length(t),length(x));

%Condición inicial: Estiramos 2 m del centro
%Modelizado como dos rectas antisimétricas que se cortan en x=5 y
%sujetas en x=0 y x=L.

U(1,:)=2*sin(x*pi/L);
Dt_U(1,:)=zeros(1,length(x))';
D2t_U=Ecuacion6(U(1,:),x,-1);

%Inicio del esquema numérico
%empieza en 2 ya que el 1 es la cond. inicial

%ESQUEMA NEWMARK

 for i=2:length(t)
    iter=[k,i]
    Ua=U(i-1,:);
    Dt_Ua=Dt_U(i-1,:);
    D2t_Ua=D2t_U;
    %Derivada segunda respecto del tiempo en el instante t, se obtiene de
    %eq. diferencial u_tt=f(x,u,u_x,u_xx)
    D2t_U=Ecuacion6(U(i-1,:),x,t(i))-A(k)*Dt_Ua;
    %Derivada primera
    Dt_U(i,:)=Dt_Ua+0.5*dt*(D2t_Ua+D2t_U);
    %Ecuacion
    U(i,:) = Ua+Dt_Ua*dt+(1/6)*(D2t_Ua+D2t_U)*dt^2;
 end


 %ENERGIA


 for i=1:length(t)
    Dx_U(i,:)=deriva(x,U(i,:),0);
 end
dE=zeros(length(t),1);
for i=1:length(x)
    dE=dE+(Dt_U(:,i).^2+Dx_U(:,i).^2)*dx;
end
 ENE(k,:)=dE;
end
 figure(1)
 hold on
 grid
 aa=find(x==L/4);
 ab=find(x==3*L/4);
 imagesc(U)


FUNCIONES:

DERIVA: 

function v = deriva(x,y,t)
n = length(x);
del=zeros(1,n);
if t==1
    i=1;
    del(1)=(y(i)-y(i+1))/(-x(i)+x(i+1));
    for i=2:n
    del(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
    end 
elseif t==-1
    for i=1:n-1
    del(i)=(y(i)-y(i+1))/(x(i)-x(i+1));
    end 
    i=n;
    del(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
elseif t==0
    i=1;
    del(1)=(y(i)-y(i+1))/(+x(i)-x(i+1));
    for i=2:n-1
    del(i)=(y(i+1)-y(i-1))/(2*(x(i)-x(i-1)));
    end 
    i=n;
    del(n)=(y(i)-y(i-1))/(x(i)-x(i-1));  
end
v=del;
end

ECUACIÓN 6:

function [ V ] = Ecuacion6( u,x,t ) 

%se guarda aquí la ecuación diferencial del tipo u_tt=f(x,u,u_x,u_xx)
%INCLUIDAS LAS CC AQUÍ!
%PARA u_tt=u_xx
%CC: u(0,t)=u(L,t)=0;

Ax=x(2)-x(1);
n=length(u);
V=zeros(1,n);
F0=(1/x(n))+0.01;
b=-2;

%Las condiciones de contorno vienen en x(1)=0 y x(Nx)=L

u(1)=0;
if(t==-1)
    u(n)=0;
else
    %u(n)=0;
    u(n)=sin(2*pi*F0*t);
    %u(n)=(1/(b*Ax-1))*u(n-1);
end
for i=2:n-1
    V(i)=(1/Ax^2)*(u(i+1)-2*u(i)+u(i-1));
end
i=1;
V(i)=(1/Ax^2)*(u(i+2)-2*u(i+1));
i=n;
V(i)=(1/Ax^2)*(u(i)-2*u(i-1)+u(i-2));
end
Gráfica en tres dimensiones de la energía de un cable cuyo extremo derecho está sometido a una vibración externa periódica.

function [ V ] = Ecuacion6( u,x,t ) 

%se guarda aquí la ecuación diferencial del tipo u_tt=f(x,u,u_x,u_xx)
%INCLUIDAS LAS CC AQUÍ!
%PARA u_tt=u_xx
%CC: u(0,t)=u(L,t)=0;

Ax=x(2)-x(1);
n=length(u);
V=zeros(1,n);
F0=(1/x(n))-0.01;
b=-2;

%Las condiciones de contorno vienen en x(1)=0 y x(Nx)=L

u(1)=0;
if(t==-1)
    u(n)=0;
else
    %u(n)=0;
    u(n)=sin(2*pi*F0*t);
    %u(n)=(1/(b*Ax-1))*u(n-1);
end
for i=2:n-1
    V(i)=(1/Ax^2)*(u(i+1)-2*u(i)+u(i-1));
end
i=1;
V(i)=(1/Ax^2)*(u(i+2)-2*u(i+1));
i=n;
V(i)=(1/Ax^2)*(u(i)-2*u(i-1)+u(i-2));
end
Gráfica en tres dimensiones de la energía de un cable cuyo extremo derecho está sometido a una vibración externa periódica.

function [ V ] = Ecuacion6( u,x,t ) 

%se guarda aquí la ecuación diferencial del tipo u_tt=f(x,u,u_x,u_xx)
%INCLUIDAS LAS CC AQUÍ!
%PARA u_tt=u_xx
%CC: u(0,t)=u(L,t)=0;

Ax=x(2)-x(1);
n=length(u);
V=zeros(1,n);
F0=(1/x(n));
b=-2;

%Las condiciones de contorno vienen en x(1)=0 y x(Nx)=L

u(1)=0;
if(t==-1)
    u(n)=0;
else
    %u(n)=0;
    u(n)=sin(2*pi*F0*t);
    %u(n)=(1/(b*Ax-1))*u(n-1);
end
for i=2:n-1
    V(i)=(1/Ax^2)*(u(i+1)-2*u(i)+u(i-1));
end

i=1;
V(i)=(1/Ax^2)*(u(i+2)-2*u(i+1));
i=n;
V(i)=(1/Ax^2)*(u(i)-2*u(i-1)+u(i-2));
end


La gráfica resultante es la siguiente:

Gráfica en tres dimensiones de la energía de un cable cuyo extremo derecho está sometido a una vibración externa periódica.


En dicha gráfica, se puede observar que la forma de la función o la energía se ve alterada por la condición impuesta de tal manera que ésta adquiere una forma también sinusoidal, y que su estabilización dependería del tiempo transcurrido entre vibración y vibración.

3.3 Sujeción a un aparato que responde a una vibración externa.

La siguiente situación que vamos a representar a modelar a través de ecuaciones diferenciales va a ser la siguiente: supondremos que el cable de la anterior situación (que recibía en su extremo derecho una vibración periódica de [math] F_0 [/math] [Herzios])envía una respuesta a la vibración que recibe. Esta situación la vamos a representar como un cambio en la condición de contorno; ésta va a pasar a ser: [math] u_x(L,t)=bu(L,t) [/math]

El objetivo, de nuevo será representar la evolución de la energía respecto del tiempo y ver si este cambio en el problema se ve reflejado en su representación. En el código numérico, no obstante, la situación de [math] u_x(L,t)=bu(L,t) [/math] la vamos a aproximar por la siguiente expresión: [math] u_x(L,t)-bu(L,t) \cong \frac{u_{N+1}(t)-u_{N-1}(t)}{2h-u_N(t)} [/math]

En la resolución del problema emplearemos de nuevo la condición inicial de sujetar el cable por el centro y elevarlo 2 m perpendicularmente: [math] u(x,0)=g(x)= 2- |\frac{2}{5}x -2| [/math]

Y además obtendremos el resultado para los dos valores siguientes de [math]b[/math] : Para [math] b=2 [/math] y [math] b=-2 [/math].

El código numérico de matlab/octave es:

PROGRAMA: 

clc;clear all;close all;

%Datos
dx=0.1;L=10;dt=0.0001;NT=40;
A=[0,1,4,10,100];
x=0:dx:L;t=0:dt:NT;
ENE=zeros(5,length(t));

%Condición de contorno: U(0,t)=U(L,t)=0;
for k=1

%Vector solución y derivada temporal de la solución
%t=0<>T(1) es donde se guardan las condiciones iniciales

U=zeros(length(t),length(x));
Dt_U=zeros(length(t),length(x));
Dx_U=zeros(length(t),length(x));    

%Condición inicial: Estiramos 2 m del centro
%Modelizado como dos rectas antisimétricas que se cortan en x=5 y
%sujetas en x=0 y x=L.

U(1,:)=2*sin(x*pi/L);
Dt_U(1,:)=zeros(1,length(x))';
D2t_U=Ecuacion7(U(1,:),x,-1);

%Inicio del esquema numérico
%empieza en 2 ya que el 1 es la cond. inicial


%ESQUEMA NEWMARK
 for i=2:length(t)
    iter=[k,i]
    Ua=U(i-1,:);
    Dt_Ua=Dt_U(i-1,:);
    D2t_Ua=D2t_U;

    %Derivada segunda respecto del tiempo en el instante t, se obtiene de
    %eq. diferencial u_tt=f(x,u,u_x,u_xx)

    D2t_U=Ecuacion7(U(i-1,:),x,t(i))-A(k)*Dt_Ua;

    %Derivada primera
    Dt_U(i,:)=Dt_Ua+0.5*dt*(D2t_Ua+D2t_U);

    %Ecuacion
    U(i,:) = Ua+Dt_Ua*dt+(1/6)*(D2t_Ua+D2t_U)*dt^2;
 end


 %ENERGIA
 for i=1:length(t)
    Dx_U(i,:)=deriva(x,U(i,:),0);
 end
dE=zeros(length(t),1);
for i=1:length(x)
    dE=dE+(Dt_U(:,i).^2+Dx_U(:,i).^2)*dx;
end
 ENE(k,:)=dE;
end

 figure(1)
 hold on
 grid
 aa=find(x==L/4);
 ab=find(x==3*L/4);
 plot(t,ENE(1,:))


FUNCIONES: function v = deriva(x,y,t) 

n = length(x);
del=zeros(1,n);
if t==1
    i=1;
    del(1)=(y(i)-y(i+1))/(-x(i)+x(i+1));
    for i=2:n
    del(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
    end 
elseif t==-1
    for i=1:n-1
    del(i)=(y(i)-y(i+1))/(x(i)-x(i+1));
    end 
    i=n;
    del(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
elseif t==0
    i=1;
    del(1)=(y(i)-y(i+1))/(+x(i)-x(i+1));
    for i=2:n-1
    del(i)=(y(i+1)-y(i-1))/(2*(x(i)-x(i-1)));
    end 
    i=n;
    del(n)=(y(i)-y(i-1))/(x(i)-x(i-1));  
end
v=del;
end


function [ V ] = Ecuacion7( u,x,t ) 

%se guarda aquí la ecuación diferencial del tipo u_tt=f(x,u,u_x,u_xx)
%INCLUIDAS LAS CC AQUÍ!
%PARA u_tt=u_xx
%CC: u(0,t)=u(L,t)=0;

Ax=x(2)-x(1);
n=length(u);
V=zeros(1,n);
F0=(1/x(n));
b=2;

%Las condiciones de contorno vienen en x(1)=0 y x(Nx)=L
u(1)=0;
if(t==-1)
    u(n)=0;
else
    %u(n)=0;
    %u(n)=sin(2*pi*F0*t);
    u(n)=(1/(b*Ax-1))*u(n-1);
end
for i=2:n-1
    V(i)=(1/Ax^2)*(u(i+1)-2*u(i)+u(i-1));
end

i=1;
V(i)=(1/Ax^2)*(u(i+2)-2*u(i+1));
i=n;
V(i)=(1/Ax^2)*(u(i)-2*u(i-1)+u(i-2));
end
Gráfica de la energía (ordenadas) de un cable frente al tiempo (abcisas) que recibe una vibración en su extremo derecho y que responde a dicha vibración.

function [ V ] = Ecuacion7( u,x,t ) 

%se guarda aquí la ecuación diferencial del tipo u_tt=f(x,u,u_x,u_xx)
%INCLUIDAS LAS CC AQUÍ!
%PARA u_tt=u_xx
%CC: u(0,t)=u(L,t)=0;

Ax=x(2)-x(1);
n=length(u);
V=zeros(1,n);
F0=(1/x(n));
b=-2;

%Las condiciones de contorno vienen en x(1)=0 y x(Nx)=L
u(1)=0;

if(t==-1)
    u(n)=0;
else
    %u(n)=0;
    %u(n)=sin(2*pi*F0*t);
    u(n)=(1/(b*Ax-1))*u(n-1);
end

for i=2:n-1
    V(i)=(1/Ax^2)*(u(i+1)-2*u(i)+u(i-1));
end

i=1;
V(i)=(1/Ax^2)*(u(i+2)-2*u(i+1));
i=n;
V(i)=(1/Ax^2)*(u(i)-2*u(i-1)+u(i-2));
end


La imagen resultante es:

Gráfica de la energía (ordenadas) de un cable frente al tiempo (abcisas) que recibe una vibración en su extremo derecho y que responde a dicha vibración.


La anterior gráfica muestra el comportamiento de la cuerda como respuesta a la vibración producida en uno de sus extremos así como el movimiento transversal del punto medio al comienzo. Podemos observar el movimiento transversal de cada punto de la cuerda a través de la expresión de la energía, así como que, a medida que transcurre el tiempo, ésta se va estabilizando. La forma de la función recuerda a una función sinusoidal debido a la forma de la condición de contorno impuesta.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_de_ondas_(grupo_2B)&oldid=13105»