-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Metodo iteración del subespacio listo
- Loading branch information
Showing
7 changed files
with
194 additions
and
23 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,37 @@ | ||
| function [V] = GSM(V, M, V_b) | ||
| % GSM: Gram - Shmidt Modificado | ||
| % Esta funcion deja normalizados y ortogonales entre si (en base a la | ||
| % matriz M) los vectores de entrada, mediante el metodo de Gram-Schmidt | ||
| % modificado, sistema mas estabe que Gram-Schmidt clasico dado que itera | ||
| % por bloques. | ||
| % | ||
| % [V] = GSM(V,M,V_b) | ||
| % | ||
| % Input: | ||
| % V Matriz cuyas columnas corresponden a los vectores a ortonormalizar | ||
| % M Matriz de base para ortonormalizacion | ||
| % V_b Matriz cuyas columnas ya estan ortonormalizadas en M | ||
|
|
||
| % En una primera etapa, si se tiene la matriz V_b se ortogonalizan todos | ||
| % los vectores en V a los vectores en esta matriz | ||
| if nargin == 3 | ||
| [~, n1] = size(V_b); | ||
| for k = 1:n1 | ||
| alpha = V_b(:, k)' * M * V; | ||
| V = V - V_b(:, k) * alpha; | ||
| end % for k | ||
| end | ||
| [~, n1] = size(V); | ||
|
|
||
| % El metodo de Gram-Schmidt modificado consiste en, asumiendo en cada | ||
| % iteracion que el vector considerado ya es ortogonal a los procesados | ||
| % anteriormente, se normaliza el vector y ortogonalizan a este con todos | ||
| % los vectores posteriores | ||
| for k = 1:n1 - 1 | ||
| V(:, k) = V(:, k) / (V(:, k)' * M * V(:, k))^0.5; | ||
| alpha = V(:, k)' * M * V(:, (k + 1):n1); | ||
| V(:, (k + 1):n1) = V(:, (k + 1):n1) - V(:, k) * alpha; | ||
| end % for k | ||
| V(:, n1) = V(:, n1) / (V(:, n1)' * M * V(:, n1))^0.5; | ||
|
|
||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,86 @@ | ||
| function [v, w] = calculoEigItSubespacio(M, K, nModos, tol) | ||
| % calculoEigItSubespacio: Calcula los valores y vectores propios del sistema | ||
| % usando el algoritmo de iteracion del subespacio | ||
| % | ||
| % [w,v] = calculoEigItSubespacio(M,K,nModos,tol) | ||
| % | ||
| % Parametros: | ||
| % M, K Matriz masa y rigidez | ||
| % nModos Numero de modos del analisis | ||
| % tol Tolerancia del calculo | ||
| % | ||
| % Salida: | ||
| % w Matriz vectores propios | ||
| % v Valores propios | ||
|
|
||
| if nargin < 4 || isempty(tol) | ||
| tol = 0.001; | ||
| end | ||
|
|
||
| % Calculos iniciales | ||
| ngdl = length(K); | ||
|
|
||
| % Calcula matriz de rigidez triangular LU | ||
| [L, U] = lu(K); | ||
| U = U'; | ||
| L = L'; | ||
|
|
||
| % Genera vector v a partir de numeros aleatorios | ||
| Vsize = min(6, ngdl); | ||
| X = (rand(ngdl, Vsize)); | ||
|
|
||
| phi = X; % Matriz que contendra los vectores finales | ||
| wfinal = zeros(nModos, 1); % Contendra las frecuencias | ||
|
|
||
| % Primera estimacion de la frecuencia | ||
| wo = sqrt(X(:, 1)'*K*X(:, 1)) / (X(:, 1)' * X(:, 1)); | ||
|
|
||
| % Genera los vectores propios iterando | ||
| for i = 1:nModos | ||
|
|
||
| err = 1000; % Reinicializa error | ||
| temp = []; %#ok<NASGU> | ||
|
|
||
| while err > tol | ||
|
|
||
| % Resuelve el sistema de ecuaciones en dos etapas | ||
| % L*D*y=b y | ||
| % L'*u=y | ||
| % donde L*D=U' (U de la descomposicion LU) | ||
| y1 = U \ (M * X); | ||
| X = L \ y1; | ||
|
|
||
| % Genera bloque de vectores u_s, rigidez y masa ortogonal | ||
| M_b = X' * M * X; | ||
| K_b = X' * K * X; | ||
| [Z] = eig(K_b, M_b); | ||
| X = X * Z; % Vectores ortogonales | ||
|
|
||
| % Calcula Gram-Schmidt para hacer X ortogonal a todos los | ||
| % demas vectores | ||
| if i ~= 1 | ||
| X = GSM(X, M, phi(:, 1:i-1)); | ||
| end | ||
| X = modonorm(X, M); % Normalizacion modal | ||
|
|
||
| % Verifio tolerancia con valores propios | ||
| Xo = X(:, 1); | ||
| w = sqrt(Xo'*K*Xo) / (Xo' * M * Xo); | ||
| err = abs((wo - w)/wo); | ||
| if err <= tol | ||
| wfinal(i) = w; | ||
| phi(:, i) = X(:, 1); | ||
| X = [X(:, 2:end), (rand(ngdl, 1))]; | ||
| else | ||
| wo = w; | ||
| end | ||
|
|
||
| end | ||
|
|
||
| end % for i | ||
|
|
||
| % Guarda valores | ||
| v = phi; | ||
| w = wfinal; | ||
|
|
||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,7 +1,44 @@ | ||
| function w = modonorm(v) | ||
| % modonorm: Normaliza un modo | ||
| function [phi] = modonorm(phi, M, opt) | ||
| % modonorm: Normaliza formas modales | ||
| % | ||
| % w = modonorm(phi,M,opt) | ||
|
|
||
| vn = norm(v); | ||
| w = v ./ vn; | ||
| if nargin < 2 | ||
| opt = 1; | ||
| end | ||
|
|
||
| if nargin < 3 || isempty(opt) | ||
| opt = 2; | ||
| end | ||
|
|
||
| % Para caso general | ||
| if opt == 0 | ||
| end | ||
|
|
||
| % Normaliza con primera linea | ||
| if opt == 1 | ||
| aux = (phi(1, :)); | ||
| if min(abs(aux)) <= eps | ||
| opt = 3; | ||
| warning('No se puede dividir debido a un valor nulo, se usa opt=3'); | ||
| else | ||
| phi = phi * diag(1./aux); | ||
| end | ||
| end | ||
|
|
||
| if opt == 2 | ||
| Mn = sqrt(diag(conj(phi)'*M*phi)); | ||
| phi = phi * diag(Mn^-1); | ||
| end | ||
|
|
||
| if opt == 3 | ||
| aux = max(phi); | ||
| phi = phi * diag(1./aux); | ||
| end | ||
|
|
||
| % Simplemente normaliza por el valor absoluto | ||
| if opt == 4 | ||
| phi = phi ./ norm(phi); | ||
| end | ||
|
|
||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters