Index in position 1 exceeds array bounds when it doesn't?

10 views (last 30 days)
Hey guys,
This is my first time posting on Matlab forums, hopefully I don't break any rules, but please don't be harsh if I do.
I am a mechanical engineering student in Hungary, and I am struggling to get a code working. I am working on a 4 bar mechanism, and I am currently trying to get my kinetostatistics values for the joints. On bar1 I have a constant angular velocity torque providing motor, so I am trying to account for that with 'Mhajto' in my code, which stands for the torque of said motor.
I'll enclose my code, so you can have a look at it. If I remove the part "+Mhajto" the code works fine. As soon as I introduce it to my code, it falls apart. Please help, I have 2 weeks to fix this :D
%Dénes Péter ÓE BGK Mechanizmus elmélet 2. házifeladat
%-----------------------PISZKOZAT--------------------
%cosinus tétel BD-re
%BD^2=AB^2+AD^2-2*AB*AD*cos(fiABD)
%BD=AB^2+AD^2-2*AB*AD*cos(fiABD);
%rBD=[BD*cos(fiABD); BD*sin(fiABD); 0]
%-----------------------PISZKOZAT VÉGE-----------------
clc
clear
%Ismert adatok
w1=60;
e1=0;
m1=3;
m2=6;
m3=4;
F0=[3; 0; 0];
%fiABC= deg2rad(180-143.42);
fiBAD= deg2rad(10);%(input("Írj be egy számot, pl 10.0: "));
%fiBC= deg2rad(46.58);
%fiADC=deg2rad(180-59.22);
%fiABD= deg2rad(166.70);
L1=0.25*1000;
L2=0.65*1000;
L3=0.6*1000;
DB=L3;
AD=1*1000;
aA=[ 0; 0; 0];
vAx=0;
vAy=0;
vDx=0;
vDy=0;
vA=[vAx;vAy;0];
aD=[0 ; 0; 0];
vD=[vDx ; vDy; 0];
AB=L1;
BA=AB;
BC=L2;
CB=BC;
CD=L3;
DC=CD;
om1=[0; 0; w1];
%KÖREGYENLETEK
"A csuklók pozíciói"
A= [0;0;0]
D= [1000;0;0];
%Innentől kéne a for ciklus
H01= [cos(fiBAD) -sin(fiBAD) 0; sin(fiBAD) cos(fiBAD) 0; 0 0 1];
H12=[1 0 L1; 0 1 0; 0 0 1];
%H23=[cos(fiABC) -sin(fiABC) 0; sin(fiABC) cos(fiABC) 0; 0 0 1];
re2=[0; 0; 1];
rB=H01*H12*re2; %Ez is jó
B=[rB(1,1); rB(2,1); 0]
syms x y
eqn1=(x-B(1,1)).^2+(y-B(2,1)).^2 == L2^2;
eqn2=(x-D(1,1)).^2+(y-D(2,1)).^2 == L3^2;
megoldas = solve([eqn1, eqn2], [x,y]);
CX = double(megoldas.x);
CY = double(megoldas.y);
C=[CX(2,1); CY(2,1); 0]
fiABC=rad2deg(acos((AB^2+BC^2-vecnorm(C)^2)/(2*AB*BC)));
fiABC_kieg= deg2rad(180 - fiABC);
D
%BDhossz=sqrt((B(1,1)-D(1,1))^2+(B(2,1)-D(2,1))^2);
%CIKLUSOZÁS
%for loop ágyazással
"Helyvektorok"
rAB=[B(1,1)-A(1,1); B(2,1)-A(2,1);0]
rBC=[C(1,1)-B(1,1); C(2,1)-B(2,1);0]
rDC=[C(1,1)-D(1,1); C(2,1)-D(2,1);0]
%rAB=[AB*cos(fiBAD); AB*sin(fiBAD); 0]
%rDC=[DC*cos(fiADC); DC*sin(fiADC); 0]
%rBC=[BC*cos(fiBC); BC*sin(fiBC); 0]
%idáig jó
%H01b= [cos(fiBAD) -sin(fiBAD) 0; sin(fiBAD) cos(fiBAD) 0; 0 0 1];
%H12b=[1 0 L1; 0 1 0; 0 0 1];
%H23b=[cos(fiABC_kieg) -sin(fiABC_kieg) 0; sin(fiABC_kieg) cos(fiABC_kieg) 0; 0 0 1];
%%re2b=[L2; 0; 1];
%rC=H01b*H12b*H23b*re2b; %Ez is jó
%%C=[rC(1,1); rC(2,1); 0];
%"Redeklaráció, javítandó!"
"Sebességek"
format bank
om1_x_rAB=cross(om1,rAB);
vB=double([vAx+om1_x_rAB(1,1); vAy+om1_x_rAB(2,1);0])
%Egyenletrendszer felírása a C pontra
syms vCx vCy om2 om3 w2 w3 om2_x_rBC om3_x_rDC
%javítás: syms változók között ismert értéket nem adunk meg
om2= [0; 0; w2];
om3= [0; 0; w3];
om2_x_rBC=cross(om2,rBC);
om3_x_rDC=cross(om3,rDC);
vCx= vB(1,1)+om2_x_rBC(1,1)==vD(1,1)+om3_x_rDC(1,1);
vCy= vB(2,1)+om2_x_rBC(2,1)==vD(2,1)+om3_x_rDC(2,1);
sol = solve([vCx,vCy],[w2, w3]);
"w2:";
w2_mo = double(sol.w2)
"w3:";
w3_mo = double(sol.w3)
om2= [0; 0; sol.w2];
om3= [0; 0; sol.w3];
om2_x_rBC=cross(om2,rBC);
om3_x_rDC=cross(om3,rDC);
vCx= vB(1,1)+om2_x_rBC(1,1);
vCy= vB(2,1)+om2_x_rBC(2,1);
vC=double([vCx; vCy; 0])
"Gyorsulás"
eps1=[0; 0; e1];
eps1_x_rAB=cross(eps1,rAB);
aBx=aA(1,1)+eps1_x_rAB(1,1)-w1^2*rAB(1,1);
aBy=aA(2,1)+eps1_x_rAB(2,1)-w1^2*rAB(2,1);
aB=double([aBx; aBy; 0])
syms e2 e3 eps2_x_rBC eps3_x_rDC eps2 eps3 aCx aCy
eps2= [0; 0; e2];
eps3= [0; 0; e3];
eps2_x_rAB= cross(eps2,rAB);
eps3_x_rDC= cross(eps3,rDC);
aCx=aD(1,1)+eps3_x_rDC(1,1)-w3_mo^2*rDC(1,1)==aBx+eps2_x_rAB(1,1)-w2_mo*rBC(1,1);
aCy=aD(2,1)+eps3_x_rDC(2,1)-w3_mo^2*rDC(2,1)==aBy+eps2_x_rAB(2,1)-w2_mo*rBC(2,1);
sol = solve([aCx,aCy],[e2,e3]);
"e2:";
e2_mo = double(sol.e2)
"e3:";
e3_mo = double(sol.e3)
eps2= [0; 0; sol.e2];
eps3= [0; 0; sol.e3];
eps2_x_rAB= cross(eps2, rAB);
eps3_x_rDC= cross(eps3, rDC);
aCx=aD(1,1)+eps3_x_rDC(1,1)-w3_mo^2*rDC(1,1);
aCy=aD(2,1)+eps3_x_rDC(2,1)-w3_mo^2*rDC(2,1);
aC=double([aCx; aCy; 0])
%inerciaerők
mA=0.25*m1;
mP=0.75*m1;
FA=0;
rAP=[2/3*rAB(1,1); 2/3*rAB(2,1); 0];
eps1_x_rAP=cross(eps1, rAP);
aPx= aA(1,1)+eps1_x_rAP(1,1)-w1^2*rAP(1,1);
aPy= aA(2,1)+eps1_x_rAP(2,1)-w1^2*rAP(2,1);
aP=[aPx; aPy; 0]
FP=-mP*aP;
Fi1=FA+FP
Fi1hossz=vecnorm(Fi1);
lambda1=(rAP(1,1)*FP(2,1)-rAP(2,1)*FP(1,1))/(rAB(1,1)*Fi1(2,1)-rAB(2,1)*Fi1(1,1));
double(lambda1);
rAT=rAP
mD=0.25*m3;
mR=0.75*m3;
FD=0 ;
rDR=[2/3*rDC(1,1); 2/3*rDC(2,1); 0];
eps3_x_rDR=cross(eps3, rDR);
aRx= aD(1,1)+eps3_x_rDR(1,1)-w3_mo^2*rDR(1,1);
aRy= aD(2,1)+eps3_x_rDR(2,1)-w3_mo^2*rDR(2,1);
aR=double([aRx; aRy; 0]);
FR=-mR*aR;
Fi3=FD+FR
Fi3hossz=vecnorm(double(Fi3));
lambda3=(rDR(1,1)*FR(2,1)-rDR(2,1)*FR(1,1))/(rDC(1,1)*Fi3(2,1)-rDC(2,1)*Fi3(1,1));
double(lambda3);
rDT=lambda3* rDC
mB=0.25*m2;
mQ=0.75*m2;
FB=[mB*aB(1,1); mB*aB(2,1); 0];
rBQ=[2/3*rBC(1,1); 2/3*rBC(2,1); 0];
eps2_x_rBQ=double(cross(eps2, rBQ));
aQx= aB(1,1)+eps2_x_rBQ(1,1)-w2_mo^2*rBQ(1,1);
aQy= aB(2,1)+eps2_x_rBQ(2,1)-w2_mo^2*rBQ(2,1);
aQ=[aQx; aQy; 0];
FQ=[aQ(1,1)*-mQ; aQ(2,1)*-mQ;0]
Fi2=FB+FQ;
Fi2hossz=vecnorm(Fi2);
lambda2=(rBQ(1,1)*FQ(2,1)-rBQ(2,1)*FQ(1,1))/(rBC(1,1)*Fi2(2,1)-rBC(2,1)*Fi2(1,1));
rBT = rBC * lambda2
"Virtuális teljesítmény"
vT1= double(vA + cross(om1,rAP))
vT2= double(vB + cross(om2,rBT))
vT3= double(vD + cross(om3,rDT))
Fi1_vT=Fi1(1,1)*vT1(1,1)+Fi1(2,1)*vT1(2,1);
Fi2_vT=Fi2(1,1)*vT2(1,1)+Fi2(2,1)*vT2(2,1);
Fi3_vT=Fi3(1,1)*vT3(1,1)+Fi3(2,1)*vT3(2,1);
F0_vC= F0(1,1)*vC(1,1);
Mhajto=double([0; 0; (Fi1_vT+Fi2_vT+Fi3_vT+F0_vC)/w1])
"Kinetostatikai számítások"
rBA=rAB*-1;
rCB=-1*rBC;
rCD=-1*rDC;
rAB
rAT
rBT1=-(rAB-rAT)
rCB
rBT2=rBT
%rDC
%rDT
%rCT3=-(rDC-rDT)
syms FAx FAy FCx FCy
FAx_mo = FAx == ((rBA(1,1)*FAy)+Mhajto)/rBA(2,1)==-Fi1(1,1)-Fi2(1,1)-FCx;
FAy_mo = FAy == ((rBA(2,1)*FAx)+Mhajto)/rBA(1,1)==-Fi1(2,1)-Fi2(2,1)-FCy;
FCx_mo = FCx == (rBT2(1,1)*Fi2(2,1)-rBT2(2,1)*Fi2(1,1)+rBC(1,1)*FAy-F0(1,1)*rBC(2,1))/rBC(2,1)==-Fi1(1,1)-Fi2(1,1)-FAx-F0(1,1);
FCy_mo = FCy == (rBT2(2,1)*Fi2(1,1)-rBT2(1,1)*Fi2(2,1)+rBC(2,1)*FAx+F0(2,1)*rBC(2,1))/rBC(1,1)==-Fi1(2,1)-Fi2(2,1)-FAy-F0(2,1);
megoldasegyenlet = solve([FAx_mo, FAy_mo, FCx_mo, FCy_mo],[FAx FAy FCx FCy]);
'Az egye cuklókban ébredő erők'
FA=double([megoldasegyenlet.FAx; megoldasegyenlet.FAy;0])
FC=double([megoldasegyenlet.FCx; megoldasegyenlet.FCy;0])
FBx=-FA(1,1)-Fi1(1,1);
FBy=-FA(2,1)-Fi1(2,1);
FB=double([FBx;FBy;0])
FDx=-FC(1,1)-Fi3(1,1);
FDy=-FC(2,1)-Fi3(2,1);
FD=double([FDx; FDy; 0])
  5 Comments
