Ecuación de ondas. Grupo C2

Trabajo realizado por estudiantes
Título	Ecuación de ondas. Grupo C2
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Ana Martínez Lorente, Javier Parras Martínez, Alfredo Pazos Arjona, Antonio Pérez Mata, Javier Siguero Ginés
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 [math] L = 10 m [/math] 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 [math] u(x,t) [/math], podemos plantear el problema de su movimiento según el siguiente sistema de ecuaciones: [math] \begin{cases} u_{tt}-u_{xx}=f(x,t), \; x∈[0,10], \; t∈[0,T],\\ u(0,t)=g_{0}(t), \; u(L,t)=g_{1}(t),\\ u(x,0)=h_{0}(x), \; u_{t}(x,0)=h_{1}(x)\\ \end{cases} [/math]

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

Vamos a tratar el problema de vibración de un cable de longitud [math] L=10 [/math], en la cual los dos extremos de la misma se encuentran fijos a lo largo del tiempo y con una desplazamiento nulo. Al inicio, en [math] t=0 [/math], sujetamos el cable desde el punto [math] x=L/3 [/math], y lo desplazamos 1 m en la dirección perpendicular a la recta que une sus extremos.

El problema viene modelizado por la siguente ecuación: [math] \begin{cases} u_{tt}-u_{xx}=0, \; x∈[0,L], \; t∈[0,40], \\ u(0,t)=0, \; u(L,t)=0, \\ u(x,0)=\begin{cases} \frac{3x}{L} \\ \frac{3}{2}-\frac{3x}{2L} \end{cases}, \; u_{t}(x,0)=0\\ \end{cases} [/math] A continuación se resuelve el problema por el método de diferencias finitas, aplicando para la resolución de la ecuación matricial que aparece los métodos del trapecio, de Euler explícito y de Heun.

2.1 Método del trapecio

Aplicando el método de de diferencias finitas, o también llamado método de líneas, podemos obtener una solución aproximada del problema propuesto. Como se ha visto en las clases de numérico, al aplicar este método obtenemos la siguiente ecuación matricial a resolver:

[math] \begin{cases} U''=-KU+F\\ U(0)=u^{0}\\ U'(0)=0\\ \end{cases} [/math]

donde [math] K [/math] es la matriz de coeficientes que multiplica a cada [math] u(t) [/math], [math] F [/math] es un vector que sirve para incluir las condiciones Dirichlet de los extremos, [math] u^{0} [/math] la condición inicial de posición y [math] U [/math] es el vector solución de los desplazamientos del cable. Al ser una ecuación diferencial ordinaria de segundo orden, para poder aplicar los métodos numéricos de resolución es necesario pasar a un sistema de ecuaciones diferenciales ordinarias equivalente. Es el siguiente:

[math] \begin{cases} U' = V \\ V' = -KU + F \\ U(0) = u^{0} \\ V(0) = 0 \\ \end{cases} [/math]

donde el nuevo vector [math] V [/math] representa la velocidad de cada punto del cable. Aplicando el método del trapecio a cada ecuación del sistema por separado se obtiene:

[math] \begin{cases} U_{n+1} = U_{n} + \frac{h}{2}(V_{n} + V_{n+1}) \\ V_{n+1} = V_{n} + \frac{h}{2}(-KU_{n} + F_{n} - KU_{n+1} + F_{n+1}) \\ \end{cases} [/math]

y despejando cada variable:

[math] \begin{cases} U_{n+1} (I + \frac{h^2}{4}K) = U_{n} + \frac{h}{2}(2V_{n} + \frac{h}{2}(-KU_{n} + F_{n} + F_{n+1})) \\ V_{n+1} (I + \frac{h^2}{4}K) = V_{n} + \frac{h}{2}(-KU_{n} + F_{n} + F_{n+1}) - \frac{h}{2}K(U_{n} + \frac{h}{2}V_{n}) \\ \end{cases} [/math]

