MATLAB的曲柄滑块和四杆机构的综合设计解析
基于matlab的四杆机构运动分析
基于matlab的四杆机构运动分析一、四杆机构基本概念四杆机构是一种通过变换连杆长度，改变机构运动形态的机械系统。
四杆机构通常由固定连杆、推动连杆、连接杆和工作连杆四个连杆组成，其中固定连杆和推动连杆固定不动，连接杆和工作连杆则沿固定轴线的方向做平动或旋转运动。
四杆机构的基本构造如下图所示：四杆机构的四个连杆的长度和构造参数，以及驱动机构的运动决定了机构的运动特性。
在进行四杆机构运动分析时，需要通过求解运动学关系式和动力学方程，得到连杆的运动规律和力学特性。
二、四杆机构运动学分析1.运动学基本方程四杆机构的运动学分析基本方程是连杆长度变化的定理，即:l₁²+l₂²-2l₁l₂cosθ₂=l₃²+l₄²-2l₃l₄cosθ₄其中，l₁，l₂分别为固定连杆和推动连杆长度；l₃，l₄分别为连接杆和工作连杆长度；θ₂，θ₄分别为推动连杆和工作连杆的夹角。
2.运动学求解方法根据四杆机构运动学基本方程，可以求解机构中任意连杆的角度和位置，从而分析机构运动规律。
在matlab程序中，运动分析可以采用分析法或图解法。
分析法通常采用向量法或坐标法，即将四杆机构中各连杆和运动副的运动量表示为向量或坐标，然后根据连杆长度变化的定理，求解四个未知角度θ₁、θ₂、θ₃、θ₄。
图解法则先通过画图确定机构的运动规律，在图上求解连杆的角度。
比如可以采用伯格（Bourgeois）图法或恰普利恩（Chaplygin）图法等。
四杆机构动力学分析基本方程包括平衡方程和力平衡方程。
平衡方程：当四杆机构处于平衡状态时，连杆的受力关系可以表示为：ΣF=0其中ΣF为各连杆受力的合力。
ΣF=m×a其中，m为每个连杆的质量，a为连杆的加速度。
四杆机构动力学求解方法以matlab为工具，可借助matlab的求解器完成求解。
具体可以利用matlab的优化工具箱、控制工具箱和系统动态学工具箱等，来实现机构模型的动态模拟、仿真和优化设计。
基于matlab的曲柄滑块机构设计与运动分析_陈长秀
变，从第 i+1 个功能块开始逐位交换。
（3）变异运算的改进
由于在每个功能块中，“1”的数目即是该题型试题的数目， 因此在变异过程中应保证整个种群所有功能块中“1”的数目不 变。可执行如下过程，首先，由变异概率决定某位取反；然后，检 查、修正字符串中“1”的数目，保证不发生变化。
（4）用全局最优解替换本次迭代的最差解 为保证好的字符串不至于流失，每次遗传操作前记录本次 迭代的最优解，若该解优于全局最优解则替换全局最优解，否 则全局最优解保持不变。此次遗传操作后，用全局最优解换本 代的最差解。
（上接第 29 页）
图 1 所示的偏置曲柄滑块机构。设 l1=50mm，l2=100mm， e=20mm，w1=2rad/s，设 φ1 的初始值为 0 , 则 φ1 变化时，杆 2 的角位移、角速度和角加速度以及滑块 3 的位移、速度和加速
>> plot(t,xc,t,vc,t,ac);
度的变化值可计算求得，曲柄转角 φ1 在 0- 360°之间变化时， 在 matlab 的计算窗口输入算式后，滑块 3 的位移、速度和加速
2012 年 1 月 第 1 期（总第 158 期）
轻工科技
LIGHT INDUSTRY SCIENCE AND TECHNOLOGY
机械与电气
基于 ｍ ａ ｔｌａ ｂ 的曲柄滑块机构设计与运动分析
陈长秀
（陕西国防工业职业技术学院，陕西 西安 ７１ ０３０２）
【摘 要】 建立了曲柄滑块机构的计算模型，并使用 matlab 对曲柄滑块机构进行了运动分析，提高了设计效率和设计精度。
图 1 偏置曲柄滑块机构 建立坐标系如图 1 所示，由曲柄滑快机构的矢量封闭图[1] 可得：
φl1 cosφ1+l2 cosφ2=xc
基于MATLAB的四杆机构运动分析
石河子大学毕业设计（论文)题目:基于MATLAB的四杆机构运动分析与动画模拟系统院（系）：机械电气工程学院专业：机械设计制造及其自动化学号: 2002071189姓名: 娄元建指导教师：葛建兵完成日期：二零零六年五月基于MATLAB的四杆机构运动分析与动画模拟系统[摘要] 本文介绍MATLAB开发机构运动分析和动画模拟系统的方法，并且利用MATLAB软件实现平面四杆机构的运动仿真。
以MATLAB程序设计语言为平台，将参数化设计与交互式相结合,设计出四杆机构仿真系统，能够实现四杆机构的参数化设计，并且能够进行机构的速度和加速度分析。
系统具有方便用户的良好界面，并给出界面设计程序,从而使机构分析更加方便、快捷、直观和形象，设计者只需输几参数就可得到仿真结果，为平面四杆机构的设计与分析提供一条便捷的途径。
[关键词] 机构；运动分析;动画模拟;仿真；参数化;MATLABAbstract：The kinematical analysis and animation method of the mechanism using MATLAB was discussed in the paper ， and the kinematic simulation of planar four-bar mechanism with software MATLAB . And emulational system was developed ， the system adopted Matlab as a design ， It combined parametic design with interactive design and had good interface for user ， that can realize parametic design of four-bar mechanism , also to make real speed and acceleration of mechanism 。
基于MATLAB的平面四连杆机构优化设计
基于 MATLAB 的四连杆机构的优化设计
陈伟斌
（汕头大学，工学院）
[摘要] 对平面四连杆机构进行数学建模，要求实现预期的传递函数运动轨迹。利用 MATLAB 强大的运算功能，快速精确地计 算出优化结果。再利用 MATLAB 编写程序检验得出的运动轨迹是否达到期望目标。 [关键词] 连杆、轨迹、优化设计、MATLAB。
Optimized design for four bar linkage mechanism of crushing machine based on MATLAB
Terry Chen (Shantou University, Engineering College)
[Abstract] Analyze the model of four bar linkage mechanism and try to satisfy the movement locus that we excepted. With the strong functions of MATLAB, we can calculate and get the best result quickly. Then write a program to simulate the movement locus of the output and examine whether it satisfy our requirement. [Key Words] Linkage, Movement locus , Optimized Design, MATLAB
l 2, l 3 两 个 独 立 变 量 。 设
l 2 x1; l 3 x 2; 可以得出本题是二维优化问题。
有志，有恒，有识，有为
matlab曲柄连杆机构分析讲课讲稿
m a t l a b曲柄连杆机构分析clear;clc;n=750;l=0.975;R=0.0381;h=0.2;omiga=n.*pi/30;tmax=2.*pi/omiga;t=0:0.001:tmax; %计算曲柄转一圈的总t值alpha1=atan((h+R.*sin(omiga.*t))./sqrt(l.*l-(h+R.*sin(omiga.*t))))+pi;alpha1p=-(R.*omiga.*cos(omiga.*t))./(l.*cos(alpha1));vb=-R.*omiga.*sin(omiga.*t)+R.*omiga.*cos(omiga.*t).*tan(alpha1);ab=-R.*omiga.^2.*cos(omiga.*t)-(R.*omiga.*cos(omiga.*t)).^2./(l.*(cos(alpha1)).^3)-R.*omiga.^2.*sin(omiga.*t).*tan(alpha1);subplot(1,2,1);plot(t,vb);title('曲柄滑块机构的滑块v-t图');xlabel('时间t(曲柄旋转一周)');ylabel('滑块速度v');grid on;subplot(1,2,2);plot(t,ab);title('曲柄滑块机构的滑块a-t图');xlabel('时间t(曲柄旋转一周)');ylabel('滑块加速度a');grid on;%下面黄金分割法求滑块的速度与加速度最大值epsilon=input('根据曲线初始区间已确定，请输入计算精度epsilon(如输入0.001):');a=0;b=0.04; %初始区间n1=0; %n1用于计算次数a1=b-0.618*(b-a);y1=-R.*omiga.*sin(omiga.*a1)+R.*omiga.*cos(omiga.*a1).*tan(alpha1);a2=a+0.618*(b-a);y2=-R.*omiga.*sin(omiga.*a2)+R.*omiga.*cos(omiga.*a2).*tan(alpha1);while abs(a-b)>=epsilonif y1<=y2b=a2;a2=a1;y2=y1;a1=b-0.618*(b-a);y1=-R.*omiga.*sin(omiga.*a1)+R.*omiga.*cos(omiga.*a1).*tan(alpha1);elsea=a1;a1=a2;y1=y2;a2=a+0.618*(b-a);y2=-R.*omiga.*sin(omiga.*a2)+R.*omiga.*cos(omiga.*a2).*tan(alpha1);endn1=n1+1;endvbm1=omiga*(a+b)/2;disp(['经过',num2str(n1),'次计算,用黄金分割法找到速度最大值对应的wt是:', num2str(vbm1),'弧度。
基于matlab的平面四连杆机构设计以及该机构的运动分析参考模板
基于matlab的平面四连杆机构设计以及该机构的运动仿真分析摘要四连杆机构因其结构方便灵活,能够传递动力并实现多种运动形式而被广泛应用于各个领域,因此对其进行运动分析具有重要的意义。
传统的分析方法主要应用几何综合法和解析综合法,几何综合法简单直观,但是精确度较低;解析法精确度较高,但是计算工作量大。
随着计算机辅助数值解法的发展,特别是MATLAB软件的引入,解析法已经得到了广泛的应用。
对于四连杆的运动分析，若应用MATLAB 则需要大量的编程，因此我们引入proe软件，我们不仅可以在此软件中建立实物图，而且还可以对其进行运动仿真并对其运动分析。
在设计四连杆时，我们利用解析综合法建立数学模型，再根据数学模型在MATLAB中编程可以求得其他杆件的长度。
针对范例中所求得的各连杆的长度，我们在proe软件中画出其三维图（如图4）并在proe软件中进行仿真分析得出CB,的角加速度的变化，从而得到CB,两接触处所受到的力是成周期性变化的，可以看出CB,两点处的疲劳断裂，我们提B,两点处极易疲劳断裂，针对C出了在设计四连杆中的一些建议。
关键字：解析法 MATLAB 软件 proe 软件 运动仿真建立用解析法设计平面四杆机构模型对于问题中所给出的连架杆AB 的三个位置与连架杆CD 的三个位置相对应，即三组对应位置为：332211,,,,,ψϕψϕψϕ，其中他们对应的值分别为： 52,45,82,90,112,135，为了便于写代数式，可作出AB 与CD 对应的关系，其图如下：图—2 AB 与CD 三个位置对应的关系通过上图我们可以通过建立平面直角坐标系并利用解析法来求解，其直角坐标系图如下：φααi θi φi图—3 平面机构直角坐标系通过建立直角坐标系OXY ，如上图所示，其中0α与0φ为AB 杆与CD 杆的初始角，各杆件的长度分别用矢量d c b a ,,,,表示，将各矢量分别在X 轴与Y 轴上投影的方程为⎩⎨⎧=++=+)sin(*)sin(*)sin(*)cos(*)cos(*)cos(*φθαφθαc b a c d b a在上述的方程中我们可以消除θ，从而可以得到α与φ之间的关系如下：)cos(2)cos(2)cos(2)(2222αφαφab ac cd b d c a +-=+-++ (1) 为便于化简以及matlab 编程我们可以令：⎪⎪⎪⎩⎪⎪⎪⎨⎧==-++=c d H a d H ac b d c a H 32222212 (2) 通过将（2）式代入（1）式中则可以化简得到如下等式： )cos()cos()cos(321αφαφH H H +-=+ （3）我们可以通过（3）式将两连架杆对应的位置带入（3）式中，我们可以得到如下方程：⎪⎩⎪⎨⎧+-=++-=++-=+)cos()cos()cos()cos()cos()cos()cos()cos()cos(333332123222211311121ϕψϕψϕψϕψϕψϕψH H H H H H H H H （4） 联立（4）方程组我们可以求得321,,H H H ，再根据（2）中的条件以及所给定的机架d 的长度，我们可以求出其它杆件的长度为：⎪⎪⎪⎩⎪⎪⎪⎨⎧-++===1222322acH d c a b H d c H d a （5）四连杆设计范例：在日常生活中，我们经常看到消防门总能自动关上，其实它是利用四连杆机构与弹簧组成的。
运用MATLAB解决四杆机构问题
MATLAB 解题1.设有如图所示四杆机构，其中→R 4为机架（常矢），→R1为主动杆，→R3为从动杆，→R 2为连杆。
设在某一工作位置时各杆的角速度和角加速度分别取如下值：ω1=20 rad/s, ε1= 0；ω2=8.5 rad/s, ε2=－10 rad /s 2；ω3=13 rad/s, ε3=－160rad /s 2.试根据上述要求确定该机构尺寸比。
根据图（2），回路闭合方程可写为：→R 1 +→R 2 +→R 3＝－→R 4 回路闭合方程对时间求导一次，利用（6）式，可得： 图2 ω1→R 1 +ω2→R 2 +ω3→R 3 ＝ 0回路闭合方程对时间求导两次，利用（7）式，可得c 1→R 1 + c 2 →R 2 + c 3→R 3 ＝ 0其中 c 1=ε1+j ω12 ， c 2=ε2+j ω22， c 3=ε3+j ω32解关于→R 1 ，→R 2 和→R 3的线性方程组：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡→→→001111321321321R R R c c c ωωω→R 4 （13） 可得 →R 1=DDx →R 4， →R 2=DDy →R 4 ， →R 3=DDz →R 4注意到上述解中含有相同的分母D,它是一个复数，不妨记为D ＝k<j α|，被它除的效果是把各杆的长度都缩小k 倍，同时方向都顺时针旋转α角，相当于机构不动，坐标轴逆时针旋转α角。
设计机构时，重要的是机构的形状与尺寸比例。
基于这种考虑，可设→R 4 / D ＝1，则有→R 1=D x =32320111c c ωω-=1230－j497.3 ； →R 2= D y =311030111c c ωω-＝－3200－j1820 ； →R 3= D z =001112121c c ωω-＝200+j1955 . 于是：→R 4 = －(→R 1 +→R 2+→R 3) = 1770+j362.3在坐标系上作出上述各杆矢量图，根据各杆矢量图作出机构的闭合矢量图，再根据实际需要选定某一杆长度，其它各杆长度按图比例相似放大。
物理-曲柄滑块机构的运动分析-matlab
子函数%子函数slider_crank文件function[theta2,s3,omega2,v3,alpha2,a3]=slider_crank(theta1,omega 1,alpha1,l1,l2,e)%计算连杆2的角位移和滑块3的线位移theta2=asin((e-l1*sin(theta1))/l2);s3=l1*cos(theta1)+l2*cos(theta2);%计算连杆2的角为速度和滑块的线速度A=[-l1*sin(theta1),1;-2*cos(theta2),0];B=[-l1*sin(theta1);l1*cos(theta1)];omega=A\(omega1*B);omega2=omega(1);v3=omega(2);%计算连杆2的角加速度和滑块3的线加速度At=[omega2*l2*cos(theta2),0;omega2*l2*sin(theta2),0];Bt=[-omega1*l1*cos(theta1);-omega1*l1*sin(theta1)];alpha=A\(-At*omega+alpha1*B+omega1*Bt);alpha2=alpha(1);a3=alpha(2);主函数%住程序slider_crank_main文件%输入已经知道的数据clear;l1=100;l2=300;e=0;hd=pi/180;du=180/pi;omega1=10;alpha1=0;%调用子函数slider_ank计算曲柄滑块机构位移，速度，加速度for n1=1:720theta1(n1)=(n1-1)*hd;[theta2(n1),s3(n1),omega2(n1),v3(n1),alpha2(n1),a3(n1)]=slider_cr ank...(theta1(n1),omega1,alpha1,l1,l2,e);end%位移，速度，加速度和曲柄滑块机构图形输出figure(l1);n1=1:720;subplot(2,2,1); %绘制位移图[AX,H1,H2]=plotyy(theta1*du,theta2*du,theta1*du,s3);set(get(AX(1),'ylabel'),'String','连杆角位移/\circ')set(get(AX(2),'ylabel'),'String','滑块位移/mm')title('位移线图');xlabel('曲柄转角\theta_1/\circ')grid on;subplot(2,2,2); %绘制速度图[AX,H1,H2]=plotyy(theta1*du,omega2,theta1*du,v3);set(get(AX(2),'ylabel'),'String','滑块速度/mm\cdots^{-1}') title('速度线图');xlabel('曲柄转角\theta_1/\circ')ylabel('连杆角速度/rad\cdots^{-1}')grid on;subplot(2,2,3); %绘制加速度图[AX,H1,H2]=plotyy(theta1*du,alpha2,theta1*du,a3);set(get(AX(2),'ylabel'),'String','滑块加速度/mm\cdots^{-2}') title('加速度线图');xlabel('曲柄转角\theta_1/\circ')ylabel('连杆加速度/rad\cdots^{-2}')grid on;subplot(2,2,4);%绘曲柄滑块机构图x(1)=0;y(1)=0;x(2)=l1*cos(70*hd);y(2)=l1*sin(70*hd);x(3)=s3(70);y(3)=e;x(4)=s3(70);y(4)=0;x(5)=0;y(5)=0;x(6)=x(3)-40;y(6)=y(3)+10;x(7)=x(3)+40;y(7)=y(3)+10;x(8)=x(3)+40;y(8)=y(3)-10;x(9)=x(3)-40;y(9)=y(3)-10;x(10)=x(3)-40;y(10)=y(3)+10;i=1:5;plot(x(i),y(i));grid on;hold on;i=6:10;plot(x(i),y(i));title('曲柄滑块机构');grid on;hold on;xlabel('mm');ylabel('mm')axis([-50 400 -20 130]); plot(x(1),y(1),'o');plot(x(2),y(2),'o');plot(x(3),y(3),'o');%曲柄滑块的仿真运动figure(2)m=moviein(20);j=0;for n1=1:5:360j=j+1;clf;%x(1)=0;y(1)=0;x(2)=l1*cos(n1*hd); y(2)=l1*sin(n1*hd); x(3)=s3(n1);y(3)=e;x(4)=(l1+l2+50);y(4)=0;x(5)=0;y(5)=0;x(6)=x(3)-40;y(6)=y(3)+10;x(7)=x(3)+40;y(7)=y(3)+10;x(8)=x(3)+40;y(8)=y(3)-10;x(9)=x(3)-40;y(9)=y(3)-10;x(10)=x(3)-40;y(10)=y(3)+10;%i=1:3;plot(x(i),y(i));grid on;hold on;i=4:5;plot(x(i),y(i));i=6:10;plot(x(i),y(i));plot(x(1),y(1),'o');plot(x(2),y(2),'o');plot(x(3),y(3),'o');xlabel('mm');ylabel('mm')axis([-50 450 -150 150]); m(j)=getframe;endmovie(m)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《计算机仿真技术》课程设计报告姓名：冯叶 / 浦合昀学号： 201410302544/ 201410302547专业班级：机械卓越141 指导教师：刘孝保2015年 6月目录目录1.仿真问题描述 ...................................................................................................................................................2.仿真问题数学模型 ...........................................................................................................................................3.Matlab实现方法 ..............................................................................................................................................4.Matlab代码 ......................................................................................................................................................5.仿真结论 ...........................................................................................................................................................6.遇到的问题和解决的方式 ...............................................................................................................................7.课程学习意见与建议 .......................................................................................................................................1.仿真问题描述已知机架AD 长为L1，曲柄AB 长为L2，连杆BC 长L3，另一机架长CD 长为L4，与AB 杆相连的是一滑块E 。
BE 杆长为L5，设计一个四杆加滑块的机构，其中L1-L5杆长可变。
并且可以通过输入的杆长，来判别，该机构到底可不可行。
L3L4L2 L5L12.仿真问题数学模型（1）四杆机构的设计：在用矢量法建立机构的位置方程时，需将构件用矢量来表示，并作出机构的封闭矢量多边形。
如图1所示，先建立一直角坐标系。
设各构件的长度分别为1L 、2L 、3L 、4L ，其方位角为1θ、2θ、3θ、 4θ。
以各杆矢量组成一个封闭矢量多边形，即ABCDA 。
其个矢量之和必等于零。
易知：2314L L L L +=+角位移方程的分量形式为：2233114422331144cos cos cos cos sin sin sin sin L L L L L L L L θθθθθθθθ+=+⎧⎫⎨⎬+=+⎩⎭ 要求th3,那么 ()'''sin cos cos cos sin sin d d dt dt d d dt dt uv vu uv θθθωθθθθωθ⎧⎫==⎪⎪⎪⎪⎪⎪=-=-⎨⎬⎪⎪⎪⎪=+⎪⎪⎩⎭在角位移方程分量形式中，由于假定机架为参考系，矢量1与x 轴重合，1θ=0，则有非线性超越方程组：1342233144234223344(,)cos cos cos 0(,)sin sin sin 0f L L L L f L L L θθθθθθθθθθ=+--=⎧⎫⎨⎬=+-=⎩⎭可以借助牛顿-辛普森数值解法或Matlab 自带的fsolve 函数求出连杆3的角位移3θ和摇杆4的角位移4θ。
A B C E D求解具有n 个未知量i x （i=1,2,…,n ）的线性方程组：111122112111221211221n n n n n n nn n n a x a x a x b a x a x a x b a x a x a x b +++=⎧⎫⎪⎪+++=⎪⎪⎨⎬⎪⎪⎪⎪+++=⎩⎭式中，系列矩阵A 是一个*n n 阶方阵：1111n m mn a a A a a ⎛⎫ ⎪= ⎪ ⎪⎝⎭A 的逆矩阵为1A -;常数项b 是一个n 维矢量：12(,,,)T n b b b b =因此，线性方程组解的矢量为：12(,,,)T T n x x x x A b == 非线性超越方程组是求解连杆3和摇杆4角速度和角加速度的依据。
（2）曲柄滑块的设计：由图可知，C 滑块的位移总是与AB,BC 和他们之间的角度存在着一定的关系，关系如下：LAC=AB ×cos (∅)+√BC 2—AB 2×cos∅2通过以上这个式子，我们就可以来求C 点的位移，速度，加速度。
3.Matlab 实现方法（1）怎么设计四杆机构：创建函数FoutBarPosition ，函数fsolve 通过他确定34,θθ，然后知道34,θθ后，来求取各个点的坐标，通过plot 命令在指定的区域内连线，取点，画图。
（2）怎么设计曲柄滑块机构：通过解方程的方法，用solve 来求取C 滑块的坐标，用diff 函数求取C 滑块的速度，加速度曲线，通过plot 命令在指定的区域内连线，取点，画图。
4.Matlab 代码（1）建新的函数在点m文件中：function t=fourbarposition(th,th2,L2,L3,L4,L1)t=[L2*cos(th2)+L3*cos(th(1))-L4*cos(th(2))-L1;…L2*sin(th2)+L3*sin(th(1))-L4*sin(th(2))];（2）主程序如下：%获取杆长l1=str2double(get(handles.edit1,'string'));l2=str2double(get(handles.edit2,'string'));l3=str2double(get(handles.edit3,'string'));l4=str2double(get(handles.edit4,'string'));l5=str2double(get(handles.edit8,'string'));%滑块和四杆机构的设计syms t s; %定义变量f=l5^2-l2^2-s^2+2*l2*s*cos(t);ff=solve(f,s);vv=diff(ff,1);aa=diff(ff,2);th2=0:pi/15:6*pi;times=length(th2);for i=1:91wyy(1,i)=eval(subs(ff(1),t,th2(i)));wyy(2,i)=eval(subs(ff(2),t,th2(i)));vyy(1,i)=eval(subs(vv(1),t,th2(i)));vyy(2,i)=eval(subs(vv(2),t,th2(i)));ayy(1,i)=eval(subs(aa(1),t,th2(i)));ayy(2,i)=eval(subs(aa(2),t,th2(i)));endfor i=1:timesif wyy(1,i)>0wy(i)=wyy(1,i);elsewy(i)=wyy(2,i);endif vyy(1,i)>0vy(i)=vyy(1,i);elsevy(i)=vyy(2,i);endif ayy(1,i)>0ay(i)=ayy(1,i);elseay(i)=ayy(2,i);endendth34=zeros(length(th2),2); %%建立一个N行2列的零矩阵options=optimset('display','off');for m=1:length(th2) %用fsove函数求解关于th3,th4的非线性超越方程，结果保存在th34中th34(m,:)=fsolve('fourbarposition',[1 1],...options,th2(m),l2,l3,l4,l1);end%求各个的坐标Ex=wy;Ey=zeros(size(th2));Cy=l2*sin(th2)+l3*sin(th34(:,1)');Cx=l2*cos(th2)+l3*cos(th34(:,1)');Bx=[l2*cos(th2)];By=[l2*sin(th2)];Ax=zeros(size(th2));Ay=zeros(size(th2));Dx=l1+zeros(size(th2));Dy=zeros(size(th2));Ev=vy;Ew=zeros(size(th2));Ea=ay;En=zeros(size(th2));%求位移,速度,加速度的范围:g=[Ax Bx Cx Dx Ex];m=[Ay By Cy Dy Ey];maxX=max(g);minX=min(g);maxY=max(m);minY=min(m);maxwy=max(Ex);minwy=min(Ex);maxvy=max(Ev);minvy=min(Ev);maxay=max(Ea+50);minay=min(Ea-50);%画动画图for i=1:timesaxes(handles.axes1);plot([Ax(i),Bx(i)],[Ay(i),By(i)],'lineWidth',3); hold onplot([Bx(i),Cx(i)],[By(i),Cy(i)],'lineWidth',3); plot([Ax(i),Dx(i)],[Ay(i),Dy(i)],'lineWidth',3); plot([Cx(i),Dx(i)],[Cy(i),Dy(i)],'lineWidth',3); plot([Bx(i),Ex(i)],[By(i),Ey(i)],'lineWidth',3); plot([-10000,Ax(i)],[0,Ay(i)],'lineWidth',3);plot([10000,Ax(i)],[0,Ay(i)],'lineWidth',3);plot([Ex(i)+10,Ex(i)-10],[Ey(i)+10,Ey(i)+10]);plot([Ex(i)+10,Ex(i)-10],[Ey(i)-10,Ey(i)-10]);plot([Ex(i)-10,Ex(i)-10],[Ey(i)+10,Ey(i)-10]);plot([Ex(i)+10,Ex(i)+10],[Ey(i)+10,Ey(i)-10]);plot(0,0,'or','lineWidth',3)plot(Bx(i),By(i),'or','lineWidth',3)plot(Cx(i),Cy(i),'or','lineWidth',3)plot(Dx(i),Dy(i),'or','lineWidth',3)plot(Ex(i),Ey(i),'or','lineWidth',3)axis equal;axis([minX,maxX,minY,maxY]);axis off;hold off;pause(0.1)%画滑块位移图axes(handles.axes3);plot(th2(1:i),wy(1:i),'lineWidth',3);axis([0,2*pi,minwy,maxwy])pause(0.1)%画滑块速度图axes(handles.axes4);plot(th2(1:i),vy(1:i),'lineWidth',3);axis([0,2*pi,minvy,maxvy])pause(0.1)%画滑块加速度图axes(handles.axes5);plot(th2(1:i),ay(1:i),'lineWidth',3);axis([0,2*pi,minay,maxay])pause(0.1)end%判断杆长的代码：l1=str2double(get(handles.edit1,'string'));l2=str2double(get(handles.edit2,'string'));l3=str2double(get(handles.edit3,'string'));l4=str2double(get(handles.edit4,'string'));lall=l1+l2+l3+l4;lmax=max([l1 l2 l3 l4]);lmin=min([l1 l2 l3 l4]);if (lmax+lmin)>(lall-lmax-lmin)set(handles.edit5,'string','不符合杆长条件，请重新输入');set(handles.edit6,'string','');set(handles.edit7,'string','');elseif l2>l4set(handles.edit6,'string','不符合A点的周转条件，请重新输入'); set(handles.edit5,'string','');set(handles.edit7,'string','');elseset(handles.edit7,'string','·符合条件，可以运动'); set(handles.edit5,'string','');set(handles.edit6,'string','');endend%关闭代码：close5.仿真结论（1）如图所示为界面。