Index exceeds the number of array elements (1) error

3 views (last 30 days)
Zafer Duyenli on 18 Dec 2018
Answered: Aoi Midori on 19 Dec 2018
greetings.
I have asked another question about that topic a time ago but I have the same issue again and I can't find the solution no matter how deep I look at into the program.
I am trying to analyze velocities, accelerations and positions of each link in a mechanism by using Newton Raphson method. I have no issue until the analysis of acceleration matrix and I have the issue stated above in line 65. How can I reconstruct that line for solution?
clc, clear;
%constants
a2=29;a3=15;a4=29;a5=27;a6=25;
%max iteration nr
Nmax=100;
%input:initial guess values th3,th5,th6,to
x=[12*pi/180,61*pi/180,323*pi/180,31];
%Error tolerance
xe=0.001*abs(x);
%System inputs
dth=5*pi/180;
th2=290*pi/180:dth:368*pi/180;
w2=-5*ones(1,length(th2));
al2=0*ones(1,length(th2));
xe=transpose(abs(xe));
kerr=1;
for k=1:1:length(th2)
for n=1:Nmax
th3(k)=x(1);th5(k)=x(2);
th6(k)=x(3);to(k)=x(4);
%Jacobian Matrix
J=zeros(4,4);
J(1,1)=-a3*sin(th3(k));J(1,2)=-a4*sin(th5(k)+210*pi/180);
J(2,1)=a3*cos(th3(k));J(2,2)=a4*cos(th5(k)+210*pi/180);
J(3,2)=-a5*sin(th5(k));J(3,3)=-a6*sin(th6(k));J(3,4)=-1;
J(4,2)=a5*cos(th5(k));J(4,3)=a6*cos(th6(k))
%position matrix
f=zeros(4,1);
f(1,1)=-(a2*cos(th2(k))+a3*cos(th3(k))+a4*cos(th5(k)+210*pi/180)-43);
f(2,1)=-(a2*sin(th2(k))+a3*sin(th3(k))+a5*sin(th5(k)+210*pi/180)+33);
f(3,1)=-(a5*cos(th5(k))+a6*cos(th6(k))-to(k));
f(4,1)=-(a5*sin(th5(k))+a6*sin(th6(k))-10);
eps=inv(J)*f;
x=x+transpose(eps);
if abs(eps)<xe
kerr=0;break
end
if kerr==1
'Error nr'
end
th3(k)=x(1);th5(k)=x(2);th6(k)=x(3);to(k)=x(4);
%velocity analysis
fv=zeros(4,1);
fv(1,1)=w2(k)*a2.*sin(th2(k));
fv(2,1)=-w2(k)*a2*cos(th2(k));
vel=inv(J)*fv;
w3(k)=vel(1);w5(k)=vel(2);w6=vel(3);vt(k)=vel(4);
%acceleration analysis
fa=zeros(4,1);
fa(1,1)=al2(k)*a2*sin(th2(k))-w2(k)^2*a2*cos(th2(k))-w3(k)^2*a3*cos(th3(k))-w5(k)^2*a4*cos(th5(k)+210*pi/180);
fa(2,1)=-al2(k)*a2*cos(th2(k))+w2(k)^2*a2*sin(th2(k))+w3(k)^2*a3*sin(th3(k))+w5(k)^2*a4*sin(th5(k)+210*pi/180);
fa(3,1)=w5(k)^2*a5*cos(th5(k))+w6(k)^2*a6*cos(th6(k)); %LINE 65
fa(4,1)=w5(k)^2*a5*sin(th5(k))+w6(k)^2*a6*sin(th6(k));
acc=inv(J)*fa;
al3(k)=acc(1);al5(k)=acc(2);
al6(k)=acc(3);at(k)=acc(4);
end
end
th2d=th2*180/pi;
th3d=th3*180/pi;
th5d=th5*180/pi;
th6d=th6*180/pi;
%--------Plots---------------
figure(1),
subplot(4,3,1),plot(th2d,th4d,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\theta_3 (^o)'),xlim([300 368])
subplot(4,3,2),plot(th2d,w4,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\omega_3 (r/s)'),xlim([300 368])
subplot(4,3,3),plot(th2d,al4,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\alpha_3 (r/s^2)'),xlim([300 368])
subplot(4,3,4),plot(th2d,th5d,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\theta_5 (^o)'),xlim([300 368])
subplot(4,3,5),plot(th2d,w5,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\omega_5 (r/s)'),xlim([300 368])
subplot(4,3,6),plot(th2d,al5,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\alpha_5 (r/s^2)'),xlim([300 368])
subplot(4,3,7),plot(th2d,s,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\theta_6 (cm)'),xlim([300 368])
subplot(4,3,8),plot(th2d,Vs,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\omega_6 (cm/s)'),xlim([300 368])
subplot(4,3,9),plot(th2d,as,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\alpha_6 (cm/s^2)'),xlim([300 368])
subplot(4,3,10),plot(th2d,to,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('t (cm)'),xlim([300 368])
subplot(4,3,11),plot(th2d,Vt,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('V_t (cm/s)'),xlim([300 368])
subplot(4,3,12),plot(th2d,at,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('a_t (cm/s^2)'),xlim([300 368])

Aoi Midori on 19 Dec 2018
I assumed that you wanted to write in LINE 59 as:
w3(k)=vel(1);w5(k)=vel(2);w6(k)=vel(3);vt(k)=vel(4);
where it's not w6=vel(3) but w6(k)=vel(3), otherwise w6 would be just a scalar and you wouldn't be able to access to w6(k) as you wrote in LINE 65.

Categories

Find more on Annotations in Help Center and File Exchange

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by