tdofss_eig.m

MATLAB PROGRAMS FOR "VIBRATION SIMULATION USING MATLAB AND ANSYS" All the M-files which are liste

	echo off
%	tdofss_eig.m	eigenvalue problem solution for tdof undamped model
%	Solves for the eigenvalues and eigenvectors in the state space 
%	form of the tdof system.

	clear all;

%	define the values of masses, springs, dampers and Forces, setting c1=c2=0 for
%	undamped normal mode calculation

	m1 = 1;
	m2 = 1;
	m3 = 1;

	c1 = 0;
	c2 = 0;

	k1 = 1;
	k2 = 1;

%	define the system matrix, a

	a = [  0      1         0			 0		   0	  0
		-k1/m1 -c1/m1     k1/m1        c1/m1       0	  0	
		   0      0         0            1         0      0
		 k1/m2  c1/m2  -(k1+k2)/m2  -(c1+c2)/m2  k2/m2  c2/m2
		   0      0         0            0         0      1
		   0      0       k2/m3        c2/m3    -k2/m3 -c2/m3];

%	solve for the eigenvalues of the system matrix

	[xm,omega] = eig(a)

%	take the diagonal elements of the generalized eigenvalue matrix omega

	omegad = diag(omega);

%	convert to hz from radians/sec

	omegahz = omegad/(2*pi);

%	reorder the eigenvalues and eigenvectors from low to high frequency, 
%	keeping track of how the eigenvalues are ordered to reorder the 
%	eigenvectors to match, using indexhz

	[omegaorder,indexhz] = sort(abs(imag(omegad)))

	for  cnt = 1:length(omegad)

		omegao(cnt,1) = omegad(indexhz(cnt));	%  reorder eigenvalues

		xmo(:,cnt) = xm(:,indexhz(cnt));		%  reorder eigenvector columns

	end

	omegao

	xmo

%	check for any eigenvectors with imaginary position elements by checking
%	the first three position entries for each eigenvector (first, third and 
%	and fifth rows) and convert to real

	for  cnt = 1:length(omegad)

		if  (real(xmo(1,cnt)) & real(xmo(3,cnt)) & real(xmo(5,cnt))) == 0

			xmo(:,cnt) = i*(xmo(:,cnt));	%  convert whole column if imaginary

		else

		end

	end

	xmo

%	check for any eigenvectors with negative position elements for the first
%	displacement, if so change to positive to that eigenvectors for the same mode
%	are complex conjugates

	for  cnt = 1:length(omegad)

		if  real(xmo(1,cnt)) < 0

			xmo(:,cnt) = -1*(xmo(:,cnt));	%  convert whole column if negative

		else

		end

	end

	xmo

%	define the mass and stiffness matrices for normalization of eigenvectors
%	and for checking values in principal coordinates

	m = [m1 0 0
		 0 m2 0
		 0  0 m3];

	k = [k1		-k1		0
		-k1	   k1+k2  -k2
		 0		-k2    k2];	

%	define the position eigenvectors by taking the first, third and fifth
%	rows of the original six rows in xmo

	xmop1 = [xmo(1,:);  xmo(3,:);  xmo(5,:)]

%	define eigenvectors for the three degrees of freedom by taking 
%	the second, fourth and sixth columns

	xmop = [xmop1(:,2)  xmop1(:,4)  xmop1(:,6)]

%	normalize with respect to mass

	for  mode = 1:3

		xn(:,mode) = xmop(:,mode)/sqrt(xmop(:,mode)'*m*xmop(:,mode));

	end

	xn

%	calculate the normalized mass and stiffness matrices for checking

	mm = xn'*m*xn

	km = xn'*k*xn

%	check that the sqrt of diagonal elements of km are eigenvalues

	p = (diag(km)).^0.5;

	[p abs(imag(omegao(1:2:5,:)))]

%	rename the three eigenvalues for convenience in later calculations

	w1 = abs(imag(omegao(1)));

	w2 = abs(imag(omegao(3)));

	w3 = abs(imag(omegao(5)));