ctsim.m

很多matlab的源代码

function y = ctsim(b,a,u,t,ic,ty)
% CTSIM	 Numerical solution of CT differential equations
%
%	Y = CTSIM(N,D,X,t,IC,TY): Numerical solution of differential equations
%	N, D = num, den of the operational transfer function.
%	X = input as STRING function of t. 
%	If X ='i', y returns the IMPULSE RESPONSE. For X=0, use X='0*t' or []
%	t = regularly spaced time vector
%	IC = initial condition vector (DEFAULT: IC = 0)
%	TY = Numerical Integration routine which can be one of the following: 
%	'rk1','rk2','rk4' (Runge-Kutta order 1,2,4), 'sim'(Simpson) (Def 'rk2')
%
%	Y returns the soultion at the time values in t
%	NOTE: All derivatives of the input assumed to equal 0 at t=0-.
%
%	CTSIM (with no input arguments) invokes the following example:
%
%	% Simulate and plot the response of [s*s+4s+3]y=[s+2]x to sin(10t) 
%	% for 6s. Assume initial conditions y(0)=1, y'(0)=2:  
%	  >>tr=0:0.05:6;   %Generate time array
%	  >>yr=ctsim([1 2],[1 4 3],'sin(10*t)',tr,[1;2]);
%	  >>plot(tr,yr),grid


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if nargin==0,help ctsim,disp('Strike a key to see results of the example')
pause,tr=0:0.05:6;yr=ctsim([1 2],[1 4 3],'sin(10*t)',tr,[1;2]);
plot(tr,yr),grid,return,end

[a,b,c,d]=ode2ss(b,a,'s');
if u=='i',x=ss2mat(a,b);else,x=ss2mat(a,b,u);end,
[mm,nn]=size(x);
if nargin<5,ty='rk2';ic=zeros(mm,1);end
if nargin==5,if isstr(ic),ty=ic;ic=zeros(mm,1);else,ty='rk2';end,end
if u=='i',ic=b;
else
r1=c;r=r1;for j=2:mm,r1=r1*a;r=[r;r1];end,ic=inv(r)*ic(:); %Transform IC
%DELETE THE ABOVE LINE TO MAKE RESULTS IDENTICAL TO LSIM
end
l=length(t);t0=t(1);tf=t(l);h=t(2)-t(1);
ty=ty(1:3);if ty=='rk1',[t,z]=oderk1(x,t0,tf,ic,h);
elseif ty=='rk2',[t,z]=oderk2(x,t0,tf,ic,h);
elseif ty=='rk4',[t,z]=oderk4(x,t0,tf,ic,h);
elseif ty=='sim',[t,z]=odesimp(x,t0,tf,ic,h);
else,error('Ungecognized integration routine');return
end
y=z(1:l,:)*c(:);
if u=='i',return,else,u=eval(u);u=u(1:l);y=y+u(:)*d.';end