Newmark法求解单自由度

%单位：N/mm/s/tonfunctionres=Newmark(alpha,C)%系统设置；T=0.1/alpha;K=(2*3.1415926/T)^2;M=1;%C=0;%定义参数h=0.0002;beta=0.25;gamma=0.5;con=zeros(1,7);con(1)=1/(beta*h^2);con(2)=gamma/(beta*h);con(3)=1/(beta*h);con(4)=1/(2*beta)-1;con(5)=gamma/beta-1;con(6)=0.5*h*(gamma/beta-2);con(7)=h*(1-gamma/(2*beta));%有效刚度Ke=K+con(1)*M+con(2)*C;%定义矩形荷载t=0:h:1;f=zeros(1,size(t,2));fori=1:size(t,2)ift(i)==0f(i)=0;elseift(i)0&&t(i)=0.1f(i)=1000*(3.1415926)^2;elsef(i)=0;endend%plot(t,f);%系统初始条件u0=0;du0=0;ddu0=0;U=zeros(3,size(t,2));%求解fori=1:(size(t,2)-1)fe=f(i+1)+M*(con(1)*u0+con(3)*du0+con(4)*ddu0)+C*(con(2)*u0+con(5)*du0+con(6)*ddu0);u1=fe/Ke;du1=con(2)*(u1-u0)-con(5)*du0+con(7)*ddu0;%计算速度和加速度；ddu1=(f(i+1)-C*du1-K*u1)/M;U(:,i+1)=[u1;du1;ddu1];u0=u1;du0=du1;ddu0=ddu1;endres=[U;t];end