tdofxfer.m

振动仿真 vibration simulation using matlab and ansys!

	z22tf = tf(z22num,den);

	z11tf_lo = tf(z11num_lo,z11den_lo);

	z11tf_hi = tf(z11num_hi,z11den_hi);

	z21tf_lo = tf(z21num_lo,z21den_lo);

	z21tf_hi = tf(z21num_hi,z21den_hi);

	z31tf_lo = tf(z31num_lo,z31den_lo);

	z31tf_hi = tf(z31num_hi,z31den_hi);

	z22tf_lo = tf(z22num_lo,z22den_lo);

	z22tf_hi = tf(z22num_hi,z22den_hi);

%	use the bode command with left hand magnitude and phase vector arguments
%	to provide values for further analysis/plotting

	[z11mag,z11phs] = bode(z11tf,w);	

	[z21mag,z21phs] = bode(z21tf,w);	

	[z31mag,z31phs] = bode(z31tf,w);	

	[z22mag,z22phs] = bode(z22tf,w);	

	[z11maglo,z11phslo] = bode(z11tf_lo,w);

	[z21maglo,z21phslo] = bode(z21tf_lo,w);

	[z31maglo,z31phslo] = bode(z31tf_lo,w);

	[z22maglo,z22phslo] = bode(z22tf_lo,w);

	[z11maghi,z11phshi] = bode(z11tf_hi,w);

	[z21maghi,z21phshi] = bode(z21tf_hi,w);

	[z31maghi,z31phshi] = bode(z31tf_hi,w);

	[z22maghi,z22phshi] = bode(z22tf_hi,w);

%	calculate the magnitude in decibels, db

	z11magdb = 20*log10(z11mag);

	z21magdb = 20*log10(z21mag);

	z31magdb = 20*log10(z31mag);

	z22magdb = 20*log10(z22mag);

	z11maglodb = 20*log10(z11maglo);
	
	z21maglodb = 20*log10(z21maglo);
	
	z31maglodb = 20*log10(z31maglo);
	
	z22maglodb = 20*log10(z22maglo);
	
	z11maghidb = 20*log10(z11maghi);
	
	z21maghidb = 20*log10(z21maghi);
	
	z31maghidb = 20*log10(z31maghi);
	
	z22maghidb = 20*log10(z22maghi);
	
%	plot the four frequency responses using log magnitude

	loglog(w,z11mag(1,:),w,z21mag(1,:),w,z31mag(1,:),w,z22mag(1,:))
	title('xfer function form, Bode z11, z21, z31 and z22 magnitude')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('magnitude')
	grid

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

%	plot the four frequency responses using db

	semilogx(w,z11magdb(1,:),w,z21magdb(1,:),w,z31magdb(1,:),w,z22magdb(1,:))
	title('xfer function form, Bode z11, z21, z31 and z22 db magnitude')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	grid

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

%	plot phase

	semilogx(w,z11phs(1,:),w,z21phs(1,:),w,z31phs(1,:),w,z22phs(1,:))
	title('xfer function form, Bode z11, z21, z31 and z22 phase')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	grid

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

	legend off;