Una vez realizado este proceso analítico, pasamos a implementar el código MatLab/OctaveUPM que resuelve el problema.

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)');

En el gráfico tridimensional podemos observar como varía el desplazamiento vertical en cada punto del cable a lo largo del tiempo. En la parte más cercana al observador podemos apreciar la posición inicial del cable, formando una especie de triángulo, estando el punto de [math] x=L/3 [/math] 1m por encima de la posición horizontal. Cuando se suelta el cable con velocidad cero desde esa posición, el cable tiende a recuperar su posición horizontal, pasando por ella, y alcanzando una posición simétrica con respecto a esta misma horizontal, en la que el punto de [math] x=L/3 [/math] tendrá una desplazamiento vertical negativo de 1m. De nuevo, la cuerda tiende a recuperar su posición horizontal, pasando por ella, y alcanzando otra vez la posición inicial. Al no existir ni amortiguamiento ni ninguna fuerza aplicada, este proceso de oscilación se repite indefinidamente a lo largo del tiempo. Por último, se puede apreciar que ambos extremos del cable tienen desplazamiento vertical nulo a lo largo del tiempo.

2.2 Método de Euler explícito

Ahora vamos a realizar el mismo problema anterior, utilizando el método de Euler explícito. El código MatLab es el siguiente:

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 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)');

Como podemos observar en la gráfica, no tiene nada que ver esta solución a la calculada anteriormente con el método del Trapecio. Esto es debido a que el método de Euler utilizado, a diferencia del método del Trapecio, es un método explícito. Nuestra gráfica aparece de este modo debido a que el método no converge. Se puede apreciar mirando en los valores de la posición de la onda, que son del orden de [math] 10^{136} [/math] , obviamente valores no razonables. Es por esto que no podemos utilizar este método para la resolución de nuestro problema.

2.3 Método de Heun

Realizamos el mismo problema, utilizando el método de Heun. El código MatLab es el siguiente:

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 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)');

Como podemos observar en la gráfica y al igual que en el caso del método de Euler explícito, no tiene nada que ver esta solución a la calculada anteriormente con el método del Trapecio. Esto es debido a que el método de Heun es también un método explícito. Nuestra gráfica aparece de este modo debido a que el método no converge. Se puede apreciar mirando en los valores de la posición de la onda, que son del orden de [math]10^{306}[/math] , obviamente valores no razonables. Es por esto que tampoco podemos utilizar este método para la resolución de nuestro problema.

2.4 Cálculo de la energía

A continuación vamos a hallar la energía de nuestra ecuación de ondas, que viene definida según la ecuación: \begin{equation} E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx \end{equation}

Para ello hemos utilizado la resolución de la ecuación de ondas por el método de diferencias finitas, añadiendo lo necesario para representar la gráfica de la energía.

clear all, clf
%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];
%Como el desplazamiento es nulo, sabemos que la velocidad en esos puntos
%también será nula
V=[UA;V;UB];
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
    %Aproximamos las ux en los extremos de la cuerda
    Ux(1)=U(2,l)/k;
    Ux(end)=-U(end-1,l)/k;
    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.
plot(t,E)
%Elegimos el intervalo de los ejes
axis([0,40,0,0.5]) 
xlabel('Tiempo (s)'); ylabel('Energía (J)');

El resultado de esto es la siguiente gráfica:

Gráfica de la posición de cada punto del cable.

De la gráfica resultante vemos que la energía es constante, está energía es la mecánica. Al principio como la única fuerza que hay es la tensión, así que la energía quue hay es la potencial. A medida que pasa el tiempo aumenta la energía cinética, y se complementa manteniéndose constante en todo tiempo.

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

Vamos a estudiar el movimiento del mismo cable, pero esta vez sumergido en un medio viscoso, como sería el caso de un cable sumergido en el mar. El problema que modeliza este comportamiento es el siguiente:

