利用高斯消元法计算各杆长，定义各杆的参数，实时模拟四杆运动轨迹，输出运动的动画matlab.zip

%获取真实动画模拟数据 mm=90; for i=1:mm theta(i)=2*i*pi/mm; bx0=la*cos(theta(i));by0=la*sin(theta(i)); Vb(i)=complex(bx0,by0); %角A theta_a(i)=angle(Vb(i)); %向量BD，lbd为BD杆瞬时长 Vbd(i)=Vb(i)-Vd; theta_bda(i)=pi-angle(Vbd(i)); lbd(i)=abs(Vbd(i)); theta_cdb(i)=acos(((lc*lc+lbd(i)*lbd(i))-lb*lb)/(2*lc*lbd(i))); %角D theta_d(i)=theta_bda(i)-theta_cdb(i); %C的向量 cx0=ld-lc*cos(theta_d(i));cy0=lc*sin(theta_d(i)); Vc(i)=complex(cx0,cy0); %求C杆的角速度 Vbc(i)=Vc(i)-Vb(i); theta_2(i)=angle(Vb(i)); theta_3(i)=angle(Vbc(i)); theta_4(i)=pi-theta_d(i); w2=w; %列矩阵方程求解w3,w4 R1=[-lb*sin(theta_3(i)) lc*sin(theta_4(i));lb*cos(theta_3(i)) -lc*cos(theta_4(i))]; R2=[w2*la*sin(theta_2(i));w2*la*sin(theta_2(i))]; Result_w=R1\R2; w3(i)=Result_w(1,1); w4(i)=Result_w(2,1); %求C杆的角加速度 Ra=[-lb*sin(theta_3(i)) lc*sin(theta_4(i));-lb*cos(theta_3(i)) -lc*cos(theta_4(i))]; Ras=[w2*w2*la*cos(theta_2(i))+w3(i)*w3(i)*lb*cos(theta_3(i))-w4(i)*w4(i)*lc*cos(theta_4(i)) w2*w2*la*sin(theta_2(i))+w3(i)*w3(i)*lb*sin(theta_3(i))-w4(i)*w4(i)*lc*sin(theta_4(i))]; Result_a=Ra\Ras; a3(i)=Result_a(1,1); a4(i)=Result_a(2,1); end