Práctica 11 - Informe

Computación II.
Práctica 11 Xulong Ye
Práctica 11. Informe

1. Planteamiento del problema
El problema presentado de los bloques de viviendas se puede a partir de unas
aproximaciones, modelizarlo como uno de osciladores acoplados. De esta forma si
aplicamos las 2º Ley de Newton para cada bloque y teniendo en cuenta que la única
fuerza que actúa es la elástica, tenemos que:
(1)
Donde m es la masa, x la posición y k la constante elástica. Los subíndices i indica a la
partícula a la que se refiere.
Teniendo en cuenta que las masas y las constantes de los mulles son iguales, de donde
m=4000 kg y k=5000 N/m, la ecuación (1) se puede escribir como:
−1
(2) donde despejando M nos queda X’’ + Kr X = 0 (3), Kr= M K
De aquí puedo hallar el valor de la matriz Kr que es:

Nota: Téngase en cuenta de la ec.(1) que el
primer bloque y el último bloque en la ec.
de movimiento (3) solo le afecta su propia
posición y la del bloque superior en el
primero y el inferior en el segundo
Ahora como nos interesa resolver este sistema de osciladores acoplados, se puede
resolver calculando la matriz Kr en una base en la que sea diagonal. De esta forma, si
encontramos los autovalores de Kr tendremos los coeficientes de X en esa base, por el
cual se encuentra el sistema desacoplado y tendremos la Ec. del MAS (Fácil de
resolver). En ella la nueva matriz Kr, tomará los valores de la frecuencia angular al
cuadrado, y si queremos pasar de esta base ortogonal a la base original solo necesitamos
multiplicárselo a la matriz de paso U, la de autovectores. Como resultado nos queda que
la ec. de movimiento de cada bloque es una combinación lineal de MAS con distintas
frecuencias correspondiente a los autovalores, y cuyo factor viene dado por los
autovectores.
2. Código C++ para diagonalizar
Primero de todo, escribimos las librerías que vamos a usar:
#include <iostream> Librería de C++
#include <cmath> Para operaciones trigonométricas.
#include <fstream> Para crear un flujo de entrada
#include<iomanip> Para presentar el resultado con los decimales que nos interese
Computación II. Práctica 11 Xulong Ye
Luego declaramos nuestra variables n, dimensión de la matriz, la tolerancia, tol=1E-10