[math] \begin{cases} u_{tt}-u_{xx}+au_{t}=0, \; x∈[0,L], \; t∈[0,40], \\ u(0,t)=0, \; u(L,t)=0, \\ u(x,0)=\begin{cases} \frac{3x}{L} \\ \frac{3}{2}-\frac{3x}{2L} \end{cases}, \; u_{t}(x,0)=0\\ \end{cases} [/math]

donde [math] a [/math] es una constante que depende del amortiguamiento que produce el medio. Se puede ver que [math] a [/math] actúa sobre la velocidad del cable, lo que será importante para la interpretación posterior. Vamos a estudiar el comportamiento del cable para [math] a=0,1,4,10,100 [/math].

Al ser una ecuación diferente, cambia la ecuación matricial que se obtiene al plantear el método de líneas para el problema propuesto. La ecuación matricial diferencial ordinaria de segundo orden que se obtiene para el caso de amortiguamiento viscoso es:

[math] \begin{cases} U''=-aU'-KU+F\\ U(0)=u^{0}\\ U'(0)=0\\ \end{cases} [/math]

y pasando a un sistema de ecuaciones matriciales diferenciales ordinarias de primer orden:

[math] \begin{cases} U' = V \\ V' = -aV -KU + F \\ U(0) = u^{0} \\ V(0) = 0 \\ \end{cases} [/math]

Aplicamos ahora el método del trapecio a cada ecuación del sistema, obteniendo que:

[math] \begin{cases} U_{n+1} = U_{n} + \frac{h}{2}(V_{n} + V_{n+1}) \\ V_{n+1} = V_{n} + \frac{h}{2}(-aV_{n} -KU_{n} + F_{n} -aV_{n+1} - KU_{n+1} + F_{n+1}) \\ \end{cases} [/math]

y despejando cada variable:

[math] \begin{cases} V_{n+1} (I + \frac{h^2}{4}K + \frac{h}{2}I) = V_{n} + \frac{h}{2}(-aV_{n} -KU_{n} + F_{n} + F_{n+1}) - \frac{h}{2}K(U_{n} + \frac{h}{2}V_{n}) \\ U_{n+1} = U_{n} + \frac{h}{2}(V_{n} + V_{n+1}) \\ \end{cases} [/math]

Una vez realizado este proceso analítico, se aplica la resolución numérica con MatLab/OctaveUPM.

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(1)=U(2,l)/k;
        Ux(end)=-U(end-1,l)/k;
        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

En el gráfico adjunto se puede observar una representación del valor de la energía para cada valor de [math] a [/math]. Cuando [math] a=0 [/math], podemos ver que la energía se mantiene sensiblemente constante, tal y como ocurría en el apartado anterior. Sin embargo, cuando [math] a=1 [/math] la energía decrece rápidamente. Este efecto disminuye según va aumentando el valor de [math] a [/math], como se puede observar por ejemplo para el caso de [math] a=100 [/math]. ¿Por qué ocurre esto? De primera mano, tal vez podríamos pensar que cuanto mayor sea el valor de [math] a [/math], antes se disipara la energía. Sin embargo, esto no es así, ya que la energía es la suma de la energía cinética y de la energía potencial, tal y como se muestra en la ecuación del apartado anterior. El valor de [math] a [/math] solo afecta a la energía cinética, no a la potencial. Es por ello que cuando el valor de [math] a [/math] es bajo, el cable tiene mayor facilidad para moverse, transformándose la energía potencial inicial en cinética, adquiriendo por lo tanto una mayor velocidad, siendo esta velocidad la que produce esta energía cinética se disipe. Por otro lado, cuando el valor de [math] a [/math] es elevado, al cable le cuesta más moverse, tranformándose menos energía potencial en cinética, y tardando por ello más en disiparse la energía total. En el caso extremo de que el valor de [math] a [/math] tendiera a infinito, el cable no se movería, siendo toda su energía potencial, y manteniéndose constante a lo largo del tiempo.

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. Por tanto, nuestro problema queda: [math] \begin{cases} u_{tt}-u_{xx}=0, \; x∈[0,10], \; t∈[0,60],\\ u(0,t)=\sin(2\pi*F0*t), \; u(10,t)=0,\\ u(x,0)=0, \; u_{t}(x,0)=0.\\ \end{cases} [/math]

