操作方法
原理如下图公式所示： 其中：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+1 t=(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); end close all figure plot(T,x(1,:)*1000 ,'--b') ylabel('{\itx}_1/mm') figure plot(T,x(2,:)*1000,'-b') ylabel('{\itx}_2/mm') figure plot(T,v(1,:),'--g') ylabel('{\itv}_1/m/s') figure plot(T,v(2,:),'-g') ylabel('{\itv}_2/m/s') figure plot(T,aa(1,:),'--r') ylabel('{\ita}_1/m^2/s') figure plot(T,aa(2,:),'-r') ylabel('{\ita}_2/m^2/s'
结果显示：在matlab命令窗口输入以上命令，运行上述命令，就可以得到系统的位移、速度、加速度纽马克法数值解，我们用MATLAB画图工具，直观的显示出系统两质量块的位移、速度、加速图时间图，如下图所示：