振动理论与应用 matlab 程序

clear global;clg; t=0; x=0; y=0; % This program solves the differential equation of the motion of the mass- % spring system using the Runge Kutta method % The force function must be in the file f.m m=input('Enter the value of the mass (kg) : ') c=input('Enter the value of the damping coefficient (N.s/m) : ') k=input('Enter the value of the spring constant (N/m) : ') x1=input('Enter the value of the initial position (m) : ') y1=input('Enter the value of the initial velocity (m/s) : ') 'Natural period' T=2*pi*(m/k)^0.5 dt=input('Enter the value of time increment in seconds ( < T as above ) : ') t1=input('Enter the value of the initial time in seconds (0) : ') tf=input('Enter the value of the final time in seconds : ') x(1)=x1; y(1)=y1; t(1)=t1; n=(tf-t(1))/dt; N=round(n); for I=1:N, T1=t(I); X1=x(I); Y1=y(I); F1=(1/m)*(f(T1)-k*X1-c*Y1); T2=t(I)+dt/2; X2=x(I)+Y1*dt/2; Y2=y(I)+F1*dt/2; F2=(1/m)*(f(T2)-k*X2-c*Y2); T3=t(I)+dt/2; X3=x(I)+Y2*dt/2; Y3=y(I)+F2*dt/2; F3=(1/m)*(f(T3)-k*X3-c*Y3); T4=t(I)+dt; X4=x(I)+Y3*dt; Y4=y(I)+F3*dt; F4=(1/m)*(f(T4)-k*X4-c*Y4); x(I+1)=x(I)+dt/6*(Y1+2*Y2+2*Y3+Y4); y(I+1)=y(I)+dt/6*(F1+2*F2+2*F3+F4); t(I+1)=dt*I; end [' Time ' ' Displ. ' ' Vel.'] [t' x' y'] plot(t,x) xlabel('Time') ylabel('Displacement') 'maximum amplitude' A=norm(x,inf)