用matlab求解方程组
l1sinθ1+l2sinθ2=l3sinθ3+l4sinθ4
l1cosθ1+l2cosθ2=l3cosθ3+l4cosθ4
θ1=ω1t
其中θ为角度，要求θ4的角速度ω4。
l1=150mm l2=34.25mm
l3=131.2mm l4=138.35mm θ20=0 θ10=0 ω1=100*pi/3 rad/s 展开
以前曾回答过相关的问题（编号2010907084650078548），虽然给了一些求解的方法，但结果始终不是很清晰，不能直观地解释为什么换了数据就不行。后来，我换了一种思路，应该可以很直观地说明不同的参数有的可以而有的不行，现简单介绍如下。
基本思想就是，既然要求解的是关于Q3、Q4的两个方程构成的方程组，那么，我们把两个方程的曲线都用隐函数绘图（ezplot）画出来，看看曲线的交点是怎样随时间t变化的。
对于第一组参数
1l1=150; l2=34.25; l3=131.2; l4=138.35;
画出的图形如下。由图可见，在一个完整的周期内（Q3、Q4均变化2*pi），两个方程的曲线始终有交点，也就是解始终存在，所以，这组参数没问题。
但换了另外一组参数后：
1l1=60; l2=50.1799; l3=38.1013; l4=100.4134;
我们看到，随着t变化，在某些条件下，两个方程的曲线不再相交，也就意味着此时方程组无解。至于说物理意义，我觉得可能类似于连杆的臂长在工作过程中不能达到要求，所以无法正常工作。
东莞大凡
2024-11-14 广告
题主几天前问到这个问题（编号1511389540175596140），我做了一些研究，但因为感觉还没有真正搞明白，所以没有及时回答。因时间精力所限，可能暂时没法再深入研究了，所以，姑且把已知的一些情况贴出来，供有意研究者参考。
1、根据题主在之前提问时所说，这个是【行星齿轮连杆机构的矢量方程】。我对这个方面的了解有限，但凭直觉猜测，在相应的学科领域可能会有一些更好的数学技巧来解决这类问题，如果确实有，再辅以MATLAB，也许会有更好的结果。
2、我现在采用的做法是，使用符号运算工具箱的solve函数求解。
需要说明一下：不同版本的MATLAB所用的符号运算内核可能存在差别，解方程得到的结果也可能会差别很大。尤其值得一提的是，MATLAB在2008a之前使用的是Maple内核，而2008b之后换用MuPad内核。很多人（包括我在内）认为Maple在多数场合的表现要优于MuPad，但到现在毕竟已经过了很多年，MuPad后来逐渐在很多方面可能要比Maple更好（是指2008a之前的Maple内核，毕竟Maple本身也是一个独立的软件，这些年也在发展，但我没有进行更多对比）。
我个人使用的几台机器中，安装了多个版本的MATLAB，对于同一个问题所写的代码，我通常会在多个不同的环境下测试并加以分析改进，以尽量减少因软件版本差异导致的结果不确定。
3、就题主这个问题而言，可以大致编写如下代码：
l1=150; l2=34.25;
l3=131.2; l4=138.35; Q20=0; Q10=0; w1=100*pi/3;
syms t
Q1=w1*t;
Q2=Q20-(Q1-Q10);
syms Q3 Q4
eq1=l1*sin(Q1)+l2*sin(Q2)-(l3*sin(Q3)+l4*sin(Q4));
eq2=l1*cos(Q1)+l2*cos(Q2)-(l3*cos(Q3)+l4*cos(Q4));
Q = solve(eq1,eq2,Q3,Q4);
T = linspace(0,0.1,50);
n=0;
for i=1:length(Q.Q4)
if any(simple(subs([eq1 eq2],{Q3 Q4},{Q.Q3(i) Q.Q4(i)})) ~= 0),
continue,
end
n = n + 1;
figure(n)
subplot(2,1,1);
plot(T,double(subs(Q.Q4(i),T)))
xlabel \itt; ylabel {\it\theta}_4
subplot(2,1,2);
plot(T,double(subs(diff(Q.Q4(i)),T)))
xlabel \itt; ylabel {\it\omega}_4
end
代码本身没有太多需要说明的地方，基本思路就是，先定义时间t为符号变量，然后把方程表达出来，再用solve求解，最后把t的值代入进行绘图。用solve求解方程组得到的是θ3、θ4，将其对时间求导数即可得到ω3、ω4。
然而，这段代码在不同版本的MATLAB中的运行情况差别很大。我最初是在2014a上编写的，得到的结果如下（图形是用2014b生成的，除外观稍有差别外，曲线是一致的）：
对于这个结果，我认为基本是可信的。值得说明的有两点：
（1）用solve解方程，可得到多个解（确切地说，是4个），但部分（其中的2个）解代回原方程之后却并不成立（疑问：不知道这种解是怎样求出来的），所以在代码的for循环中有个if语句对此进行判断。
（2）θ曲线的周期性跳变是由于atan函数的值域引起的，实际情况应该是连续的多周转动。
然而，把这段代码放在其它一些版本中运行时，又遇到了一些意想不到的情况：在2007b和2008a中得到的结果如下（结果相同，但2008a的速度似乎快一些——有可能是其它因素导致，因为两个版本在不同机器上，硬件和操作系统都不同）：
而在6.5版中，solve求出的只有两个解并且都是满足方程的，也就是说，其实可以不需要代回原方程进行判断。结果如下：
至于在2008b中（这是第一个正式使用MuPad内核的版本），情况比较惨——它干脆解方程失败（Explicit solution could not be found），后续的代码当然也就无法执行了。
总结
我编写的代码在6.5、2007b、2008a、2008b、2014a、2014b等版本中做了测试，至少见到了四种不同的结果。至于哪种更可信，我个人更倾向于相信2014a/b的结果，但更深入的原因我说不清，其它版本为什么会出现那样的结果？这个留待其他高手来解答吧。
我目前能给的建议就是，符号运算在很多时候非常方便，但毕竟只是一个工具，不要过分迷信，使用时需要保持一份谨慎，可能的话换个版本多做一做测试，遇到问题再具体分析。
theta20=0;theta10=0;omega1=100*pi/3;
syms theta2 theta1 theta3 theta4 t
eq1=theta2-theta20-(theta1-theta10);
eq2=l1*sin(theta1)+l2*sin(theta2)-l3*sin(theta3)-l4*sin(theta4);
eq3=l1*cos(theta1)+l2*cos(theta2)-l3*cos(theta3)-l4*cos(theta4);
eq4=theta1-omega1*t;
s=solve(eq1,eq2,eq3,eq4,theta1,theta2,theta3,theta4);
omega4=s.theta4/t