y la función con argumentos de entrada la matriz a diagonalizar A[n][n], su dimensión y
la tolerancia.
Posteriormente, vamos a escribir nuestra función main, leyendo la matriz a diagonalizar
a partir de un fichero externo (“kr.dat”) y con la matriz ya determinada, llamamos a la
función diag_Jacobi.
2.1. Explicación de la función diag_Jacobi
El algoritmo que vamos a emplear es el de la rotación de Jacobi, siguiendo las
instrucciones dadas.
Paso 1: Construir una matriz U(nxn) = In (matriz identidad)
for(int i=0;i<n;i++){
for (int j=0;j<n;j++){
if(i==j){
U[i][i]=1;
}
Paso 2: Declaramos un contador, las coordenadas del elemento con mayor valor
absoluto y su valor.
int contador_rotacion=0, imax, jmax; //Contador de rotaciones
double max_v;
Paso 3: Abrimos bucle do while, y hacemos contador de rotación aumente en una
unidad cada vez que se repita el bucle
do{
contador_rotacion=contador_rotacion+1;
Paso 4: Buscamos el elemento máximo por encima de la diagonal. Para ello recorremos
las filas desde el primero al último y las columnas desde la fila+1 hasta el último.
Sería equivalente mirar los elementos por debajo de la diagonal puesto que es simétrica
y ahora sería empezando por columnas y para las filas se empezaría desde la
columna+1. Guardamos las coordenadas del elemento con mayor valor absoluto y su
valor.
for (int ib=0; ib < n ;ib++) {
for (int jb=ib+1; jb < n; jb++)
{ if (fabs(A[ib][jb]) > max_v){
max_v=fabs(A[ib][jb]);
imax=ib;
jmax=jb; } } }
Paso 5: Hallamos theta de la fórmula derivada en clase y llamamos al seno y al coseno s
y c, respectivamente para facilitar la escritura:
(4)
double theta;
if (fabs(A[imax][imax]-A[jmax][jmax]) < tol*fabs(A[imax][imax]))
{ theta = atan(1.0);} // pi/4
else
{ theta = 0.5*atan(2.0*A[imax][jmax]/(A[imax][imax]-A[jmax][jmax]));}
double s=sin(theta);
double c=cos(theta);
Paso 6: Reemplazar la matriz A por RT AR. Para ello: a) desde k=0 hasta k=n-1 (k≠i,j),
calcular los elementos Aki (=Aik ) y Akj(=Ajk) mediante las ec. (2) y (3) ; b) calcular
los elementos Aii y Ajj mediante las (4 y (5); y finalmente, c) imponer que Aij =Aji=0
La matriz , U va rotando al igual que la matriz A
for(int k=0;k<n;k++){
double uki=U[k][imax]; //Daría errores si uso directamente U[k][imax]
U[k][imax]=uki*c+U[k][jmax]*s;
U[k][jmax]=-uki*s+U[k][jmax]*c;
if(k != imax && k != jmax){
A[k][imax]=A[imax][k]*c+A[jmax][k]*s;
A[k][jmax]=-A[imax][k]*s+A[jmax][k]*c;
A[imax][k]=A[k][imax];
A[jmax][k]=A[k][jmax];
}
}
double aii=A[imax][imax]; //Aquí igual, no se usa directamente A[imax][imax]
A[imax][imax]=aii*c*c+A[jmax][jmax]*s*s+2*A[imax][jmax]*s*c;
A[jmax][jmax]=aii*s*s+A[jmax][jmax]*c*c-2*A[imax][jmax]*s*c;
A[imax][jmax]=0;
A[jmax][imax]=0;
El motivo por el que hemos creado ‘uki’ y ‘aii’ es que U[k][imax] y A[imax][imax] lo
estamos cambiando, por lo que si en la siguiente línea lo volvemos a usar, se usaría el
nuevo valor. Sin embargo, debemos de trabajar con un mismo valor en cada iteración de
nuestro bucle do while.
Paso 7: El bucle acaba cuando el máximo valor de la diagonal es menor que la
tolerancia, es decir:
do{….}while(max_v>tol)
Paso 8: Mostramos en pantalla los autovalores, que son los elementos de la diagonal de
A y la matriz de autovectores que es la matriz U después de que solo quede elementos
en la diagonal para la matriz A. Mostramos tanto decimales como tolerancia hayamos
puesto
cout<< "EL numero de rotaciones realizadas es de"<< contador_rotacion<< endl;
cout<< "Los autovalores son: "<<endl;
cout <<fixed<<setprecision(10)<< A[i][i];
cout<< endl;
}
cout<< "La matriz de autovectores es: "<< endl;
for(int k=0;k<n;k++){
cout<<fixed<<setprecision(10)<<U[i][k]<< " ";
}
cout<< endl;
}}
3. Comprobación en Matlab
Para obtener la matriz de autovectores y autovalores en Matlab, es bastante sencillo,
solo necesitamos usar la función eig() de Matlab y nos devuelve la matriz de
autovalores y autovectores.
Imponiendo que las amplitudes de la base ortogonal son A=[15,4.5,7,6]m y

suponiendo una solución sinusoidal que varía con un ángulo phi0 que va desde 0 a pi/2
equiespaciados. Si uso la matriz de cambio de base U para obtener el resultado en la

base original, podemos llegar a esta gráfica.
Código Matlab:
clear all
close all
clc
%%Cálculo de los autovalores y autovectores de una matriz en
Matlab
kr=load('Kr.dat');
[V J]=eig(kr); %%De donde obtengo que V es la matriz de
autovectores y J la de autovalores
%% Suponiendo una condiciones iniciales que donde phi0=0 y donde
las amplitudes de cada una es igual
for i=1:length(J(1,:))
omega(i)=sqrt(J(i,i));
end
A=[15,4.5,7,6]; % Amplitud en metro imponiendo unas condiciones
iniciales
t=[linspace(0,100,10000);linspace(0,100,10000);linspace(0,100,10
000);linspace(0,100,10000)];
phi0=linspace(0,pi/2,4)'; %Condicion que imponemos
xp=A'.*sin(omega'.*t+phi0);
x=V*xp;
plot(t(1,:),x(1,:),t(2,:),x(2,:),t(3,:),x(3,:),t(4,:),x(4,:))
legend('Masa 1', 'Masa 2', 'Masa 3', 'Masa 4')
title('Posición de las distintas masas respecto del tiempo')
xlabel('Tiempo (s)')
ylabel('Posición (m)')
grid on

Práctica 11 - Informe

Cargado por

Copyright:

Formatos disponibles

Práctica 11 - Informe

Cargado por

Información del documento

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Práctica 11 - Informe

Cargado por

Copyright:

Formatos disponibles

Computación II.

Práctica 11. Informe

De aquí puedo hallar el valor de la matriz Kr que es:

Luego declaramos nuestra variables n, dimensión de la matriz, la tolerancia, tol=1E-10

Imponiendo que las amplitudes de la base ortogonal son A=[15,4.5,7,6]m y

equiespaciados. Si uso la matriz de cambio de base U para obtener el resultado en la

También podría gustarte