|
|
| (No se muestran 3 ediciones intermedias del mismo usuario) |
| Línea 2: |
Línea 2: |
| | | | |
| | =Cable tensado y ondas= | | =Cable tensado y ondas= |
| − |
| |
| − | Se nos presenta un cable utilizado para estructura civiles de longitud '''L=10m''', con una sección que se nos permite despreciar con respecto a su longitud. Además, las vibraciones que este sufre se pueden modelizar mediante la ecuación de ondas.
| |
| − | Expresaremos como el desplazamiento vertical la función <math>u(x,t)</math> que sufren los distintos puntos del cable en función de su posición a lo largo del mismo (discretizados en <math>x ∈ [0,L] </math>) y en distintos instantes de tiempo '''t'''.
| |
| − |
| |
| − | Inicialmente, el cable se encuentra en un medio de viscosidad nula. Por ello, en la ecuación ya mencionada no aparece ningún término que amortigüe el desplazamiento de los puntos del cable. El cable se encuentra sujeto por ambos extremos, lo que nos facilita las primeras condiciones de frontera (condiciones Dirichlet). Concretamente, establecemos que tanto el desplazamiento del extremo izquierdo (x=0) como el extremo derecho (x=10) son nulos.
| |
| − |
| |
| − | Se considera que el cable es perfectamente homogéneo en toda su longitud. Por tanto, no existen imperfecciones que alteren su comportamiento. Consecuentemente, la función a la que igualamos nuestra ecuación de ondas '''f(x)''' es nula. REVISAR: En el caso de que nuestro cable presentara alguna anomalía, ésta debería ser reflejada en la función anterior.
| |
| − |
| |
| − | Si las vibraciones que sufre el cable son lo suficientemente pequeñas, la función <math>u(x,t)</math> satisface la ecuación de ondas <math>u_tt - u_xx=f(x)</math>. El sistema completo que muestra la situación del cable es el siguiente:
| |
| − |
| |
| − | <math>
| |
| − | u_tt - u_xx=f(x)</math>
| |
| − | <math>u(0,t)=0</math>
| |
| − | <math>u(L,t)=0</math>
| |
| − | <math>u(x,0)=h_0(x)
| |
| − | </math>
| |
| − |
| |
| − | ==Euler explícito==
| |
| − | {{matlab|codigo=
| |
| − | clear all
| |
| − | a=0; b=10;
| |
| − | h=0.1;
| |
| − | x=a:h:b;
| |
| − | N=round(b-a)/h;
| |
| − | %Condiciones de Contorno
| |
| − | ua=0; ub=0;
| |
| − | xx=x(2:N); %Extraemos los elementos interiores del vector x
| |
| − | xx=xx';%Trasponemos el vector para poder utilizarlo en el sistema
| |
| − | U0=zeros(size(xx)); %Generamos un vector de ceros al que luego introduciremos los valores de la función según el punto del cable en el que nos encontremos
| |
| − | n=length(xx);
| |
| − | for i=1:n
| |
| − | if xx(i)<=10/3
| |
| − | U0(i)=(3/10)*xx(i); %Condición Inicial Posición: Función a ejecutar si x<=10/3
| |
| − | else
| |
| − | U0(i)=1.5-(3/20)*xx(i); %Condición Inicial Posición: Función a ejecutar si x>10/3
| |
| − | end
| |
| − | end
| |
| − | V0=0*xx; %Condición Inicial Velocidad
| |
| − | K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1)); %Matriz de coeficientes del sistema
| |
| − | F=0*xx;
| |
| − | %Añadimos el primer y último elemento de las condiciones iniciales
| |
| − | F(1)=F(1)+ua/h^2; %Extremo izquierdo (x=0)
| |
| − | F(end)=F(end)+ub/h^2; %Extremo derecho (x=10)
| |
| − | %Procedemos a la discretización del tiempo
| |
| − | t0=0; tM=40;
| |
| − | k=h;
| |
| − | t=t0:k:tM;
| |
| − | M=(tM-t0)/k;
| |
| − | %El sistema a resolver es W'=LW+G, W(0)=W(super)0
| |
| − | W0=[U0;V0]; %Valor Inicial Sistema de orden 1
| |
| − | G=[zeros(size(xx));F]; %Vector término independiente Sistema de orden 1
| |
| − | L=[zeros(size(K)),eye(size(K));-K,zeros(size(K))];
| |
| − | %RESOLUCIÓN DEL SISTEMA (Por Euler Explícito)
| |
| − | sol(:,1)=W0;
| |
| − | for j=1:M
| |
| − | sol(:,j+1)=sol(:,j)+k*(L*sol(:,j)+G);
| |
| − | end
| |
| − | sol=sol(1:N-1,:);
| |
| − | %Tomamos las N-1 primeras filas e incluimos las condiciones de contorno
| |
| − | UA=ua*ones(1,length(t));
| |
| − | UB=ub*ones(1,length(t));
| |
| − | sol=[UA;sol;UB];
| |
| − | %Hay que tomar la parte que afecta a los primeros L/3 del cable de sol1 y
| |
| − | %la parte que afecta al resto del cable de sol2
| |
| − | [Mt,Mx]=meshgrid(t,x);
| |
| − | set(gcf,'renderer','painters');
| |
| − | mesh(Mx,Mt,sol)
| |
| − | }}
| |
| − |
| |
| − | ==Método de Heun==
| |
| − | {{matlabcodigo=
| |
| − | clear all
| |
| − | a=0; b=10;
| |
| − | h=0.1;
| |
| − | x=a:h:b;
| |
| − | N=round(b-a)/h;
| |
| − | %Condiciones de Contorno
| |
| − | ua=0; ub=0;
| |
| − | xx=x(2:N); %Extraemos los elementos interiores del vector x
| |
| − | xx=xx';%Trasponemos el vector para poder utilizarlo en el sistema
| |
| − | U0=zeros(size(xx)); %Generamos un vector de ceros al que luego introduciremos los valores de la función según el punto del cable en el que nos encontremos
| |
| − | n=length(xx);
| |
| − | for i=1:n
| |
| − | if xx(i)<=10/3
| |
| − | U0(i)=(3/10)*xx(i); %Condición Inicial Posición: Función a ejecutar si x<=10/3
| |
| − | else
| |
| − | U0(i)=1.5-(3/20)*xx(i); %Condición Inicial Posición: Función a ejecutar si x>10/3
| |
| − | end
| |
| − | end
| |
| − | V0=0*xx; %Condición Inicial Velocidad
| |
| − | K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1)); %Matriz de coeficientes del sistema
| |
| − | F=0*xx;
| |
| − | %Añadimos el primer y último elemento de las condiciones iniciales
| |
| − | F(1)=F(1)+ua/h^2; %Extremo izquierdo (x=0)
| |
| − | F(end)=F(end)+ub/h^2; %Extremo derecho (x=10)
| |
| − | %Procedemos a la discretización del tiempo
| |
| − | t0=0; tM=40;
| |
| − | k=h;
| |
| − | t=t0:k:tM;
| |
| − | M=(tM-t0)/k;
| |
| − | %El sistema a resolver es W'=LW+G, W(0)=W(super)0
| |
| − | W0=[U0;V0]; %Valor Inicial Sistema de orden 1
| |
| − | G=[zeros(size(xx));F]; %Vector término independiente Sistema de orden 1
| |
| − | L=[zeros(size(K)),eye(size(K));-K,zeros(size(K))];
| |
| − | %RESOLUCIÓN DEL SISTEMA (Por Euler Explícito)
| |
| − | sol(:,1)=W0;
| |
| − | for j=1:M
| |
| − | K1=L*sol(:,j)+G;
| |
| − | K2=sol(:,j)+K1*k;
| |
| − | sol(:,j+1)=sol(:,j)+k/2*(K1+K2);
| |
| − | end
| |
| − | sol=sol(1:N-1,:);
| |
| − | %Tomamos las N-1 primeras filas e incluimos las condiciones de contorno
| |
| − | UA=ua*ones(1,length(t));
| |
| − | UB=ub*ones(1,length(t));
| |
| − | sol=[UA;sol;UB];
| |
| − | %Hay que tomar la parte que afecta a los primeros L/3 del cable de sol1 y
| |
| − | %la parte que afecta al resto del cable de sol2
| |
| − | [Mt,Mx]=meshgrid(t,x);
| |
| − | set(gcf,'renderer','painters');
| |
| − | mesh(Mx,Mt,sol)
| |
| − | }}
| |