El código MatLab para resolver dicho problema es el siguiente:

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
F1=0*xx;
F2=F1;
F1(end)=F1(end)+ub/h^2;
F2(end)=F2(end)+ub/h^2;
%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
    F1(1)=sin(2*pi*F0*t(i))/h^2;
    F2(1)=sin(2*pi*F0*t(i+1))/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)+F1+F2)));
    V(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(V(:,i)+0.5*k*(-K*U(:,i)+F1+F2)-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=2*pi*F0*cos(2*pi*F0*t);
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(1)=U(2,l)/k;
    Ux(end)=-U(end-1,l)/k;
    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)
grid on
xlabel('Tiempo (s)'); ylabel('Energía (J/kg)');

Gráfica de la posición del cable en cada punto y en cada instante.

Gráfica de la energía del cable. F0=0.11 Hz.

Gráfica de la posición del cable en cada punto y en cada instante.

Gráfica de la energía del cable. F0=0.09 Hz.

Gráfica de la posición del cable en cada punto y en cada instante.

Gráfica de la energía del cable. F0=0.10 Hz.

Como se puede observar en las gráficas la energía del cable no es constante. Debido a que la estructura sufre unas vibraciones periódicas, esta fuerza se le transmite al cable en su extremo izquierdo en el primer instante. Como las vibraciones de dicha estructura varían a lo largo del tiempo (pues viene definida por la función [math]f(t)=\sin(2\pi*F0*t)[/math] ) , la energía transmitida va cambiando, y por eso se pueden observar las oscilaciones en la gráfica. Además también se observa, para una frecuencia de 0.11 y 0.09 Hz, que la energía va aumentando durante aproximadamente los primeros 40 segundos. Esto es debido a que en la gráfica se está obteniendo la energía del cable entero y, a medida que pasa el tiempo, la energía del extremo izquierdo se va trasmitiendo a su vez a lo largo del cable. A partir de los 40 primeros segundos, se podría decir que la energía que había sido transmitida por la vibración de la estructura al extremo izquierdo del cable al principio ha "llegado" al final del cable. Esto se traduce gráficamente a que la curva de la energía no sigue ascendiendo, si no que habría alcanzado su tope y a partir de ahí oscilaría en esos intervalos de energía durante un período corto de tiempo y después, la contribución de la fuerza, en lugar de sumarse, se contrarresta por el efecto de la vibración de la onda. Lo podemos representar gráficamente aumentando el intervalo de tiempo.

Gráfica de la posición del cable en cada punto y en cada instante.

Gráfica de la energía del cable en los primeros 200 segundos. F0=0.11 Hz.

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

Vamos a plantear una nueva situación, en la que el extremo izquierdo del cable ya no se encuentra fijo, sino que su posición y la pendiente del cable en ese punto vienen determinadas por una condición Robin, como se puede ver en la siguiente modelización del problema: [math] \begin{cases} u_{tt}-u_{xx}=0, \; x∈[0,L], \; t∈[0,40], \\ u_{x}(0,t)=b*u(0,t), \; u(10,t)=0, \\ u(x,0)=\begin{cases} \frac{3x}{L} \\ \frac{3}{2}-\frac{3x}{2L} \end{cases}, \; u_{t}(x,0)=0\\ \end{cases} [/math] donde [math] b=2,-2 [/math] es una constante. Para aplicar el método de líneas, es necesario realizar antes un planteamiento analítico como el que se realizó en la primera modelización de este trabajo, con la salvedad de que habrá que plantear una nueva ecuación para el primer nodo del bucle, que viene dada por: [math] u''_{0}(t)+\frac{(2hb+2)u_{0}(t)-2u_{1}(t)}{h^{2}}=0 [/math] lo que a su vez se traduce en un cambio de los valores [math] (1,1) [/math] y [math] (1,2) [/math] en la matriz [math] K [/math] dentro del programa de MatLab/OctaveUPM, cuyo código es el que sigue.

clear all, clf
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
b2=input('Inserte un valor para b: ');
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(1: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))-diag(ones(1,N-1),-1) -diag(ones(1,N-1),1));
%Condición Robin
K(1,1)=(1/h^2)*(2*h*b2+2);
K(1,2)=(1/h^2)*(-2);
%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=[U;UB];
%La velocidad en el extremo derecho es nula a lo largo del timepo
V=[V;UB];
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
    %Aproximamos las ux en los extremos de la cuerda
    Ux(1)=U(2,l)/k;
    Ux(end)=-U(end-1,l)/k;
    Ux=Ux';
    E(l)=trapz(x,V(:,l).^2)+trapz(x,Ux.^2); %Cálculo de la energía.
