2 views (last 30 days)
gheorghe on 8 Dec 2012
I am trying to calculate the acceleration of a mechanism and I am stuck with this error:
here are the code lines that describe the acceleration and the aD42 variable:
aE=[0 0 0];
ep5z = sym('ep5z','real');
ep5 = [ 0 0 ep5z ];
Ep5 = [0 0 ep5zs];
dot(Omega5,Omega5)*(rD-rE);
##### 2 CommentsShowHide 1 older comment
Walter Roberson on 8 Dec 2012
Put a breakpoint after the solve() and see what fields were returned for solaD
Note: you should not eval() a symbolic variable. symbolic variables are not MATLAB expressions, they are just very similar to MATLAB expressions. You should subs() or double() the symbolic variable.

Matt Fig on 8 Dec 2012
Edited: Matt Fig on 8 Dec 2012
Use this line:
Instead of the one you use for solaD. The problem is you have three unknowns and two equations, so SOLVE chose to leave aD42 as the independent variable in the solutions.
Also, when you call FPRINTF with a symbolic variable you need to convert in with CHAR first.
and so for the rest.
gheorghe on 8 Dec 2012
Thank you so very much!

gheorghe on 8 Dec 2012
clear all % sterge toate variabilele
clc % sterge fereastra command si duce cursorul in prima pozitie
close all % inchide toate ferestrele grafice
AB=0.25; % marimea segmentului AB in m
BC=223.94*0.0025;
phi=130*(pi/180); % definirea unghiului phi1
%pozitia cuplei A
xA=0; yA=0;
%pozitia cuplei B
xB=AB*cos(phi);
yB=AB*sin(phi);
%pozitia cuplei C
yC=0;
eqnD1='(xB-xCsol)^2+(yB-yC)^2=BC^2'; % distanta intre 2 pct
solC=solve(eqnD1,'xCsol'); %rezolvarea ec pt xC
xC1=eval(solC(1));
xC2=eval(solC(2));
if xC1 > xB xC = xC1;
else xC = xC2; end %conditia ca cupla sa se afle in cadranul 1
phi2 = atan((yB-yC)/(xB-xC));%unghiul care il face bara BC cu axa xx
% vectori pozitie in coordonate spatiale x y z
rA=[xA yA 0];
rB=[xB yB 0];
rC=[xC yC 0];
AE=0.3;
BD=90.42*0.0025;
ED=149.76*0.0025;
CF=357.73*0.0025; % din C pana in capatul barei
DF=43.37*0.0025;
phi0=pi;
xE=0;
yE=0.3;
% punctele B C D sunt pe aceeasi linie au acceasi panta
% m = (y2 -y1) / (x2 -x1)
eqnD1='(xB-xDsol)^2+(yB-yDsol)^2=BD^2'; % distanta intre 2 puncte
eqnD2=' (yB-yC) / (xB -xC)= (yDsol - yB) / (xDsol -xB)';
solD = solve(eqnD1, eqnD2, 'xDsol, yDsol');
xDpositions = eval(solD.xDsol);
yFpositions = eval(solD.yDsol);
% first component of the vector xDpositions
xD1 = xDpositions(1);
% second component of the vector xDpositions
xD2 = xDpositions(2);
% first component of the vector yDpositions
yD1 = yFpositions(1);
% second component of the vector yDpositions
yD2 = yFpositions(2);
if xD1 <= xC
xD = xD1; yD=yD1;
else
xD = xD2; yD=yD2;
end
rD = [xD yD 0]; % vector pozitie D
eqnF1='(xD-xFsol)^2+(yD-yFsol)^2=DF^2'; % distanta intre 2 puncte
eqnF2=' (yB-yC) / (xB -xC)= (yFsol - yB) / (xFsol -xB)';
solF = solve(eqnF1, eqnF2, 'xFsol, yFsol');
xFpositions = eval(solF.xFsol);
yFpositions = eval(solF.yFsol);
% first component of the vector xDpositions
xF1 = xFpositions(1);
% second component of the vector xDpositions
xF2 = xFpositions(2);
% first component of the vector yDpositions
yF1 = yFpositions(1);
% second component of the vector yDpositions
yF2 = yFpositions(2);
if xF1 <= xC
xF = xF1; yF=yF1;
else
xF = xF2; yF=yF2;
end
%sau mai simplu xF = xD + DF*cos(phi2+pi) ; yF = yD + DF*sin(phi2+pi) ;
rF = [xF yF 0]; % pozitie vector F
%trasarea graficului
plot([xA,xB],[yA,yB],'r-o',[xB,xC],[yB,yC],'b-o',[xB,xD],[yB,yD],'g-o',...
[xE,xD],[yE,yD],'y-o',[xD,xF],[yD,yF],'k-o') % trasarea mecanismului
%cu proprietaile axis([xMIN xMAX yMIN yMAX]); r = red; - = solid line;
%o = cerculet in fiecare punct; b = blue; g = green; c = cyan; m = magenta
xlabel('x (m)') % definirea axelor
ylabel('y (m)')
title('pozitie mecanism pentru \phi = 130 (deg)') % titlu grafic
text(xA,yA,' A'),... % virgula si cele 3 puncte reprezinta asteptarea
% urmatoarei comenzi fara a se executa instructiunea
text(xB,yB,' B'),...
text(xC,yC,' C'),...
text(xD,yD,' D'),...
text(xE,yE,' E'),...
text(xE,yE,' E'),...
axis([-0.6 0.6 -0.2 0.5]),... % marimea axelor
grid
% mai avem comenzi ca hold on - retine graficul
% ’LineWidth’,1.5 grosime linie
% grid on
% text(xA,yA,'\leftarrow A = ground',...'HorizontalAlignment','left'),...
phi4 = phi2;
phi5 = atan((yD-yE)/(xD-xE));
rE = [xE yE 0]; % pozitie vector E
%afisarea valorilor in fereastra command
fprintf('xA = %g (m) \n', xA)
fprintf('yA = %g (m) \n', yA)
fprintf('xB = %g (m) \n', xB)
fprintf('yB = %g (m) \n', yB)
fprintf('xC = %g (m) \n', xC)
fprintf('yC = %g (m) \n', yC)
fprintf('xD = %g (m) \n', xD)
fprintf('yD = %g (m) \n', yD)
fprintf('xE = %g (m) \n', xE)
fprintf('yE = %g (m) \n', yE)
fprintf('xF = %g (m) \n', xF)
fprintf('yF = %g (m) \n', yF)
fprintf('phi0 = %g (degrees) \n', phi0*180/pi)
fprintf('phi1 = %g (degrees) \n', phi*180/pi)
fprintf('phi2 = %g (degrees) \n', phi2*180/pi)
fprintf('phi4 = %g (degrees) \n', phi4*180/pi)
fprintf('phi5 = %g (degrees) \n', phi5*180/pi)
fprintf('rA = [ %g, %g, %g ] (m) \n', rA)
fprintf('rB = [ %g, %g, %g ] (m) \n', rB)
fprintf('rC = [ %g, %g, %g ] (m) \n', rC)
fprintf('rD = [ %g, %g, %g ] (m) \n', rD)
fprintf('rE = [ %g, %g, %g ] (m) \n', rE)
fprintf('rF = [ %g, %g, %g ] (m) \n', rF)
% vB = vB1 = vA+vBA = vA+?1×rAB = ?1×rB
omega1=[0 0 20]; % rad/s constant
vA=[0 0 0]; %punctul A se afla pe batiu
vB1=vA+cross(omega1,rB); % viteza punctului B1
vB2=vB1;
vB = norm(vB1); % norm() is the vector norm
fprintf('omega1 = [ %g, %g, %g ] (rad/s)\n', omega1)
fprintf('vB=vB1=vB2 = [ %g, %g, %g ] (m/s)\n', vB1)
fprintf('|vB|= %g (m/s)\n', vB)
% dot(u,v) produs scalar dintre u si v
% cross(u,v) produs vectorial dintre u si v
% vC = vC2 = vB+vCB = vB2 +?2×rBC = vB+?2×(rC ?rB)
omega2z = sym('omega2z','real');% creare variabile simbolice sym
vCx = sym('vCx','real');
omega2 = [ 0 0 omega2z ];
vC = [ vCx 0 0 ];
eqvC = vC - (vB2 + cross(omega2,rC-rB));
eqvCx = eqvC(1); % equation component on x-axis
eqvCy = eqvC(2); % equation component on y-axis
solvC = solve(eqvCx,eqvCy);
% solutiile
omega2zs = eval(solvC.omega2z);
vCxs = eval(solvC.vCx);
% rezulta
Omega2 = [0 0 omega2zs];
VC = [vCxs 0 0];
vCB = cross(Omega2,rC-rB); %viteza lui C in functie de B
% pentru a arata corect expresiile ecuatiilor eqvCx si y
qvCx = vpa(eqvCx,6);
fprintf('x-axis: %s = 0 \n', char(qvCx))
qvCy = vpa(eqvCy,6);
fprintf('y-axis: %s = 0 \n', char(qvCy))
%vpa(S,D) compute each element of S to D decimal digits of accuracy
% char() creates a character array (string)
% fprintf('|vC|= %g (m/s)\n', norm(VC))
fprintf('=>\n')
fprintf('vCx = %g (m/s)\n', vCxs)
fprintf('\n')
fprintf('omega2 = [ %g, %g, %g ] (rad/s)\n', Omega2)
fprintf('vC = [ %g, %g, %g ] (m/s)\n', VC)
fprintf('vCB = [ %g, %g, %d ] (m/s)\n', vCB)
fprintf('\n')
% Acceleratie punct B
aA = [0 0 0 ]; % (m/sˆ2)
ep1=[0 0 0]; % rad/s^2 omega1 este constant
aB1 = aA + cross(ep1,rB) - dot(omega1,omega1)*rB;
aB2 = aB1;
aBn = - dot(omega1,omega1)*rB;
aBt = cross(ep1,rB);
fprintf('|aB1|=|aB2|= %g (m/s^2)\n', norm(aB1))
fprintf('|aBn|= %g (m/s^2)\n', norm(aBn))
fprintf('|aBt|= %g (m/s^2)\n', norm(aBt))
ep2z=sym('ep2z','real');
aCx=sym('aCx','real');
ep2 = [ 0 0 ep2z ];
aC = [aCx 0 0 ];
eqaC=aC-(aB1+cross(ep2,rC-rB)-...
dot(Omega2,Omega2)*(rC-rB));
eqaCx = eqaC(1); % equation component on x-axis
eqaCy = eqaC(2); % equation component on y-axis
solaC = solve(eqaCx,eqaCy);
ep2zs=eval(solaC.ep2z);
aCxs=eval(solaC.aCx);
Ep2 = [0 0 ep2zs];
aCs = [aCxs 0 0];
aCB=cross(Ep2,rC-rB)-dot(Omega2,Omega2)*(rC-rB);
aCBn=-dot(Omega2,Omega2)*(rC-rB);
aCBt=cross(Ep2,rC-rB);
fprintf('aC = aB + omega2 x rBC - (omega2.omega2)rBC =>\n')
qaCx=vpa(eqaCx,6);
fprintf('x-axis: %s = 0 \n', char(qaCx))
qaCy=vpa(eqaCy,6);
fprintf('y-axis: %s = 0 \n', char(qaCy))
fprintf('=>\n')
fprintf('aCx = %g (m/sˆ2)\n', aCxs)
fprintf('\n')
fprintf('alpha2 = [%g, %g, %g] (rad/sˆ2)\n', Ep2)
fprintf('aC = [ %g, %g, %g ] (m/sˆ2)\n', aCs)
fprintf('aCB = [ %g, %g, %d ] (m/sˆ2)\n', aCB)
fprintf('aCBn = [ %g, %g, %d ] (m/sˆ2)\n', aCBn)
fprintf('|aCBn| = %g (m/sˆ2)\n', norm(aCBn))
fprintf('aCBt = [ %g, %g, %d ] (m/sˆ2)\n', aCBt)
fprintf('|aCBt| = %g (m/sˆ2)\n', norm(aCBt))
fprintf('\n')
omega5z = sym('omega5z','real');
omega5 = [ 0 0 omega5z ];
vD42 = sym('vD42','real');
vD4D2 = [ vD42*cos(phi4) vD42*sin(phi4) 0 ];
vE = [0 0 0 ];
vD5 = vE + cross(omega5,rD-rE);
vD4 = vD5;
eqvD = vD4 - ( vB2 + vD4D2);
eqvDx = eqvD(1); eqvDy = eqvD(2);
solvD = solve(eqvDx,eqvDy);
omega5zs = eval(solvD.omega5z);
vD42s = eval(solvD.vD42);
Omega5 = [0 0 omega5zs];
VD42 = vD42s*[cos(phi4) sin(phi4) 0];
VD5 = vE + cross(Omega5,rD-rE);
fprintf('vD4 = vD2 + vD4D2 => \n')
qvDx=vpa(eqvDx,6);
fprintf('x-axis:\n')
fprintf('%s=0\n',char(qvDx))
qvDy=vpa(eqvDy,6);
fprintf('y-axis:\n')
fprintf('%s=0\n',char(qvDy))
fprintf('=>\n')
fprintf('vD42 = %g (m/s)\n\n', vD42s)
fprintf('omega5 = [%g, %g, %g](rad/s)\n', Omega5)
fprintf('vD4D2 = [ %g, %g, %d ] (m/s)\n', VD42)
fprintf('\n')
fprintf('vD5 = [ %g, %g, %d ] (m/s)\n', VD5)
aE=[0 0 0];
ep5z = sym('ep5z','real');
ep5 = [ 0 0 ep5z ];
% parallel to the sliding direction
Ep5 = [0 0 ep5zs];
dot(Omega5,Omega5)*(rD-rE);
% print the equations for calculating
% alpha3z and aB21
fprintf...
fprintf('\n')
fprintf...
fprintf('x-axis:\n')
fprintf('y-axis:\n')
fprintf('=>\n')
fprintf('\n')
fprintf...
('ep5 = [ %g, %g, %g ] (rad/sˆ2)\n', Ep5)
fprintf('\n')
omega45 = omega4 - Omega5;
ep45 = ep4 - Ep5;
fprintf...
('omega45 = [ %g, %g, %g ] (rad/s)\n', omega45)
fprintf...
('ep45 = [ %g, %g, %g ] (rad/sˆ2)\n', ep45)
% end of program

### Categories

Find more on MuPAD in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by