电脑、MATLAB
原理如下图公式所示：其中：Newmark-β正确选择两个参数β和γ很重要；当γ>=1/2, β>=γ/2时，Newmark-β法是无条件稳定的。
代入系统的增量运动方程式，得到Newmark-β法的一系列等效参数，见下图公式所示：
算例演示：问题：下图所示：一无阻尼两自由度受迫振动系统，已知：m1=1kg、m2=2kg,1=k2=k3=1000N/m,系统在m2上作用有一阶跃激振力f(t)=10N，求解该系统的动力响应。
MATLAB程序编写如下：m=[1,0;0,2]; %% 或m=diag([1,2]),输入质量矩阵c=[0,0;0,0]; %% 或c=zeros(2),输入质量矩阵k=[2000,-1000;-1000,2000]; %%输入刚度矩阵f0=[0;0]; %%输入激励x0=[0;0]; %%输入初始位移v0=[0;0]; %%输入初始速度aa0=inv(m)*[f0-c*v0-k*x0]; %%计算初始加速度tm=2;dt=0.01;nt=tm/dt;x=[x0,zeros(2,nt)];v=[v0,zeros(2,nt)];aa=[aa0,zeros(2,nt)];T=[0:dt:tm];alfa=0.25;beta=0.5;a0=1/alfa/dt/dt;a1=beta/alfa/dt;a2=1/alfa/dt;a3=1/2/alfa-1;a4=beta/alfa-1;a5=dt/2*(beta/alfa-2);a6=dt*(1-beta);a7=dt*beta;for i=2:nt+1t=(i-1)*dt;fi=[0;10];f=fi;ke=k+a0*m+a1*c;fe=f+m*(a0*x(:,i-1)+a2*v(:,i-1)+a3*aa(:,i-1))+c*(a1*x(:,i-1)+a4*v(:,i-1)+a5*aa(:,i-1));x(:,i)=ke\fe;aa(:,i)=a0*(x(:,i)-x(:,i-1))-a2*v(:,i-1)-a3*aa(:,i-1);v(:,i)=v(:,i-1)+a6*aa(:,i-1)+a7*aa(:,i);endclose allfigureplot(T,x(1,:)*1000 ,'--b')ylabel('{\itx}_1/mm')figureplot(T,x(2,:)*1000,'-b')ylabel('{\itx}_2/mm')figureplot(T,v(1,:),'--g')ylabel('{\itv}_1/m/s')figureplot(T,v(2,:),'-g')ylabel('{\itv}_2/m/s')figureplot(T,aa(1,:),'--r')ylabel('{\ita}_1/m^2/s')figureplot(T,aa(2,:),'-r')ylabel('{\ita}_2/m^2/s'
结果显示：在matlab命令窗口输入以上命令，运行上述命令，就可以得到系统的位移、速度、加速度纽马克法数值解，我们用MATLAB画图工具，直观的显示出系统两质量块的位移、速度、加速图时间图，如下图所示：
命令语法一定要正确
参数β和γ合理选择很重要