John D'Errico
John D'Errico on 10 May 2020
I want to upvote Steve's comment at least a hundred times. And then a hundred more.
Peter Denes
Peter Denes on 10 May 2020
Thank you Steve, although I already wrote this code, and I am not going to implement functions into it this time, I already started working on another project in the same subject, and your ideas came in very handy, so I used your suggestions. I even contemplated building a GUI on it, but unfortunately my upcoming exams are more important.

Sign in to comment.

Answers (2)

Ameer Hamza
Ameer Hamza on 3 May 2020
Edited: Ameer Hamza on 3 May 2020
Explicitly considering the error, you mentioned in your question. It happens because 'Mhajto' is a 3x1 vector. And when you add it in the expression of FAx_mo, then FAx_mo also has dimensions of 3x1. FCx_mo is not dependent on Mhajto on Mhajto, so its dimensions are 1x1. I don't understand the logic of your code so I am not sure what is the correct thing to do here but changing the value of Mhajto to the following line will correct this error
Mhajto=double([(Fi1_vT+Fi2_vT+Fi3_vT+F0_vC)/w1]) % Mhajto is 1x1
However, note that with the new of Mhajto, the solve() function returns an empty vector, i.e., no solution exist
megoldasegyenlet = solve([FAx_mo, FAy_mo, FCx_mo, FCy_mo],[FAx FAy FCx FCy]);
which will result in a different error on the line
FBx=-FA(1,1)-Fi1(1,1);
FBy=-FA(2,1)-Fi1(2,1);
because FA has dimensions of 1x1.
An alternate solution is to keep Mhajto as a 3x1 vector and use their single elements in the following equations
FAx_mo = FAx == ((rBA(1,1)*FAy)+Mhajto(1))/rBA(2,1)==-Fi1(1,1)-Fi2(1,1)-FCx;
FAy_mo = FAy == ((rBA(2,1)*FAx)+Mhajto(2))/rBA(1,1)==-Fi1(2,1)-Fi2(2,1)-FCy;
Which of the above solution is correct depends on the formula you are trying to implement.
  12 Comments