end
%Dibujamos el gráfico del cable
figure(1)
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');
%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 para b=2.

Gráfica de la energía del cable para b=2.

Gráfica de la posición del cable en cada punto y en cada instante para b=-2.

Gráfica de la energía del cable para b=-2.

Si observamos las gráficas, y especialmente las escalas de los ejes, se puede ver que existe una diferencia notable entre la energía para el caso de [math] b=2 [/math] y la del caso de [math] b=-2 [/math]. Esto se debe a que en el caso de b=2, la posición del extremo izquierdo del cable se mantiene entre unos valores próximos a 0. Aunque la energía no se mantiene constante para este caso, sí que se mantiene más o menos dentro de un mismo rango de valores.

Sin embargo, para el caso de [math] b=-2 [/math], una vez transcurrido un cierto tiempo desde el comienzo, se produce una subida brusca del extremo izquierdo del cable, debido a que el dispositivo que se encuentra en ese extremo (y que materializa la condición Robin), amplifica el efecto de la onda, produciéndose lo que podríamos llamar como un fallo, ya que al estar el otro extremo fijo, el cable debería deformarse notablemente, para poder cumplir las concidiones que se presentan.. Es por ello que también la energía se dispara, alcanzando valores de orden muy elevado.

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

Procedemos ahora a resolver el problema del cable vibrante del apartado 2, cuya interpretación podemos encontrar en el propio apartado, mediante el método de Fourier. Lo que buscamos ahora es comprobar que realmente la solución hallada en el ejercicio dos es la correcta mediante el uso de aproximaciones dadas por el método de Fourier.

clear all
a=0;b=10; %espacio
h=0.1;%En x--------Paso espacial=h
x=a:h:b;
t=0:0.1:40;
[Mx,Mt]=meshgrid(x,t);
u0=zeros(size(x)); %primera función valor inicial
  for i=1:length(x)
      if x(i)<10/3 
       u0(i)=(3*x(i))/10;
      else 
       u0(i)=1.5-1.5*x(i)/10;
      end
  end
u0t=0; %segunda función valor inicial
Q=input('Introduzca el número de autofunciones a tratar: ');
U=0;
 for j=1:Q
    p=sin(j*pi/10*x);
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    bj=1/(j*pi)*trapz(x,u0t.*p)/trapz(x,p.*p);
    U=U+(aj.*cos(j*pi*Mt/10)+bj.*cos(j*pi*Mt/10)).*sin(j*pi*Mx/10);
 end
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');
%Cadena de texto con las autofunciones tomadas, para el título de la gráfica.
Aut=sprintf('Autofunciones tomadas: %d',Q); 
title(Aut);

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.

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.

Como se puede apreciar en las diferentes gráficas, cada una corresponde a un número diferente de autofunciones utilizadas; conforme aumenta el número de autofunciones usadas se comprueba que la solución se aproxima más a la solución aproximada utilizando el método de líneas del ejercicio dos, confirmándonos que ambas aproximaciones son correctas.