Diferencia entre revisiones de «Ecuación de ondas (grupo 2B)»

De MateWiki
Saltar a: navegación, buscar
(→‎Energía del Cable.)
(→‎Sumersión en un medio viscoso.)
Línea 350: Línea 350:
  
 
==== Sumersión en un medio viscoso. ====
 
==== Sumersión en un medio viscoso. ====
 +
clc;clear all;
 +
L=10;dx=0.1;N=L/dx;
 +
x=dx:dx:L-dx;
 +
xtot=0:dx:L;
 +
T=10;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
 +
%BORRAR ESTO!!!!!! HASTA AQUI CREO QUE ESTA MAL EL PROGRAMA
 +
 +
a=0;
 +
sol2=sol+a*sol;
 +
%DERIVAMOS LA SOLUCION EN x Y EN t:
 +
dsol_dx=ones(size(sol2));
 +
for i=1:length(t)
 +
    dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);
 +
end
 +
dsol_dt=ones(size(sol2));
 +
for i=1:length(xtot)
 +
    dsol_dt(:,i)=deriva(sol2(:,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
 +
figure(1)
 +
plot(dE,t)
 +
 +
a=1;
 +
sol2=sol+a*sol;
 +
%DERIVAMOS LA SOLUCION EN x Y EN t:
 +
dsol_dx=ones(size(sol2));
 +
for i=1:length(t)
 +
    dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);
 +
end
 +
dsol_dt=ones(size(sol2));
 +
for i=1:length(xtot)
 +
    dsol_dt(:,i)=deriva(sol2(:,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
 +
figure(2)
 +
plot(dE,t)
 +
 +
a=4;
 +
sol2=sol+a*sol;
 +
%DERIVAMOS LA SOLUCION EN x Y EN t:
 +
dsol_dx=ones(size(sol2));
 +
for i=1:length(t)
 +
    dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);
 +
end
 +
dsol_dt=ones(size(sol2));
 +
for i=1:length(xtot)
 +
    dsol_dt(:,i)=deriva(sol2(:,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
 +
figure(3)
 +
plot(dE,t)
 +
 +
a=10;
 +
sol2=sol+a*sol;
 +
%DERIVAMOS LA SOLUCION EN x Y EN t:
 +
dsol_dx=ones(size(sol2));
 +
for i=1:length(t)
 +
    dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);
 +
end
 +
dsol_dt=ones(size(sol2));
 +
for i=1:length(xtot)
 +
    dsol_dt(:,i)=deriva(sol2(:,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
 +
figure(4)
 +
plot(dE,t)
 +
 +
a=100;
 +
sol2=sol+a*sol;
 +
%DERIVAMOS LA SOLUCION EN x Y EN t:
 +
dsol_dx=ones(size(sol2));
 +
for i=1:length(t)
 +
    dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);
 +
end
 +
dsol_dt=ones(size(sol2));
 +
for i=1:length(xtot)
 +
    dsol_dt(:,i)=deriva(sol2(:,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
 +
figure(5)
 +
plot(dE,t)
 +
}}
  
 
==== Sujeción a una estructura de vibración periódica. ====
 
==== Sujeción a una estructura de vibración periódica. ====

Revisión del 16:09 17 may 2014

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 Deslizamiento Vertical del Cable.

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 a la ecuación de ondas. Se caracteriza al cable de una masa constante por unidad de volumen, es decir, será homogéneo. Éste, flexible e inextensible a la tracción, únicamente ofrece resistencias en su dirección longitudinal, tangenciales, pero no a esfuerzos de flexión o cortes. 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. Al estar el cable sujeto por ambos extremos, las condiciones de frontera serán homogéneas.

Deslizamiento vertical del cable.


[math] \begin{cases} u_tt- u_xx=0\\ u(0,t)=0\\ u(10,t)=0\\ u(x,0)=g(x)= 2- |\frac{2}{5}x -2|\\ u_t(x,0)=0 \end{cases} [/math]

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

%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=(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
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);


Aproximación por el medio del Trapecio.

1.2 Aproximació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.
%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.

%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/2*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)/(k*pi/L);
    sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
figure (1)
surf(xx,tt,sol)

%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/2*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)/(k*pi/L);
    sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
figure (2)
surf(xx,tt,sol)

%5 iteracion
clear all
L=10;
T=40;
Q=5;
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/2*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)/(k*pi/L);
    sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
figure (3)
surf(xx,tt,sol)

%10 iteracion
clear all
L=10;
T=40;
Q=10;
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/2*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)/(k*pi/L);
    sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
figure (4)
surf(xx,tt,sol)

%20 iteracion
clear all
L=10;
T=40;
Q=20;
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/2*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)/(k*pi/L);
    sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
figure (5)
surf(xx,tt,sol)


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.

clc;clear all; L=10;dx=0.1;N=L/dx; x=dx:dx:L-dx; xtot=0:dx:L; T=10;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 %BORRAR ESTO!!!!!! HASTA AQUI CREO QUE ESTA MAL EL PROGRAMA

a=0; sol2=sol+a*sol; %DERIVAMOS LA SOLUCION EN x Y EN t: dsol_dx=ones(size(sol2)); for i=1:length(t) 

   dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);

end dsol_dt=ones(size(sol2)); for i=1:length(xtot) 

   dsol_dt(:,i)=deriva(sol2(:,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 figure(1) plot(dE,t)

a=1; sol2=sol+a*sol; %DERIVAMOS LA SOLUCION EN x Y EN t: dsol_dx=ones(size(sol2)); for i=1:length(t) 

   dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);

end dsol_dt=ones(size(sol2)); for i=1:length(xtot) 

   dsol_dt(:,i)=deriva(sol2(:,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 figure(2) plot(dE,t)

a=4; sol2=sol+a*sol; %DERIVAMOS LA SOLUCION EN x Y EN t: dsol_dx=ones(size(sol2)); for i=1:length(t) 

   dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);

end dsol_dt=ones(size(sol2)); for i=1:length(xtot) 

   dsol_dt(:,i)=deriva(sol2(:,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 figure(3) plot(dE,t)

a=10; sol2=sol+a*sol; %DERIVAMOS LA SOLUCION EN x Y EN t: dsol_dx=ones(size(sol2)); for i=1:length(t) 

   dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);

end dsol_dt=ones(size(sol2)); for i=1:length(xtot) 

   dsol_dt(:,i)=deriva(sol2(:,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 figure(4) plot(dE,t)

a=100; sol2=sol+a*sol; %DERIVAMOS LA SOLUCION EN x Y EN t: dsol_dx=ones(size(sol2)); for i=1:length(t) 

   dsol_dx(i,:)=deriva(sol2(i,:),xtot,0);

end dsol_dt=ones(size(sol2)); for i=1:length(xtot) 

   dsol_dt(:,i)=deriva(sol2(:,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 figure(5) plot(dE,t) }}

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

3.3 Sujeción a un aparato condicionado.

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