Peter Denes
Peter Denes on 9 May 2020
Also sorry for the handwriting, visuals are not really my strongsuit, especially not when they are just drafts. However, I did not know that I would be posting them online, so if you need, I can rewrite them tomorrow!

Sign in to comment.


Peter Denes
Peter Denes on 10 May 2020
  3 Comments
Peter Denes
Peter Denes on 11 May 2020
I just realized I messed up the equtations I inserted here as a picture yesterday...ouch.
Anyway, a friend of mine was able to help me, so I'll attach the code section in case someone needs it in the future.
syms FAx FAy FCx FCy FBx FBy FDx FDy i j k
FA = [FAx; FAy; 0];
FC = [FCx; FCy; 0];
crossA = (cross(rBA,FA)+Mhajto);
crossC = (cross(rBT2,Fi2)+cross(rBC,FC));
FA_mox = crossA(3,1)/rBA(1,1) == 0;
FA_moy = crossA(3,1)/rBA(2,1) == 0;
FC_mox = crossC(3,1)/rBC(1,1) == 0;
FC_moy = crossC(3,1)/rBC(1,1) == 0;
sumFx = FA(1,1)+Fi1(1,1)+Fi2(1,1)+FC(1,1) == 0;
sumFy = FA(2,1)+Fi1(2,1)+Fi2(2,1)+FC(2,1) == 0;
megoldasegyenlet = solve([FA_mox, FA_moy, FC_mox, FC_mox, sumFx, sumFy],[FAx FAy FCx FCy]);
fprintf('Az egyes csuklókban ébredő erők')
FA=double([megoldasegyenlet.FAx; megoldasegyenlet.FAy;0])
FC=double([megoldasegyenlet.FCx; megoldasegyenlet.FCy;0])
FB = [FBx; FBy; 0];
FD = [FDx; FDy; 0];
eqn2 = FA + Fi1 + FB == 0;
eqn3 = +F0 - FC + Fi3 + FD == 0;
FB_FD_mo = solve([eqn2, eqn3], [FBx FBy FDx FDy]);
FB = double([FB_FD_mo.FBx; FB_FD_mo.FBy; 0])
FD = double([FB_FD_mo.FDx; FB_FD_mo.FDy; 0])
Also, thanks for all the help and feedback, it was nice talking to you!
Ameer Hamza
Ameer Hamza on 11 May 2020
I am glad that the issue is resolved. Sometimes, when implementing mathematical equations, even a trivial mistake can make the code very difficult to debug.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!