多自由度系统的wislson method求解方法

博大广阔

%wislson method
%the wislon method assume that the acceleration of the systemvaries
%linearly between two instants of time.A stability analysis of the wislon
%method shows that q>=1.37.thissection ,we shall consider the wilso method
%for a multidegree of freedom system.
%since:ddxt is assmue to vary linearly between ti and t(i+q),we can predict
%the value of ddxt at ti+h,0<h<q*h1
% ddxt(i+h)=ddx(i)+h/(q*h1)*erro; erro=ddxt(i+h)-ddx(i) （1） %依据梯形关系
%by intergrating to h ,and obtain dxt(i+h) and xt(i+h)  
%BY substituting h=q*h1;
%wei get :
% ddXiq= 6/(q*h1)^2*(Xiq-Xi)-6/(q*h1)*dXi-2*ddXi;
% dXiq=3/(q*h1)(Xiq-Xi)-(q*h1)*ddXi/2-2*dXi;
%at last we can get the eqution:
% a=(6/(q*h)^2)*M+(3/(q*h))*C+K;b=(6/(q*h)^2)*M+(3/(q*h))*C;c=(6/(q*h))*M+2*C;d=2*M+(q*h/2)*C;
% a*XXnew=Fold+h*(Fnew-Fold)+b*Xold+c*dXold+d*ddXold;      (2)振动方程的带入后的结果
%程序步骤
%1）from the know initial conditions X0 dx0 ddX0;
%2)select a suitable time step h and a suitable value of q(q usually taken as 1.4)
%3)find the new diaplacement vector XXnew by (2)equation
%4)calculate the acceleration ,velcocity and displacement vectors at next time :Xnew
%5)a=(6/(q^3*h^2));b=(6/(q^2*h));c=(1-3/q);
%ddXnew=(6/(q^3*h^2))*(XXnew-Xold)-(6/(q^2*h))*dXold+(1-3/q)*ddXold;
%dXnew=dXold+(h/2)*(ddXnew+ddXold);
%Xnew=Xold+h*Xold+((h^2)/6)*(ddXnew+2*ddXold);
function Newmark_method(ff,M,K,C,T)
clc
clear all
%需在内部设置函数————系统指定
M=[1 0;0 2];K=[6 -2 ;-2 8];C=[0 0;0 0];
t=0:h:Tmax;
FF=[1;1]*sin(10*t);
Xold=[0;0];dXold=[0;0];ddXold=[0;0];Fold=[0;0];
%%%%%%%%
accuracy=0.001;erro=10^4;h=0.02;Tmax=10;n=1;q=1.4;
a=(6/(q^2*h^2))*M+(3/(q*h))*C+K; b=(6/(q*h)^2)*M+(3/(q*h))*C; c=(6/(q*h))*M+2*C; d=2*M+(q*h/2)*C;
for t=0:h:Tmax
       Fnew=FF(n);                                  %新值计算
        %%%计算中间位移
   XXnew=inv(a)*(Fold+h*(Fnew-Fold)+b*Xold+c*dXold+d*ddXold);
   %%%计算新的加速度和速度及位移
   ddXnew=(6/(q^3*h^2))*(XXnew-Xold)-(6/(q^2*h))*dXold+(1-3/q)*ddXold;
    dXnew=dXold+(h/2)*(ddXnew+ddXold);
    Xnew=Xold+h*dXold+((h^2)/6)*(ddXnew+2*ddXold);
    %%%%更新加速度，位移及速度
    ddXold=ddXnew;    dXold=dXnew;   Xold=Xnew;Fold=Fnew;
  %%保存
   XX(n,:)=Xnew;dXX(n,:)=dXnew; TT(n,:)=t;n=n+1;
end
figure(1)
%plot(XX(1:end,1),dXX(1:end,1),'r');grid on
plot(TT,XX)
%figure(2)
%plot(XX(1:end,2),dXX(1:end,2),'r');grid on
end
！！！！！！如果有误敬请告知

VibrationMaster

Wilson Theta Method??