wilson_wilson_Wilson-θ法matlab_

%提取参数 m1=input('请输入质量块1的质量：') m2=input('请输入质量2的质量：') c=input('请输入阻尼系数：') k1=input('请输入弹簧1的刚度系数：') k2=input('请输入弹簧2的刚度系数：') f=input('请输入激励力（简谐激励）：') x1=input('请输入质量块1的初始位移：') x2=input('请输入质量块2的初始位移：') v1=input('请输入质量块1的初始速度：') v2=input('请输入质量块2的初始速度：') theta=input('请输入算法系数θ（参考值1.4）：') timespan=input('请输入模拟的时长：') dt=input('请输入求解步长：') %存储各量 M = [m1 0 ; 0 m2 ]; C = [c -c ; -c c]; K = [k1+k2 -k2;-k2 k2]; X = zeros(2,timespan/dt+1); X(:,2) = [x1;x2]; V = zeros(2,timespan/dt+1); V(:,2) = [v1;v2]; A = zeros(2,timespan/dt+1); F = zeros(2,timespan/dt+1); F(1,:) = 0; F0 = [0;f]; F1 = subs(F0,0); A(:,1) = inv(M)*(F1-C*V(:,1)-K*X(:,1)); %定义系数矩阵和系数 d = 6/(theta*dt)^2; e = 3/(theta*dt); f = theta*dt/2; KS = K+d*M+e*C; for i=1:timespan nt=i*dt; F1=subs(F0,nt); F(:,i+1)=F1; F2=F(:,i)+theta*(F(:,i+1)-F(:,i))+M*(d*X(:,i)+2*e*V(:,i)+2*A(:,i))+C*(e*X(:,i)+2*V (:,i)+f*A(:,i)); h=inv(KS)*F2; A(:,i+1)=d/theta*(h-X(:,i))-d*dt*V(:,i)+(1-3/theta)*A(:,i); V(:,i+1)=V(:,i)+dt/2*(A(:,i+1)+A(:,i)); X(:,i+1)=X(:,i)+V(:,i)*dt+dt^2/6*(A(:,i+1)+2*A(:,i)); end t = (0:1:timespan/dt); plot(t,X(1,:),t,X(2,:)); legend('物块1','物块2'); xlabel('时间/s'); ylabel('位移/m'); title('位移时程曲线');