%	plot the four frequency responses separately with overlaid
%   asymptotes, in a 2x2 subplot form

	subplot(2,2,1)
	semilogx(w,z11magdb(1,:),'k',w,z11maglodb(1,:),'k--',w,z11maghidb(1,:),'k--')
	title('xfer function form, Bode z11, z33 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,2)
	semilogx(w,z21magdb(1,:),'k',w,z21maglodb(1,:),'k--',w,z21maghidb(1,:),'k--')
	title('xfer function form, Bode z21, z12, z23, z32 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,3)
	semilogx(w,z31magdb(1,:),'k',w,z31maglodb(1,:),'k--',w,z31maghidb(1,:),'k--')
	title('xfer function form, Bode z31, z13 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,4)
	semilogx(w,z22magdb(1,:),'k',w,z22maglodb(1,:),'k--',w,z22maghidb(1,:),'k--')
	title('xfer function form, Bode z22 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

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

	subplot(2,2,1)
	semilogx(w,z11phs(1,:),'k',w(1),z11phslo(1),'ko',w(200),z11phshi(200),'ko')
	title('xfer function form, Bode z11, z33 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	axis([.1 10 -400 -150])
	grid

	subplot(2,2,2)
	semilogx(w,z21phs(1,:),'k',w(1),z21phslo(1),'ko',w(200),z21phshi(200)-360,'ko')
	title('xfer function form, Bode z21, z12, z23, z32 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	axis([.1 10 -400 -150])
	grid

	subplot(2,2,3)
	semilogx(w,z31phs(1,:),'k',w(1),z31phslo(1),'ko',w(200),z31phshi(200),'ko')
	title('xfer function form, Bode z31, z13 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	axis([.1 10 -400 -150])
	grid

	subplot(2,2,4)
	semilogx(w,z22phs(1,:),'k',w(1),z11phslo(1),'ko',w(200),z11phshi(200),'ko')
	title('xfer function form, Bode z22 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	axis([.1 10 -400 -150])
	grid

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

%	plot frequency response asymptotes for notes and for hand plotting

	subplot(2,2,1)
	semilogx(w,z11maglodb(1,:),w,z11maghidb(1,:),'--')
	title('Asymptotes - Bode z11, z33 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,2)
	semilogx(w,z21maglodb(1,:),w,z21maghidb(1,:),'--')
	title('Asymptotes - Bode z21, z12, z23, z32 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,3)
	semilogx(w,z31maglodb(1,:),w,z31maghidb(1,:),'--')
	title('Asymptotes - Bode z31, z13 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,4)
	semilogx(w,z11maglodb(1,:),w,z11maghidb(1,:),'--')
	title('Asymptotes - Bode z22 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	clear all;

	clf;

	legend off;

	subplot(1,1,1);

%	using Matlab's "bode" plotting capability, defining the transfer
%	functions in "zero/pole/gain" form by column vectors of poles and zeros 
%   and defining output vectors for magnitude and phase as well as a
%	defined range of radian frequencies

%	assign values for masses, damping, and stiffnesses

	m1 = 1;
	m2 = 1;
	m3 = 1;
	c1 = 0;
	c2 = 0;
	k1 = 1;
	k2 = 1;

	m = m1;
	k = k1;

%	define column vectors of poles and zeros from previous derivation
%
%	there are three ways to make a column vector:
%
%	1)	define a row vector and then transpose it:
%
%		p = [0 0 1*j -1*j sqrt(3*k/m)*j -sqrt(3*k/m)*j]';
%
%	2)	define a column vector by using semi-colons between elements:
%
%		p = [0; 0; 1*j; -1*j; sqrt(3*k/m)*j; -sqrt(3*k/m)*j];
%
%	3)	define a column vector directly:
%
%		p = [	0
%			 	0
%		  	   1*j
%		 	  -1*j
%		   sqrt(3*k/m)*j
%		  -sqrt(3*k/m)*j];

%	zeros for z1/f1; quartic so four zeros

	z11_1 = -sqrt((-3*k-sqrt(5)*k)/(2*m));
	z11_2 = sqrt((-3*k-sqrt(5)*k)/(2*m));

	z11_3 = -sqrt((-3*k+sqrt(5)*k)/(2*m));
	z11_4 = sqrt((-3*k+sqrt(5)*k)/(2*m));

%	zeros for z2/f1; quadratic so two zeros

	z21_1 = -sqrt(-k/m);
	z21_2 = sqrt(-k/m);

%	zeros for z3/f1; no zeros, so use empty brackets

	z31_1 = [];

%	zeros for z2/f2: quadratic so two zeros

	z22_1 = -sqrt(-k/m);
	z22_2 = sqrt(-k/m);

%

	z11 = [z11_1 z11_2 z11_3 z11_4]';

	z21 = [z21_1 z21_2]';

	z31 = z31_1;

	z22 = [z22_1 z22_2]';

	p = [0 0 1*j -1*j sqrt(3*k/m)*j -sqrt(3*k/m)*j]';

	gain = 1;

%	use the zpk command to define the four pole/zero/gain systems

	sys11pz = zpk(z11,p,gain);

	sys21pz = zpk(z21,p,gain);

	sys31pz = zpk(z31,p,gain);

	sys22pz = zpk(z22,p,gain);

%	Define a vector of frequencies to use, radians/sec.  The logspace command uses
%	the log10 value as limits, i.e. -1 is 10^-1 = 0.1 rad/sec, and 1 is
%	10^1 = 10 rad/sec.  The 200 defines 200 frequency points.

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

%	use the bode command with left hand magnitude and phase vector arguments
%	to provide values for further analysis/plotting

	[z11mag,z11phs] = bode(sys11pz,w);	

	[z21mag,z21phs] = bode(sys21pz,w);	

	[z31mag,z31phs] = bode(sys31pz,w);	

	[z22mag,z22phs] = bode(sys22pz,w);	

%	calculate the magnitude in decibels, db

	z11magdb = 20*log10(z11mag);

	z21magdb = 20*log10(z21mag);

	z31magdb = 20*log10(z31mag);

	z22magdb = 20*log10(z22mag);

%	plot the four frequency responses using log magnitude

	loglog(w,z11mag(1,:),w,z21mag(1,:),w,z31mag(1,:),w,z22mag(1,:))
	title('Zero/Pole/Gain, Bode z11, z21, z31 and z22 magnitude')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('magnitude')
	grid

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

%	plot the four frequency responses using db

	semilogx(w,z11magdb(1,:),w,z21magdb(1,:),w,z31magdb(1,:),w,z22magdb(1,:))
	title('Zero/Pole/Gain, Bode z11, z21, z31 and z22 db magnitude')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	grid

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

%	plot phase

	semilogx(w,z11phs(1,:),w,z21phs(1,:),w,z31phs(1,:),w,z22phs(1,:))
	title('Zero/Pole/Gain, Bode z11, z21, z31 and z22 phase')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	grid

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

	legend off;

%	plot the four frequency responses separately, in a 2x2 subplot form

	subplot(2,2,1)
	semilogx(w,z11magdb(1,:))
	title('Zero/Pole/Gain, Bode z11, z33 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,2)
	semilogx(w,z21magdb(1,:))
	title('Zero/Pole/Gain, Bode z21, z12, z23, z32 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,3)
	semilogx(w,z31magdb(1,:))
	title('Zero/Pole/Gain, Bode z31, z13 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,4)
	semilogx(w,z22magdb(1,:))
	title('Zero/Pole/Gain, Bode z22 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

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

	subplot(2,2,1)
	semilogx(w,z11phs(1,:))
	title('Zero/Pole/Gain, Bode z11, z33 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	axis([.1 10 -400 -150])
	grid

	subplot(2,2,2)
	semilogx(w,z21phs(1,:))
	title('Zero/Pole/Gain, Bode z21, z12, z23, z32 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	axis([.1 10 -400 -150])
	grid

	subplot(2,2,3)
	semilogx(w,z31phs(1,:))
	title('Zero/Pole/Gain, Bode z31, z13 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	axis([.1 10 -400 -150])
	grid

	subplot(2,2,4)
	semilogx(w,z22phs(1,:))
	title('Zero/Pole/Gain, Bode z22 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	axis([.1 10 -400 -150])
	grid

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