ndof_numzeros.m

振动仿真 vibration simulation using matlab and ansys!

	echo off
%	ndof_numzeros.m	 number of zeros in an arbitrarily sized series spring/mass
%	model transfer function.  Calculates and plots poles/zeros and
%	transfer functions for user selected input/output locations on a
%	"n" dof series spring/mass model.  Shows that poles of "constrained"
%	structures to left and right of input/output degrees of freedom are 
%	the zeros of the unconstrained structure.  

	clf;

	subplot(1,1,1);

	legend off;

	n = input('Input number of masses (dof) in model, default = 10 ...  ');

	if (isempty(n))
   		n = 10;
   	else
	end

	i = input('Input dof number for left side, 1:n, default = 1 ...  ');

	if (isempty(i))
   		i = 1;
   	else
	end

	j = input(['Input dof number for right side, 1:n, default = ',num2str(n),' ...  ']);

	if (isempty(j))
   		j = n;
   	else
	end

%	define random values between 1 and 2 for m and k

   	m = 1+rand(1,n);

	k = 1+rand(1,n+1);

%	build the mass and stiffness matrices

	M = diag(m);

	K = zeros(n);

	for  row = 1:n

		K(row,row) = k(row) + k(row+1)

	end

	for  row = 1:n-1

		colr = row + 1;

		K(row,colr) = -k(row+1);

	end

	for  row = 2:n

		coll = row - 1;

		K(row,coll) = -k(row);

	end

%	now form the "a" matrix

	a = -inv(M)*K

%	now form b vector

	b = zeros(n,1);

	b(i) = 1;

	b = inv(M)*b

%	now change to state-space form

	ass = zeros(2*n);

%	fill in the ones on odd rows

	for  row = 1:2:2*n
	
		ass(row,row+1) = 1;
		
	end
	
%	fill in the zeros on the even rows

	for  row = 2:2:2*n

		rowa = row/2;

		for  col = 1:2:2*n

			cola = ceil(col/2);

			ass(row,col) = a(rowa,cola);

		end

	end

	ass

%	expand the b vector for state space

	bss = zeros(2*n,1);

	for  row = 2:2:2*n

		rowb = row/2;

		bss(row) = b(rowb);

	end

	bss
	
%	define c, output vector, for transfer function dof(j)/dof(i)

	css = zeros(1,2*n);

	css(2*j-1) = 1;

	css

%	plot transfer function

%	define a vector of frequencies to use, radians/sec

	w = logspace(-1,1,800);

	[mag,phs] = bode(ass,bss,css,[0],1,w);

	magdb = 20*log10(mag);

	semilogx(w,magdb,'k-')
	title(['transfer function, ',num2str(n),' dof, input at ',num2str(i),', output at ',num2str(j)]);
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

%	plot poles/zeros

	[p,z] = pzmap(ass,bss,css,[0]);
	maxp = max(abs(p));

	plot(p,'k*')
	title('poles/zeros of system')
	hold on

%	check to see if there are any zeros

	if  ((i >= 2) | (j <= n-1))

	plot(z,'ko')

	else
	end

	axis([-maxp maxp -maxp maxp]);
	axis('square')
	grid on
	disp('execution paused to display figure, "enter" to continue'); pause
	hold off
	
	subplot(1,3,1)
	plot(p,'k*')
	title('poles/zeros of system')
	hold on

%	check to see if there are any zeros

	if  ((i >= 2) | (j <= n-1))

	plot(z,'ko')

	else
	end

	axis([-maxp maxp -maxp maxp]);
	axis('square')
	grid on
	hold off
	
%	now define the constrained system to the left of the "i" dof as long as 
%	"i" is greater than 1

%	check to see if there are any zeros

	if  i >= 2

		assl = ass(1:2*(i-1),1:2*(i-1));

		[evecl,evall] = eig(assl);

		subplot(1,3,2)
		plot(diag(evall),'k*')
		title('poles of lhs')
		axis([-maxp maxp -maxp maxp]);
		axis('square')
		grid on

	else

	end
	
%	check to see if there are any zeros

	if  j <= n-1

		assr = ass(2*(j+1)-1:2*n,2*(j+1)-1:2*n);

		[evecr,evalr] = eig(assr);

		subplot(1,3,3)
		plot(diag(evalr),'k*')
		title('poles of rhs')
		axis([-maxp maxp -maxp maxp]);
		axis('square')
		grid on

	else

	end