how to solve trigonometric equations
10 次查看(过去 30 天)
显示 更早的评论
Here I have 5 equations:
eq1= DE*cos(alpha)-AC*cos(beta)-CE*cos(sigma)==-AD;
eq2 = DE*sin(alpha)-AC*sin(beta)-CE*sin(sigma)==0;
eq3 = CE*cos(alpha)+CB*cos(beta)-FG*cos(theta)-GE*sin(lambda)-BF*cos(delta)==0;
eq4= CE*sin(alpha)-CB*cos(beta)-FG*sin(theta)+GE*cos(lambda)-BF*sin(delta)==0;
eq5 = sigma - lambda == 1.815;
I've already specified other parameters as follows:
DE = 35.39;
AC = [35.51;40.51];
BF = [20.82;24.52];
CE = 4.31;
GE = 15.29;
FG = 3.5;
AD = 4.75;
CB = 6.5;
Where AC and BF are lengths that are variable with upper and lower bounds.
The only known angle will be beta which is also an input angle with range from 0 deg to 80 deg. I've tried using fsolve but this function can only deal with inputs that are constant not 'range'. Are there any methods to obtain the rest of five unknown angles as results of also a range of angle on the plot?
Something similar to this:
5 个评论
David Goodmanson
2017-8-10
Yes it looks like the results can be pretty strange. What is being varied to give the different plot lines in the graphs? Also, you have a CB*cos(beta) term in both eqn. 3 and eqn. 4. Should one of those be a sin(beta) term?
回答(2 个)
John D'Errico
2017-8-10
You have insufficient information to solve those equations.
There are 5 equations. Unknowns? There are 6 unknowns: {alpha, beta, sigma, lambda, delta, theta}.
So at best there will be infinitely many solutions, an entire manifold of solutions. Worse, two of the coefficients are themselves unknowns, as all you are willing to say is they live in intervals.
Anyone who tells you to set AC and BF to their limits, then solve the problem is simply wrong, because you can never solve the equations for a solution even if AC and BF are fixed and known. And no, interval arithmetic will also fail you here, again because no solution exists.
You can never solve this system of equations as it is. This is not a problem that can be solved EVER using fsolve, nor will solve or vpasolve succeed. Sorry.
Were you to derive a 6th independent equation, independent of the other 5, then you could start to try solving the problem.
6 个评论
John D'Errico
2017-8-10
编辑:John D'Errico
2017-8-10
Yes. There are an infinite number of points. An infinite manifold of points. In fact, there will be more than just a simple surface. There will be a complete volume of points. Beta is NOT known. It ranges over a range that you will allow. But so are two of your coefficients. So you have three parameters that are allowed to range over a THREE dimensionalvolume. Fix those three to ANY number in their ranges, and you will POSSIBLY get some solution. But those three parameters mean that there is a 3-dimensional volume of solutions. That volume will live in the 5 dimensional space of the 5 unknown parameters. It is not as simple as the plot that you show.
My point is that for ANY combination of those three values, you will get 5 other parameters. So you essentially need to plot alpha as a function of (beta,AC,BF). Then plot sigma(beta,AC,BF), theta(beta,AC,BF), etc. But each of these plots is a FOUR dimensional plot.
So, you will need to loop over the possible values of AC, BF, beta. Now you can solve for the other 5 parameters. For ANY value of beta, compute the min and max of alpha. Do this for all AC and BF in their ranges. That gives you two curves, each as a function of beta. It does not give you the complete understanding of that volume. But it will give you the 2-d plot you are looking for.
John BG
2017-8-14
Hi Jenny
1.
The equations in page 409, 3rd page, of the document you have supplied, just above Figure 4, are not complete.
The authors supply the main equations, but there are missing boundary equations needed to narrow down results, for instance the rigid (constant) angle µ.
The key comment is
.
' .. a MATLAB model is created ..'
.
The equations (2) to (5) seem to be used in the model, but fur such unknown model to supply graphs 5 to 7, it has to take more equations into account that (2) to (5).
Probably a MATLAB SIMSCAPE model has been built, or with Autodesk, and Figure 4
.
% 002
.
shows precisely a couple views of the simulated hardware of model, which prodces Figures 5 to 7.
If one attempts to work just with supplied equations only, as already commented by different readers, the direct assault doesn't work.
2.
One can then try combining equations, like (2) and (3)
with AC spring extended (ABC=AC=40.51) and BF spring compressed (BF=20.82)
close all; clear all; clc
DE = 35.39;
% AC = [35.51;40.51]; % spring compressed 35.51 spring extended 40.51
AC=40.51;
AB=35.51;
% BF = [20.82;24.52]; % again, spring
BF=20.82;
CE = 4.31;
GE = 15.29;
FG = 3.5;
AD = 4.75;
CB = 6.5;
% AB=35;DE=35;CE=5;AC=40;AD=5;
da=1/10*pi/180*2;a=[0:da:pi/2];
b=a;
[A,B]=meshgrid(a,b);
Y=acos(-1/CE*(AC*cos(B)-AD-DE*cos(A)))-asin(1/CE*(DE*sin(A)-AC*sin(B)));
Yimag=imag(Y);
hs=surf(A,B,Yimag)
hs.EdgeColor='none'
campos([.8 .8 5])
.
% 101
.
Now with AC spring compressed (ABC=AB=AC=35.51) and BF spring extended (BF=24.52)
.
close all; clear all; clc
DE = 35.39;
% AC = [35.51;40.51]; % spring compressed 35.51 spring extended 40.51
AC=35.51;
AB=35.51;
% BF = [20.82;24.52]; % again, spring
BF=24.52;
CE = 4.31;
GE = 15.29;
FG = 3.5;
AD = 4.75;
CB = 6.5;
% AB=35;DE=35;CE=5;AC=40;AD=5;
da=1/10*pi/180*2;a=[0:da:pi/2];
b=a;
[A,B]=meshgrid(a,b);
Y=acos(-1/CE*(AC*cos(B)-AD-DE*cos(A)))-asin(1/CE*(DE*sin(A)-AC*sin(B)));
Yimag=imag(Y);
hs=surf(A,B,Yimag)
hs.EdgeColor='none'
campos([.8 .8 5])
.
% 102
.
One notes that combining (2) and (3), although there's a clear narrow diagonal strip of null imaginary values of the resulting angle sigma (just considering levers ABC and DE, joined with CE), such angles are possible, but the complete model has a far more narrower imag(sigma)=0 gap between alpha and beta, really close but not null, can be appreciated around mid range graph Figure 5.
.
% 003
.
3.
So, if in order to generate the graphs you want to build the model with Simscape https://uk.mathworks.com/products/simscape.html
a start point could from scratch, with an empty Simscape model.
Or open already existing Simscape examples like
4.
There are people who have attempted to model similar moving linkage systems direrctly with MATLAB only, like
attached function wattLinkage.m
attempting to modify these kind of examples can be difficult because there are no blocks to play with, all the equations have to be available.
5.
Another example, 4 linkages, with similar approach
This example, attached script Four_bar_linkage.m has the same configuration (horizontal) of the section ABC and DE.
6.
Yet another of interest
Attached Linkage_mechanism.zip
the test lines
xxx = mnism(3,5,6,6,8)
xxx = animism(3,5,6,7,8,1,10)
show again a similar configuration to a partial of the sought model.
.
% 103
.
7.
Jed's Linkage solver helps build the equations
Linkage_solver.m
8.
Another option is Do It Yourself:
The linkage diagram to reproduce, animate and then read the angles of interest, to obtain the 3 graphs, is the following
.
% 001
.
The 1st section, related to Figure 5 graph, with 4 sample points
clear all;close all;clc;
% format long;
% stepper motor angle = input angle = beta = b: angle between beam ABC and vertical knucle plane
%
% angles that tell whether good grip can be achieved:
% Proximal Phalanx angle = alpha = a
% angle MP = beta in diagram = b : metacarpophalangeal or 1st phalanx. Spelling error in text, the word metacarphopahalangeal
doesn't exist, page 1
% PIP = sigma in diagram = s: proximal phalanx or 2nd phaanx
% DIP = delta in diagram = d: distal phalanx or 3rd phalanx
% alpha, a: angle
AB=35;DE=35;CE=5;AC=40;AD=5; % beams, test values, slightly simplified from supplied figures
dv=.5;v1=[0:dv:180];circ1=[sind(v1);cosd(v1)]; % r=1 centre [0 0] reference circle points [x1 x2 x3 ..; y1 y2 y3 ..]
v2=[0:dv:360];circ2=[sind(v2);cosd(v2)];
cntrABC=[0;50];cntrDE=[0;50+AD]; % static pivot points on vertical knucle plane
% closest point that point E can get to knucle plane, with cosines trigonometric law, checked it matches calculating beta swapping
beams not shown
% S1/sin(a1)=S2/sin(a2)=S3/sin(a3)
% A1^2=S2^2+S3^2-2*S2*S3*cos(a1)
% translated to the diagram: S2=S3=DE S1=CE a1=b0
a0=acosd(-1/(2*DE^2)*(CE^2-2*DE^2))
cntrCE1=DE*[sind(a0);cosd(a0)]+cntrDE % [4.9872 89.6429]
cntrCE2=[27.77;76.31];cntrCE3=[34.78;58.96];cntrCE4=[25.6;31.13];cntrCE5=[0;20]; % sample points along circle DE radius
circAB=AB*circ1+cntrABC;
circAC=AC*circ1+cntrABC;
circDE=DE*circ1+cntrDE;
circCE1=CE*circ2+cntrCE1; % sample circles along circle BE radius
circCE2=CE*circ2+cntrCE2;
circCE3=CE*circ2+cntrCE3;
circCE4=CE*circ2+cntrCE4;
circCE5=CE*circ2+cntrCE5;
hf=figure;ax=gca
plot(zeros(1,2),[0 120],'k--','LineWidth',1) % knucle or vertical plane
hold all
plot(circAB(1,:),circAB(2,:),'Color',[.1 .1 .8],'LineWidth',1) % beam ABC circle
plot(circAC(1,:),circAC(2,:),'Color',[.1 .1 .8],'LineWidth',1)
plot(circDE(1,:),circDE(2,:),'Color',[.8 .1 .1],'LineWidth',1) % beam DE circle
plot([cntrDE(1) cntrCE1(1)],[cntrDE(2) cntrCE1(2)],'Color',[.8 .5 .5],'LineWidth',4) % beam DE, sample position
plot([cntrDE(1) cntrCE2(1)],[cntrDE(2) cntrCE2(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot([cntrDE(1) cntrCE3(1)],[cntrDE(2) cntrCE3(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot([cntrDE(1) cntrCE4(1)],[cntrDE(2) cntrCE4(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot([cntrDE(1) cntrCE5(1)],[cntrDE(2) cntrCE5(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot(circCE1(1,:),circCE1(2,:),'Color',[.8 .5 .5],'LineWidth',1) % sample circles
plot(circCE2(1,:),circCE2(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(circCE3(1,:),circCE3(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(circCE4(1,:),circCE4(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(circCE5(1,:),circCE5(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(cntrABC(1),cntrABC(2),'bo') % beam ABC circle centre
plot(cntrDE(1),cntrDE(2),'ro') % beam DE circle centre
text(cntrABC(1)-5,cntrABC(2),'A') % text point A
text(cntrDE(1)-5,cntrDE(2),'D') % text point D
plot(cntrCE1(1),cntrCE1(2),'ro') % sample circles CE centres
plot(cntrCE2(1),cntrCE2(2),'ro')
plot(cntrCE3(1),cntrCE3(2),'ro')
plot(cntrCE4(1),cntrCE4(2),'ro')
plot(cntrCE5(1),cntrCE5(2),'ro')
% the spring pushes out, expands AC towards AD, not AD backwards to AC
% 1st sample position beta=0 then alpha is not zero
[x_CE1_AC,y_CE1_AC]=intersections(circCE1(1,:),circCE1(2,:),circAC(1,:),circAC(2,:)); % when motor input angle b=0, a is already a
few degree ahead, the Proximal Phalanx angle in 1st graph, Y axis in 1st graph, effect of the spring pushing. a passive device like a
spring substantially limits the grimping capacity, thus the top weight that can be gripped.
[xmin nxmin]=min(x_CE1_AC); % choose intersection point closer to knucle plane
x_CE1_AC=x_CE1_AC(nxmin);y_CE1_AC=y_CE1_AC(nxmin);
plot(x_CE1_AC,y_CE1_AC,'bd') % joint CE AC
plot(cntrCE1(1),cntrCE1(2),'go') % join CE DE
plot([x_CE1_AC cntrCE1(1)],[y_CE1_AC cntrCE1(2)],'Color',[.8 .5 .5],'LineWidth',4) % beam CE
plot([cntrABC(1) x_CE1_AC],[cntrABC(2) y_CE1_AC],'Color',[.5 .5 .8],'LineWidth',4) % beam ABC, 1st sample position
% when beta=0 0, alpha=MP angle is not zero.
b1=0
a1=a0 % angle when all up
% 2nd sample position
[x_CE2_AC,y_CE2_AC]=intersections(circCE2(1,:),circCE2(2,:),circAC(1,:),circAC(2,:));
[xmin nxmin]=min(x_CE2_AC); % choose intersection point closer to knucle plane
x_CE2_AC=x_CE2_AC(nxmin);y_CE2_AC=y_CE2_AC(nxmin);
plot(x_CE2_AC,y_CE2_AC,'bd') % joint CE AC
plot(cntrCE2(1),cntrCE2(2),'go') % joint CE DE
plot([x_CE2_AC cntrCE2(1)],[y_CE2_AC cntrCE2(2)],'Color',[.8 .5 .5],'LineWidth',4) % beam CE
plot([cntrABC(1) x_CE2_AC],[cntrABC(2) y_CE2_AC],'Color',[.5 .5 .8],'LineWidth',4) % beam ABC, 1st sample position
b2=atand((x_CE2_AC-cntrABC(1))/(y_CE2_AC-cntrABC(2)))
a2=atand((cntrCE2(1)-cntrDE(1))/(cntrCE2(2)-cntrDE(2)))
b2-(a2+a0) % just shy of half degree different alpha-beta, when beta ~ 45° matching green blue graph lines in 1st graph page 4
% =
% -0.4198
% 3rd sample position
[x_CE3_AC,y_CE3_AC]=intersections(circCE3(1,:),circCE3(2,:),circAC(1,:),circAC(2,:));
[xmin nxmin]=min(x_CE3_AC); % choose intersection point closer to knucle plane
x_CE3_AC=x_CE3_AC(nxmin);y_CE3_AC=y_CE3_AC(nxmin);
plot(x_CE3_AC,y_CE3_AC,'bd') % joint CE AC
plot(cntrCE3(1),cntrCE3(2),'go') % joint CE DE
plot([x_CE3_AC cntrCE3(1)],[y_CE3_AC cntrCE3(2)],'Color',[.8 .5 .5],'LineWidth',4) % beam CE
plot([cntrABC(1) x_CE3_AC],[cntrABC(2) y_CE3_AC],'Color',[.5 .5 .8],'LineWidth',4) % beam ABC, 1st sample position
b3=atand((x_CE3_AC-cntrABC(1))/(y_CE3_AC-cntrABC(2)))
a3=atand((cntrCE3(1)-cntrDE(1))/(cntrCE3(2)-cntrDE(2)))
b3-(a3+a0)
hf.Position=[200 200 1000 600]
ax.DataAspectRatio=[1 1 1]
ax.XLim=[-10 150];
ax.YLim=[0 100];
ax.PlotBoxAspectRatio=[1 1 1]
.
% 010
.
gripping with beta>90° not plotted in document, it may either damage some fragile objects or the spring may reach full run and not pushing so no grip, too small items may drop
.
% MIP
% EG=15 but because in linkage diagram, for visual convenience EG length
% similar to DE length, I am choosing EG=DE
EG=DE
% BF=[20 24] % again BF has another spring, up to 24
BF1=40
%FG=4 % short beam, again from linkage diagram, despite FG value < CE, in diagram look same size so FG=CE
FG=CE
% rigid angle u not supplied, but from linkage diagram, counting pixels,
% triangle EQG complementary angle to u is ~ 17° thus u=107°
u=109
% [x_CE1_AC,y_CE1_AC] rotate CW around [cntrCE1(1),cntrCE1(2)]
G01=[cosd(-u) -sind(-u);sind(-u) cosd(-u)]*[x_CE1_AC-cntrCE1(1);y_CE1_AC-cntrCE1(2)]+[cntrCE1(1);cntrCE1(2)]
plot(G01(1),G01(2),'go')
m1=(G01(2)-cntrCE1(2))/(G01(1)-cntrCE1(1))
G1=[cntrCE1(1)+EG*cosd(atand(m1));cntrCE1(2)+EG*sind(atand(m1))]
plot([cntrCE1(1) G1(1)],[cntrCE1(2) G1(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot(G1(1),G1(2),'ro') % 1st sample point G acquired
% BF also has another spring pushing forward, but there's no spring
% compression point like B on ABC, assuming BF fully extended
% [x_CE2_AC,y_CE2_AC] rotate CW around [cntrCE2(1),cntrCE2(2)]
G02=[cosd(-u) -sind(-u);sind(-u) cosd(-u)]*[x_CE2_AC-cntrCE2(1);y_CE2_AC-cntrCE2(2)]+[cntrCE2(1);cntrCE2(2)]
plot(G02(1),G02(2),'go')
m1=(G02(2)-cntrCE2(2))/(G02(1)-cntrCE2(1))
G2=[cntrCE2(1)+EG*cosd(atand(m1));cntrCE2(2)+EG*sind(atand(m1))]
plot([cntrCE2(1) G2(1)],[cntrCE2(2) G2(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot(G2(1),G2(2),'ro')
% [x_CE3_AC,y_CE3_AC] rotate CW around [cntrCE3(1),cntrCE3(2)]
G03=[cosd(-u) -sind(-u);sind(-u) cosd(-u)]*[x_CE3_AC-cntrCE3(1);y_CE3_AC-cntrCE3(2)]+[cntrCE3(1);cntrCE3(2)]
plot(G03(1),G03(2),'go')
m1=(G03(2)-cntrCE3(2))/(G03(1)-cntrCE3(1))
G3=[cntrCE3(1)+EG*cosd(atand(m1));cntrCE3(2)+EG*sind(atand(m1))]
plot([cntrCE3(1) G3(1)],[cntrCE3(2) G3(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot(G3(1),G3(2),'ro')
ax.YLim=[0 140];
.
% 011
.
% calculate 1st set of sigmas s, mCE: slope section between points C and E
% reviewing the diagram: if m_angle<0 s=90°+abs(m_angle); if m_angle>0 s=90-m_angle+180
a1CE=atand((cntrCE1(2)-y_CE1_AC)/(cntrCE1(1)-x_CE1_AC))
if a1CE<0 && abs(a1CE)<90 s1=90+abs(a1CE); end
if a1CE<0 && abs(a1CE)==90 s1=90; end
if a1CE>0 s1=180+(90-abs(a1CE)); end
a2CE=atand((cntrCE2(2)-y_CE2_AC)/(cntrCE2(1)-x_CE2_AC))
if a2CE<0 && abs(a2CE)<90 s2=90+abs(a2CE); end
if a2CE<0 && abs(a2CE)==90 s2=90; end
if a2CE>0 s2=180+(90-abs(a2CE)); end
a3CE=atand((cntrCE3(2)-y_CE3_AC)/(cntrCE3(1)-x_CE3_AC))
if a3CE<0 && abs(a3CE)<90 s3=90+abs(a3CE); end
if a3CE<0 && abs(a3CE)==90 s3=90; end
if a3CE>0 s3=180+(90-abs(a3CE)); end
s1
s2
s3
s1 =
94.0961
s2 =
153.5055
s3 =
218.0812
% now going for 1st set of F points
circFG1=FG*circ2+G1; % small beam circles
circFG2=FG*circ2+G2;
circFG3=FG*circ2+G3;
plot(circFG1(1,:),circFG1(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(circFG2(1,:),circFG2(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(circFG3(1,:),circFG3(2,:),'Color',[.8 .5 .5],'LineWidth',1)
v3=[0:dv:45];v4=[45:dv:130];v5=[135:dv:170]; % avoiding overcrowding schematic, now just arcs
circ3=[sind(v3);cosd(v3)]; % support arcs
circ4=[sind(v4);cosd(v4)];
circ5=[sind(v5);cosd(v5)];
B1=AB*[sind(b1);cosd(b1)]+cntrABC % sample B points
plot(B1(1),B1(2),'b*')
B2=AB*[sind(b2);cosd(b2)]+cntrABC
plot(B2(1),B2(2),'b*')
B3=AB*[sind(b3);cosd(b3)]+cntrABC
plot(B3(1),B3(2),'b*')
circBF1=BF1*circ3+B1 % support circles
circBF2=BF1*circ4+B2
circBF3=BF1*circ5+B3
plot(circBF1(1,:),circBF1(2,:),'Color',[.7 .7 .7])
plot(circBF2(1,:),circBF2(2,:),'Color',[.7 .7 .7])
plot(circBF3(1,:),circBF3(2,:),'Color',[.7 .7 .7])
[x_FG1_BF,y_FG1_BF]=intersections(circBF1(1,:),circBF1(2,:),circFG1(1,:),circFG1(2,:));
[ymax nymax]=max(y_FG1_BF); % choose intersection point with higher y
x_FG1_BF=x_FG1_BF(nymax);y_FG1_BF=y_FG1_BF(nymax); % F point acquired
plot(x_FG1_BF,y_FG1_BF,'bd') % joint FG BF
plot([x_FG1_BF B1(1)],[y_FG1_BF B1(2)],'Color',[.5 .5 .8],'LineWidth',4) % beam BF
.
% 012
.
I have exhausted the time allocated to your question without feed back from you yet.
The above lines attached test_123.m
Jenny would you please comment how would you like to proceed in order to obtain graphs 5 to 7?
Do you take it from here, or would you like to see how to complete the points of interest, assess angles, and compare to graph?
Therefore
if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?
To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link
thanks in advance
John BG
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Bodies 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!