下载文档原格式
基于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的四杆机构运动分析与动画模拟系统院（系）：机械电气工程学院专业：机械设计制造及其自动化学号: 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 。
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 则需要大量的编程，因此我们引入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 解题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在坐标系上作出上述各杆矢量图，根据各杆矢量图作出机构的闭合矢量图，再根据实际需要选定某一杆长度，其它各杆长度按图比例相似放大。
子函数%子函数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)。
matlab曲柄滑块机构课程设计一、课程目标知识目标：1. 理解并掌握曲柄滑块机构的基本原理与运动特性；2. 学会使用MATLAB软件进行曲柄滑块机构的运动仿真；3. 掌握运用MATLAB分析曲柄滑块机构的运动数据及性能参数。
技能目标：1. 能够运用所学知识，设计简单的曲柄滑块机构；2. 熟练操作MATLAB软件，进行曲柄滑块机构的运动分析与仿真；3. 能够通过MATLAB处理数据，优化曲柄滑块机构的设计。
情感态度价值观目标：1. 培养学生的团队协作精神，提高沟通与表达能力；2. 激发学生探索科学、技术问题的兴趣，培养创新意识和实践能力；3. 增强学生对机械工程领域的认识和热爱，提高专业认同感。
课程性质：本课程为机械工程专业课程设计，旨在通过实践操作，使学生掌握曲柄滑块机构的设计与分析方法。
学生特点：学生已具备一定的机械原理、力学和MATLAB基础，具有较强的动手能力和学习兴趣。
教学要求：结合实际工程案例，以实践为主，注重培养学生的实际操作能力、分析问题和解决问题的能力。
通过课程学习，使学生能够独立完成曲柄滑块机构的设计与分析任务。
二、教学内容1. 曲柄滑块机构原理及运动特性分析：- 曲柄滑块机构基本组成与工作原理- 曲柄滑块机构的运动学分析- 运动方程的建立及求解2. MATLAB软件在曲柄滑块机构中的应用：- MATLAB软件的基本操作与常用函数- MATLAB曲线拟合、数值计算等功能在曲柄滑块机构分析中的应用- MATLAB/Simulink环境下曲柄滑块机构的运动仿真3. 曲柄滑块机构设计及优化：- 设计原则与步骤- 参数化设计方法- 基于MATLAB的曲柄滑块机构设计优化4. 实践操作与案例分析：- 实际工程案例介绍与分析- 曲柄滑块机构设计及运动分析的实践操作- 数据处理与结果分析教学内容安排与进度：1. 第一周：曲柄滑块机构原理及运动特性分析2. 第二周：MATLAB软件在曲柄滑块机构中的应用3. 第三周：曲柄滑块机构设计及优化4. 第四周：实践操作与案例分析教材章节：1. 《机械原理》中曲柄滑块机构相关章节2. 《MATLAB基础与应用》中相关章节3. 《机械设计》中机构设计及优化相关章节教学内容注重理论与实践相结合，通过系统性的教学，使学生掌握曲柄滑块机构的设计与分析方法，并能够运用MATLAB软件进行实际操作。
KUNMING UNIVERSITY OF SCIENCE AND TECHNOLOGY《计算机仿真技术》课程设计报告冯叶/ 浦合旳201410302544/ 201410302547刘孝保2015年6月姓名: 学号: 专业班级: 指导教师:机械卓目录©区肌理乂殳申KUMWBG sngn OF SCIENCE MO TCCWlOGr目录1 •仿真问题描述.........................................................................2•仿真问题数学模型......................................................................3. Mat lab实现方法 .....................................................................4・Mat lab代码..........................................................................5•仿真结论..............................................................................6.遇到的问题和解决的方式.................................................................7 •课程学习意见与建议...................................................................《计算机仿真技术》课程设计报告 艮咽疗N 乂孝 ItnVH ； WmJSTY ：f SCOtCE MP TOCtXCCf 1 •仿真问题描述已知机架AD 长为L1,曲柄AB 长为L2,连杆BC 长L3,另一机架长CD 长为L4,与AB 杆相 连的是一滑块E 。
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),'弧度。
一、课程设计容及要求：1．对连杆机构运动工作原理及运动参数有一定理解2．掌握MATLAB基本命令3．了解MATLAB编程的基本知识，并能编写简单M文件4．了解MATLAB图形界面设计的基本知识5．课程设计说明书：应阐述整个课程设计容，要突出重点和特色，图文并茂，文字通畅。
应有目录、摘要及关键词、正文、参考文献等容，字数一般不少于6000字。
二、主要参考资料有关复杂刀具参数计算及结构设计、机械制造工艺与设备的手册与图册。
三、课程设计进度安排指导教师（签名）：时间：教研室主任(签名)：时间：院长（签名）：时间：目录1平面连杆机构的运动分析 (1)1.1 机构运动分析的任务、目的和方法 (1)1.2 机构的工作原理 (1)1.3 机构的数学模型的建立 (1)1.3.1建立机构的闭环矢量位置方程 (1)1.3.2求解方法 (2)2 基于MATLAB程序设计 (4)2.1 程序流程图 (4)2.2 M文件编写 (6)2.3 程序运行结果输出 (7)3 基于MATLAB图形界面设计 (11)3.1界面设计 (11)3.2代码设计 (12)4 小结 (17)参考文献 (18)1平面连杆机构的运动分析1.1 机构运动分析的任务、目的和方法曲柄摇杆机构是平面连杆机构中最基本的由转动副组成的四杆机构，它可以用来实现转动和摆动之间运动形式的转换或传递动力。
对四杆机构进行运动分析的意义是：在机构尺寸参数已知的情况下，假定主动件（曲柄）做匀速转动，撇开力的作用，仅从运动几何关系上分析从动件（连杆、摇杆）的角位移、角速度、角加速度等运动参数的变化情况。
还可以根据机构闭环矢量方程计算从动件的位移偏差。
上述这些容，无论是设计新的机械，还是为了了解现有机械的运动性能，都是十分必要的，而且它还是研究机械运动性能和动力性能提供必要的依据。
机构运动分析的方法很多，主要有图解法和解析法。
当需要简捷直观地了解机构的某个或某几个位置的运动特性时，采用图解法比较方便，而且精度也能满足实际问题的要求。
基于matlab GUI的平面四杆机构的运动分析一、目的通过matlab对平面四杆机构进行运动仿真，并以GUI界面方式实现输入输出的参数化，对平面四杆机构进行位置分析、速度分析、加速度分析和静力学分析。
此外，通过动画演示，更加形象直观地观察机构的运动过程。
最后，将程序编译成.exe独立可执行文件，可以在其它没有安装matlab的机器上运行。
二、设计思路通过matlab的GUI功能模块，创建一个图形用户界面，在自动生成的代码框架中对初始化函数和回调函数等进行编辑，建立与控件相关联的程序：控件属性、位置分析、速度分析、加速度分析、静力学分析、动画演示等。
图1是平面四杆机构的示意图，输入角q的运动规律为q=pi/50*t^2+q0，r1、r2是从动角。
对t时刻沿着杆长距离原点A的任意一点进行分析。
注意：输入输出角的单位为度，时间t的取值范围为0：0.05：10，任意点lx的取值范围为0～a1+a2+a3，估算的从动角r1、r2的迭代初始值不能偏离平衡位置太大。
图1、平面四杆机构示意图三、设计流程1、通过GUI模块创建图形用户界面命令方式：在Matlab命令窗口键入>>guide；菜单方式：在Matlab的主窗口中，选择File>New>GUI命令，就会显示GUI的设计模板。
如图1所示。
图2、创建图形界面2、设计图形界面在创建之后的图形界面中插入坐标轴axes，静态文本框static text，编辑文本框edit text，按钮push button等等。
如图所示。
图3、图形界面设计3、编辑回调函数1）位置分析：输入角的函数为：q=pi/50*t^2+q0。
在时间t=0~10s内，每一个时间点估算两个初始从动角，根据牛顿-拉普森迭代得到准确的机构位置。
10s刚好主动角经历了360度，记录每一时刻的位置，便可以动画演示。
2）速度分析：输入角速度为：dq=pi/25*t。
选择杆件上的任意一点（坐标表示为质点沿着杆件到原点A的距离）做分析，正确表达出角速度系数和速度系数，便可以求出质点的速度。
⽤matlab分析四杆机构⾸先创建函数FoutBarPosition，函数fsolve通过他确定。
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))];主程序如下：disp ' * * * * * * 平⾯四杆机构的运动分析* * * * * *'L1=304.8;L2=101.6;L3=254.0;L4=177.8; %给定已知量，各杆长L1,L2,L3,L4th2=[0:1/6:2]*pi; %曲柄输⼊⾓度从0⾄360度，步长为pi/6th34=zeros(length(th2),2); %建⽴⼀个N⾏2列的零矩阵，第⼀列存放options=optimset('display','off'); %θ_3,第⼆列存放θ_3 for m=1:length(th2) %建⽴for循环，求解θ_3，θ_4th34(m,:)=fsolve('fourbarposition',[1 1],…%调⽤fsove函数求解关于θ_3，θ_4options,th2(m),L2,L3,L4,L1); %的⾮线性超越⽅程，结果保存在th34中endy=L2*sin(th2)+L3*sin(th34(:,1)'); %连杆3的D端点Y坐标值x=L2*cos(th2)+L3*cos(th34(:,1)'); %连杆3的D端点X坐标值xx=[L2*cos(th2)]; %连杆3的C端点X坐标值yy=[L2*sin(th2)]; %连杆3的C端点Y坐标值figure(1)plot([x;xx],[y;yy],'k',[0 L1],[0 0],…%绘制连杆3的⼏个位置点'k--^',x,y,'ko',xx,yy,'ks')title('连杆3的⼏个位置点')xlabel('⽔平⽅向')ylabel('垂直⽅向')axis equal %XY坐标均衡th2=[0:2/72:2]*pi; %重新细分曲柄输⼊⾓度θ_2，步长为5度th34=zeros(length(th2),2);options=optimset('display','off');for m=1:length(th2)th34(m,:)=fsolve('fourbarposition',[1 1],…options,th2(m),L2,L3,L4,L1);endfigure(2)plot(th2*180/pi,th34(:,1),th2*180/pi,th34(:,2)) %绘制连杆3的⾓位移关于曲柄2的⾓位移图plot(th2*180/pi,th34(:,1)*180/pi,…th2*180/pi,th34(:,2)*180/pi) %绘制摇杆4的⾓位移关于曲柄2的⾓位移图axis([0 360 0 170]) %确定XY边界值grid %图形加⽹格xlabel('主动件转⾓\theta_2(度)')ylabel('从动件⾓位移(度)')title('⾓位移线图')text(120,120,'摇杆4⾓位移')text(150,40,'连杆3⾓位移')w2=250; %设定曲柄⾓速度for i=1:length(th2)A=[-L3*sin(th34(i,1)) L4*sin(th34(i,2));…L3*cos(th34(i,1)) -L4*cos(th34(i,2))];B=[w2*L2*sin(th2(i)); -w2*L2*cos(th2(i))];w=inv(A)*B;w3(i)=w(1);w4(i)=w(2);endfigure(3)plot(th2*180/pi,w3,th2*180/pi,w4); %绘制⾓速度线图axis([0 360 -175 200])text(50,160,'摇杆4⾓速度(\omega_4)')text(220,130,'连杆3⾓速度(\omega_3)')gridxlabel('主动件转⾓\theta_2(度)')ylabel('从动件⾓速度(rad\cdot s^{-1})')title('⾓速度线图')for i=1:length(th2)C=[-L3*sin(th34(i,1)) L4*sin(th34(i,2));…L3*cos(th34(i,1)) -L4*cos(th34(i,2))];D=[w2^2*L2*cos(th2(i))+w3(i)^2*L3*cos(th34(i,1))-w4(i)^2*L4*cos(th34(i,2));... w2^2*L2*sin(th2(i))+w3(i)^2*L3*sin(th34(i,1))-w4(i)^2*L4*sin(th34(i,2))];a=inv(C)*D;a3(i)=a(1);a4(i)=a(2);endfigure(4)plot(th2*180/pi,a3,th2*180/pi,a4); %绘制⾓加速度线图axis([0 360 -70000 65000])text(50,50000,'摇杆4⾓加速度(\alpha_4)')text(220,12000,'连杆3⾓加速度(\alpha_3)')gridxlabel('从动件⾓加速度')ylabel('从动件⾓加速度(rad\cdot s^{-2})')title('⾓加速度线图')disp '曲柄转⾓连杆转⾓-摇杆转⾓-连杆⾓速度-摇杆⾓速度-连杆加速度-摇杆加速度' ydcs=[th2'*180/pi,th34(:,1)*180/pi,th34(:,2)*180/pi,w3',w4',a3',a4'];disp(ydcs)>> * * * * * * 平⾯四杆机构的运动分析* * * * * *曲柄转⾓连杆转⾓-摇杆转⾓-连杆⾓速度-摇杆⾓速度-连杆加速度-摇杆加速度1.0e+004 *0 0.0044 0.0097 -0.0125 -0.0125 -0.5478 4.84580.0005 0.0042 0.0094 -0.0126 -0.0107 0.2300 5.5630 0.0010 0.0039 0.0092 -0.0124 -0.0086 0.8946 6.0520 0.0015 0.0037 0.0091 -0.0119 -0.0065 1.4143 6.2982 0.0020 0.0034 0.0090 -0.0114 -0.0043 1.7801 6.3174 0.0025 0.0032 0.0089 -0.0107 -0.0021 2.0027 6.1467 0.0030 0.0030 0.0089 -0.0100 0.0000 2.1046 5.8339 0.0035 0.0028 0.0089 -0.0093 0.0020 2.1134 5.4272 0.0040 0.0026 0.0090 -0.0085 0.0038 2.0566 4.9687 0.0045 0.0025 0.0091 -0.0078 0.0054 1.9578 4.4918 0.0050 0.0023 0.0092 -0.0072 0.0069 1.8356 4.0198 0.0055 0.0022 0.0093 -0.0065 0.0082 1.7040 3.5680 0.0060 0.0021 0.0095 -0.0060 0.0094 1.5725 3.1450 0.0065 0.0019 0.0097 -0.0055 0.0104 1.4474 2.7545 0.0070 0.0018 0.0099 -0.0050 0.01131.33282.3968 0.0075 0.0017 0.0102 -0.0045 0.0121 1.2307 2.0702 0.0080 0.0017 0.0104 -0.0041 0.0128 1.1425 1.7716 0.0085 0.0016 0.0107 -0.0037 0.0134 1.0687 1.4971 0.0090 0.0015 0.0110 -0.0034 0.0138 1.0095 1.2426 0.0095 0.0014 0.0112 -0.0030 0.0142 0.9653 1.0035 0.0100 0.0014 0.0115 -0.0027 0.0145 0.9364 0.7752 0.0105 0.0013 0.0118 -0.0024 0.0148 0.9232 0.5530 0.0110 0.0013 0.0121 -0.0020 0.0149 0.9269 0.3319 0.0115 0.0013 0.0124 -0.0017 0.0150 0.9485 0.1069 0.0120 0.0012 0.0127 -0.0014 0.0150 0.9899 -0.1276 0.0125 0.0012 0.0130 -0.0010 0.0149 1.0530 -0.3773 0.0130 0.0012 0.0133 -0.0006 0.0147 1.1404 -0.6481 0.0135 0.0012 0.0136 -0.0002 0.0145 1.2544 -0.9455 0.0140 0.0012 0.0139 0.0002 0.0141 1.3967 -1.2743 0.0145 0.0012 0.0142 0.0008 0.0136 1.5677 -1.63680.0150 0.0012 0.0144 0.0013 0.0129 1.7648 -2.0314 0.0155 0.0012 0.0147 0.0020 0.0121 1.9807 -2.4495 0.0160 0.0013 0.0149 0.0027 0.0112 2.2018 -2.8735 0.0165 0.0013 0.0151 0.0035 0.0101 2.4071 -3.2754 0.0170 0.0014 0.0153 0.0044 0.0089 2.5697 -3.6186 0.0175 0.0015 0.0155 0.0053 0.0076 2.6616 -3.8650 0.0180 0.0016 0.0156 0.0063 0.0063 2.6609 -3.9849 0.0185 0.0018 0.0157 0.0072 0.0049 2.5591 -3.9674 0.0190 0.0019 0.0158 0.0080 0.0035 2.3638 -3.8244 0.0195 0.0021 0.0159 0.0088 0.0022 2.0959 -3.5866 0.0200 0.0023 0.0159 0.0095 0.0010 1.7823 -3.2931 0.0205 0.0025 0.0159 0.0100 -0.0001 1.4487 -2.9815 0.0210 0.0027 0.0159 0.0105 -0.0011 1.1152 -2.6809 0.0215 0.0029 0.0159 0.0108 -0.0020 0.7942 -2.4103 0.0220 0.0031 0.0158 0.0111 -0.0028 0.4916 -2.1794 0.0225 0.0033 0.0158 0.0112 -0.0035 0.2086 -1.9913 0.0230 0.0036 0.0157 0.0112 -0.0042 -0.0565 -1.8450 0.0235 0.0038 0.0156 0.0111 -0.0048 -0.3071 -1.7375 0.0240 0.0040 0.0155 0.0110 -0.0054 -0.5475 -1.6650 0.0245 0.0042 0.0154 0.0108 -0.0060 -0.7817 -1.6233 0.0250 0.0044 0.0153 0.0104 -0.0065 -1.0139 -1.6089 0.0255 0.0046 0.0151 0.0100 -0.0071 -1.2479 -1.6181 0.0260 0.0048 0.0150 0.0096 -0.0077 -1.4868 -1.6480 0.0265 0.0050 0.0148 0.0090 -0.0082 -1.7336 -1.6955 0.0270 0.0052 0.0146 0.0084 -0.0088 -1.9905 -1.7574 0.0275 0.0054 0.0145 0.0076 -0.0095 -2.2588 -1.8304 0.0280 0.0055 0.0143 0.0068 -0.0101 -2.5391 -1.9100 0.0285 0.0056 0.0141 0.0058 -0.0108 -2.8305 -1.9910 0.0290 0.0057 0.0138 0.0048 -0.0115 -3.1300 -2.06600.0295 0.0058 0.0136 0.0037 -0.0122 -3.4326 -2.1255 0.0300 0.0059 0.0133 0.0024 -0.0130 -3.7297 -2.1572 0.0305 0.0059 0.0131 0.0011 -0.0137 -4.0091 -2.1451 0.0310 0.0059 0.0128 -0.0004 -0.0145 -4.2538 -2.0696 0.0315 0.0059 0.0125 -0.0019 -0.0152 -4.4419 -1.9079 0.0320 0.0058 0.0122 -0.0035 -0.0158 -4.5473 -1.6352 0.0325 0.0058 0.0119 -0.0051 -0.0163 -4.5411 -1.2273 0.0330 0.0056 0.0115 -0.0066 -0.0166 -4.3954 -0.6661 0.0335 0.0055 0.0112 -0.0081 -0.0167 -4.0889 0.0551 0.0340 0.0053 0.0109 -0.0095 -0.0166 -3.6129 0.9243 0.0345 0.0051 0.0105 -0.0106 -0.0161 -2.9781 1.9058 0.0350 0.0049 0.0102 -0.0115 -0.0152 -2.2178 2.9395 0.0355 0.0047 0.0099 -0.0122 -0.0140 -1.3857 3.9473 0.0360 0.0044 0.0097 -0.0125 -0.0125 -0.5478 4.8458 图形输出：图2 连杆3的⼏个位置点图4 ⾓加速度线图图5 ⾓加速度线图。
Matlab在四杆机构运动分析中的应用MATLAB软件由美国MathWorks公司于1982年推出，经过十几年的发展和竞争，现已成为国际公认的最优秀的科技应用软件之一。
MATLAB提供了强大的矩阵处理和绘图功能。
它主要包括两部分内容：核心函数和工具箱。
Matlab编程代码接近数学推导公式，简洁直观，与科技人员的思维方式和书写习惯相适应，操作简易，人机交互性能好，且可以方便迅速地用三维图形、图像、声音、动画等表达计算结果，拓展思路。
编制相应的M函数文件。
Pos.m用于求解位置、角度和角加速度。
th1为初始角度，w1为杆1角速度，其余为杆长。
将课本P35（i）、（ii）、（iii）式用MATLAB语言表述，即可编制为四杆机构求解函数文件pos.m。
function f=pos(th1,w1,l1,l2,l3,l4)syms x21 x31 x22 x32x1=th1*pi/180;x11=cos(x1);x12=sin(x1);eq1=l1*x11+l2*x21-l3*x31-l4;eq2=l1*x12+l2*x22-l3*x32;eq3=x21^2+x22^2-1;eq4=x31^2+x32^2-1;s=solve(eq1,eq2,eq3,eq4,x21,x22,x31,x32);s1=double(s.x21); s2=double(s.x22);s3=double(s.x31); s4=double(s.x32);x2=(acos(s1(1,1)))/pi*180;x3=(acos(s3(1,1)))/pi*180;A=[l2*s2(1,1),-l3*s4(1,1);l2*s1(1,1),-l3*s3(1,1)];B=[-l1*x12;-l1*x11];w=A\(w1*B);w2=w(1,1);w3=w(2,1);C=[-l2*w(1,1)*s1(1,1),l3*w(2,1)*s3(1,1);l2*w(1,1)*s2(1,1),-l3*w(2,1)*s4(1,1)];D=[w(1,1);w(2,1)];E=[-l1*w1*x11;l1*w1*x12];F=[l2*s2(1,1),-l3*s4(1,1);l2*s1(1,1),-l3*s3(1,1)];t=F\(C*D+w1*E);a2=t(1,1);a3=t(2,1);l1=cat(1,th1,w1,0);l2=cat(1,x2,w2,a2);l3=cat(1,x3,w3,a3);f=(cat(2,l1,l2,l3))';subplot(2,2,1);plot(th1,w2,'r-');hold on;title('连杆2角速度分析');subplot(2,2,2);plot(th1,a2,'b-');hold on;title('连杆2角加速度分析');subplot(2,2,3);plot(th1,w3,'r-');hold on;title('连杆3角速度分析');subplot(2,2,4);plot(th1,a3,'b-');hold on;title('连杆3角加速度分析');return在MA TLAB命令窗口输入命令：>> th1=60;w1=pi/3;l1=20;l2=30;l3=40;l4=45;>> function f=pos(th1,w1,l1,l2,l3,l4)60°，ω=pi/3时的值即可得到theta=采用实时动画的方法编写draw.m文件：function dw = draw(l1,l2,l3,l4,x1,x2,x3)figure(2);th1=x1*pi/180;th2=x2*pi/180;th3=x3*pi/180;x=0:0.001:l4;plot(x,0,'r-');axis([-25,70,-25,60]);hold onfor i=0:0.1:l1;s=i*cos(th1); c=i*sin(th1);plot(s,c,'b-');hold on;end;for x=0:0.1:l3;s2=x*cos(th3)+l4;c2=x*sin(th3);plot(s2,c2);hold on;endfor x=0:0.1:l2;s3=l1*cos(th1)+x*cos(th2); c3=l1*sin(th1)+x*sin(th2);plot(s3,c3);hold on;endreturn代入pos.m中所得的结果，输入命令：>>draw(l1,l2,l3,l4,y(1,1),y(2,1),y(3,1))即可得到当前位置下的四杆图形。
·制造业信息化·Vol．22．,200922009机电产品开发与创新Development ＆Innovation of M achinery ＆E lectrical P roducts3,No．1Jan ．,2010第23卷第1期2010年1月收稿日期：2009－11－06作者简介：李建霞，女，河南信阳人，教授。
主要从事计算机辅助设计教学与研究工作。
已发表核心期刊论文20余篇；王良才，男，山东济南人，副教授。
主要从事机械设计教学与研究工作。
已发表核心期刊论文10余篇。
0引言Matlab 语言是一种非常强大的工程计算语言，提供了非常丰富的Matlab 优化工具箱。
其优化工具箱有许多常用的优化算法，广泛应用于线性规划、二次规划、非线性优化、最小二乘法问题、非线性方程求解、多目标决策等问题，其函数表达简洁，优化算法选择灵活，参数设置自由，相比于其它很多成熟的优化程序具有明显的优越性。
机械优化设计是在现代机械设计理论发展基础上产生的一种新的设计方法，在机械零部件设计中取得了广泛的应用。
机械优化设计是在进行某种机械产品设计时，根据规定的约束条件，优选设计参数，使某项或几项设计指标获得最优值。
在进行优化设计时，首先要建立工程设计问题的数学模型，然后按照数学模型的特点选择优化方法及其计算程序，求得最优设计方案。
1四杆机构优化设计的数学模型四杆机构是机械传动的重要组成部分。
设计四杆机构通常使用作图法与实验法，这两种方法简单易行，但误差较大，运用解析法，能获得所要求的精度，然而用人工进行，整个设计过程是一项繁琐、冗长的工作，甚至可能无法实现。
采用Matlab 优化工具箱对四杆机构进行优化设计，不仅参数输入简单，而且编程工作量小，可更快捷准确地达到设计要求。
四杆机构种类很多，这里只以曲柄摇杆机构为例说明四杆机构优化设计方法。
Matlab 优化工具箱是一套功能强大的工程计算软件，集数值计算、符号运算、可视化建模、仿真和图形处理等多种功能于一体，被广泛应用于机械设计、自动控制和数理统计等工程领域。
matlab四杆机构优化算法四杆机构是机械工程中常用的一种机构，广泛应用于机械传动、机械臂以及其他工程领域。
四杆机构由四个连杆组成，其中两个连杆为主动连杆，另外两个连杆为从动连杆。
四杆机构的运动特性和性能对于机械系统的设计和优化具有重要意义。
本文将介绍一种基于Matlab的四杆机构优化算法，以提高机构的运动性能和效率。
需要明确优化算法的目标。
在四杆机构的优化中，常见的目标包括最小化机构的运动误差、最大化机构的工作范围、最小化机构的能量消耗等。
根据具体应用的需求，可以选择不同的优化目标。
本文以最小化机构的运动误差为例进行讨论。
在四杆机构中，运动误差是指机构在理想情况下应有的运动轨迹与实际运动轨迹之间的差异。
运动误差的大小直接影响到机构的精度和稳定性。
为了最小化运动误差，可以采用以下步骤进行优化。
需要建立四杆机构的运动学模型。
四杆机构的运动学模型描述了机构各个连杆之间的几何关系和运动规律。
根据四杆机构的结构和运动约束条件，可以建立四杆机构的运动学方程。
在Matlab中，可以利用符号运算工具箱求解运动学方程。
需要确定优化的变量和约束条件。
四杆机构的运动性能受到各个连杆的长度、角度等参数的影响。
在优化中，可以将这些参数作为优化的变量，通过调整这些变量来最小化运动误差。
同时，还需要考虑到机构的实际制造和使用条件，设置相应的约束条件，如连杆长度的范围、运动角度的限制等。
接下来，可以采用优化算法对四杆机构进行优化。
常用的优化算法包括遗传算法、粒子群算法、模拟退火算法等。
这些算法可以在给定的变量范围内搜索最优解，通过迭代的方式逐步优化机构的运动性能。
在Matlab中，可以利用优化工具箱中的相应函数进行优化计算。
需要对优化结果进行评估和验证。
优化算法得到的结果可能并不是唯一的最优解，因此需要对优化结果进行评估，判断其是否满足设计要求。
可以通过仿真和实验验证的方式对优化结果进行验证，进一步改进和优化机构的设计。
基于Matlab的四杆机构优化算法可以帮助改善机构的运动性能和效率。
《计算机仿真技术》课程设计报告姓名：冯叶 / 浦合昀学号： 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）如图所示为界面。
