finite element method using matlab 錯誤
本帖最后由 khyiu0415 于 2013-1-26 19:46 编辑請問大家有看過這本書嗎,當中的m-file'fediresp'是否有錯誤嗎?請指教
function =fediresp(M,K,F,u,t,C,q0,dq0,a,b);
%----------------------------------------------------------------------------------
%Purpose:
% The function subroutine fediresp.m calulates impulse response
% for a damped structural system using modal analysis. It uses modal
% coordinate equations to evaluate modal responses anaytically, then
% convert modal coordinates into physical responses
%
%Synopsis:
% =fediresp(M,K,F,u,t,C,q0,dq0,a,b)
%
%Variable Description:
% Input parameters - M, K - Mass and stiffness matrices
% F - Input or forcing influence matrix
% t - Time of evaluation
% u - Index for excitation
% C - Output matrix
% q0, dq0 - Initial conditions
% a, b - Parameters for proportional damping =a+b
% Outpur parameters - eta - modal coordinate response
% y - physical coordinate response
%---------------------------------------------------------------------------------
disp('')
disp('Please wait!! - The job is being performed.')
%---------------------------------------------------------------
%Solve the eigenvalue problem and normalize the eigenvectors
% --------------------------------------------------------------
=size(M);=size(F);
nstep=size(t');
=eig(K,M);
=sort(diag(D));% Sort the eigenvaules and eigenvectors
V=V(:,k);
Factor=diag(V'*M*V);
Vnorm=V*inv(sqrt(diag(Factor))); %Eigenvectors are normailzed
omega=diag(sqrt(Vnorm'*K*Vnorm));% Natural frequencies
Fnorm=Vnorm'*F;
%----------------------------------------------------------------------------
%Compute modal damping matrix from the proportional damping matrix
%----------------------------------------------------------------------------
Modamp=Vnorm'*(a*M+b*K)*Vnorm;
zeta=diag((1/2)*Modamp*inv(diag(omega)))
if (max(zeta) >= 1),
disp('Warning - Your maximu damping ratio is grater than or equal to 1')
disp('You have to reselect a and b ')
else
pause
disp('If you want to continue, type return key')
end
%----------------------------------------------------------------------
%Find out impulse response of each modal coordinate analytically
%-------------------------------------------------------------------
eta0=Vnorm'*M*q0; % Initial conditions for modal coordinates
deta0=Vnorm'*M*dq0; % - both displacement and velocity
eta=zeros(nstep,n);
for i=1:n %Responses are obtained for n modes
omegad=omega(i)*sqrt(1-zeta(i)^2);
phase=omegad*t;
tcons=zeta(i)*omega(i)*t;
eta(:,i)=exp(-tcons)'.*(eta0(i)*(cos(phase')+zeta(i)/sqrt(1-zeta(i)^2)*...
sin(phase'))+deta0(i)*sin(phase')/omegad+sin(phase')*Fnorm(i,u)/omegad);
end
%-----------------------------------------------------------------------
%Convert modal coordinate responses to physical coordinate responses
%-----------------------------------------------------------------------
yim=C*Vnorm*eta';
%-----------------------------------------------------------------------
具体错误请给出。
页:
[1]