Hi everyone,
I'm trying to make a position analysis of a mechanism via Newton-Raphson algorithm. When end of analysis ı can see graphics but ı see "Error nr" at command window. which part is wrong in my code?resim.JPG
%CLEAR: Variables and command window in MATLAB
clc,clear
% INPUT: PHYSICAL PARAMETERS of MECHANISM (meter)
a1=81;a2=20;a3=143;a4=67;a5=95;a6=25;a7=16;b=93;th2=18*pi/180
% INPUT: Maximum Iteration Number Nmax
Nmax=100;
% INPUT: INITIAL GUESS VALUES for th3, th4, th6, s respectively
x=[50*pi/180,40*pi/180,253*pi/180,5];
% INPUT: ERROR TOLERANCE
xe=0.000000001*abs(x);
% INPUT: SYSTEM INPUTS (th2,w2,al2)
dth=5*pi/180;
th2=-30*pi/180:dth:180*pi/180;
w2=-18*ones(1,length(th2));
al2=0*ones(1,length(th2));
%----------------------------------------------
xe=transpose(abs(xe));
kerr=1; %If kerr=1, results are not converged
for k=1:1:length(th2);
for n=1:Nmax
%----------------------------------------------
%Assign initial guess to unknowns
th3(k)=x(1);
th4(k)=x(2);
th6(k)=x(3);
s(k)=x(4);
% INPUT: JACOBIAN Matrix
J=zeros(4,4);
J(1,1)=-a3*sin(th3(k));J(1,2)=a4*sin(th4(k));
J(2,1)=a3*cos(th3(k));J(2,2)=-a4*cos(th4(k));
J(3,1)=-a5*sin(th3(k))-a6*sin(th3(k)+pi/2);J(3,3)=-a7*sin(th6(k)-pi/2)+s(k)*sin(th6(k));J(3,4)=-cos(th6(k));
J(4,1)=a5*cos(th3(k))+a6*cos(th3(k)+pi/2);J(4,3)=a7*cos(th6(k)-pi/2)-s(k)*cos(th6(k));J(4,4)=-sin(th6(k));
% INPUT: Function f
f=zeros(4,1);
f(1,1)=-(a2*cos(th2(k))+a3*cos(th3(k))-a4*cos(th4(k))-a1);
f(2,1)=-(a2*sin(th2(k))+a3*sin(th3(k))-a4*sin(th4(k)));
f(3,1)=-(a2*cos(th2(k))+a5*cos(th3(k))+a6*cos(th3(k)+pi/2)+a7*cos(th6(k)-pi/2)-s(k)*cos(th6(k)));
f(4,1)=-(a2*sin(th2(k))+a5*sin(th3(k))+a6*sin(th3(k)+pi/2)+a7*sin(th6(k)-pi/2)-s(k)*sin(th6(k))-b);
%----------------------------------------------
eps=inv(J)*f;x=x+transpose(eps);
if abs(eps)<xe
kerr=0;break
end
end
if kerr==1
'Error nr'
end
th3(k)=x(1);
th4(k)=x(2);
th6(k)=x(3);
s(k)=x(4);
%---velocity---------------------------
fv(1,1)=-w2(k)*a2.*sin(th2(k));
fv(2,1)=w2(k)*a2.*cos(th2(k));
fv(3,1)=-w2(k)*a2.*sin(th2(k));
fv(4,1)=-w2(k)*a2.*cos(th2(k));
vel=inv(J)*fv;
w3(k)=vel(1);
w4(k)=vel(2);
w6(k)=vel(3);
Vd(k)=vel(4);
%---acceleration---------------------------
fa(1,1)=-al2(k)*a2*sin(th2(k))-w2(k)^2*a2*cos(th2(k))-w3(k)^2*a3*cos(th3(k))+w4(k)^2*a4*cos(th4(k));
fa(2,1)=al2(k)*a2*cos(th2(k))-w2(k)^2*a2*sin(th2(k))-w3(k)^2*a3*sin(th3(k))+w4(k)^2*a4*sin(th4(k));
fa(3,1)=-al2(k)*a2*sin(th2(k))-w2(k)^2*a2*cos(th2(k))-a5*w3(k)^2*cos(th3(k))-a6*w3(k)^2*cos(th2(k)+pi/2)+2*Vd(k)*w6(k)*sin(th6(k))+s(k)*w6(k)^2*cos(th6(k));
fa(4,1)=al2(k)*a2*cos(th2(k))-a2*w2(k)^2*sin(th2(k))-a5*w3(k)^2*sin(th3(k))-a6*w3(k)^2*sin(th2(k)+pi/2)+a7*w6(k)^2*sin(th6(k)-pi/2)-2*Vd(k)*w6(k)*cos(th6(k))-s(k)*w6(k)^2*sin(th6(k));
acc=inv(J)*fa;
al3(k)=acc(1);
al4(k)=acc(2);
al6(k)=acc(3);
ad(k)=acc(4);
end
% Angle: radian --> degree
th2d=th2*180/pi;
th3d=th3*180/pi;
th4d=th4*180/pi;
th6d=th6*180/pi;
%--------Plots---------------
figure(1),
subplot(4,3,1),plot(th2d,th3d,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\theta_3(^o)'),xlim([0 180])
subplot(4,3,2),plot(th2d,w3,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\omega_3(r/s)'),xlim([0 180])
subplot(4,3,3),plot(th2d,al3,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\alpha_3(r/s^2)'),xlim([0 180])
subplot(4,3,4),plot(th2d,th4d,'r','linewidth',2),xlabel('\theta_2(^o)'),ylabel('\theta_4 (^o)'),xlim([0 180])
subplot(4,3,5),plot(th2d,w4,'r','linewidth',2),xlabel('\theta_2(^o)'),ylabel('\omega_4 (r/s)'),xlim([0 180])
subplot(4,3,6),plot(th2d,al4,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\alpha_4(r/s^2)'),xlim([0 180])
subplot(4,3,7),plot(th2d,th6,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\theta_6(^o)'),xlim([0 180])
subplot(4,3,8),plot(th2d,w6,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\omega_6(r/s)'),xlim([0 180])
subplot(4,3,9),plot(th2d,al6,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('\alpha_6(r/s^2)'),xlim([0 180])
subplot(4,3,10),plot(th2d,s,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('s(cm)'),xlim([0 180])
subplot(4,3,11),plot(th2d,Vd,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('V_d(cm/s)'),xlim([0 180])
subplot(4,3,12),plot(th2d,ad,'r','linewidth',2),xlabel('\theta_2 (^o)'),ylabel('a_d(cm/s^2)'),xlim([0 180])

답변 (1개)

KALYAN ACHARJYA
KALYAN ACHARJYA 2019년 8월 17일
편집: KALYAN ACHARJYA 2019년 8월 17일

0 개 추천

Its not a Matlab Error.Just display the string as 'Error nr', when ker = 1
Diaplay from this part of the code
if kerr==1
'Error nr'
end
if you delete that part of the code, still, there will be no issue on results. Here it just used as flag to inform when ker=1;

카테고리

도움말 센터File Exchange에서 Numerical Integration and Differential Equations에 대해 자세히 알아보기

제품

릴리스

R2016b

태그

질문:

2019년 8월 17일

댓글:

2019년 8월 17일

Community Treasure Hunt

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

Start Hunting!

Translated by