How to get an output off a vectorized function? Error : Subscript indices must either be real positive integers or logicals

1 view (last 30 days)
Dave
Dave on 20 Oct 2014
Commented: Dave on 21 Oct 2014
I'm a new user, so bare with me. I have a function known as Mxy and I want it to return a value for a specified X, such as X=2. What code should I use to achieve this. So far, I get an error (shown below) Here is the relevent code :
%*Plane X-Z
x=0:.001:a+b+c+d
Mxz(0<=x & x<=a)=(Tt+Ts).*x(0<=x & x<=a)
Mxz(a<=x & x<=a+b)=(Tt+Ts).*a+Raz.*(x(a<=x & x<=a+b)-a)
Mxz(a+b<=x & x<=a+b+c+d)=(Tt+Ts).*a+Raz.*((a+b)-a)+Pgas.*(x(a+b<=x & x<=a+b+c+d)-(a+b))
%*X-z Moment Plot
figure(4)
plot(x,Mxz)
grid
title('Bending Moment Diagram : X-Z Plane');
xlabel('X (m)')
ylabel('Bending Moment(Nm)')
hold on
x=.2
Mxz(x)
Here is the error code :
Subscript indices must either be real positive integers or logicals. Error in Crankshaft_1_2 (line 232) Mxz(0.2)
Here is my full code, so far (bare with me... I know it's bad, but as I said, I need a program and it's my first attempt ever) :
%*Inputs given by not yet implemented programs Pengine=7457 %Engine Power Output (W) Theta1dot=377 %Crankshaft angular velocity (rad/s) Rpulley=0.045 %Assumed Crankshaft Pulley Radius (m) P3=4468*10^(3) %Maximum Operating Pressure (Pa) Dbore=0.066 %Cylinder Bore Diameter(m) Rc=0.1015 %Connecting Rod pin-to-pin length (m) Finertia=1000 %Assuming the Crankpin is balanced and case study position
%Assumed Dimensions of not yet determined portions of the crankshaft Lp=0.06 %Pulley Thickness-Assumed (m) Tsr=0.01 %Seal Retainer Thickness (m) Ls=Tsr %Seal Retainer Thickness (Flywheel side) Lch=0.002 %Chamfer Thickness (m) Lfw=0.06 %Flywheel Thickness-Assumed (m) Lmb=0.025 %Journal bearing length (m) Lcr=Lmb %Journal bearing length (m) Dpin=Lcr %Journal/Pin diameter c=0.10 %Cylinder Spacing Axis-to-Axis - Assumed (m) Lcw=c-Lcr %Middle Web based Cylinder Spacing and Journal Length (m) Lw=0.3 %Web Thickness-Assumed
%*Dimentions Short-cut a=0.5*Lp+1.2*Tsr+0.5*Lmb+Lch b=0.5*Lmb+Lw+0.5*Lcr+2*Lch
d=0.5*Lcr+Lw+0.5*Lmb+2*Lch
e=0.5*Lmb+Lch+Ls+0.5*Lfw
Length=a+b+c+d+e
%*Established Main Security Factor SF=2
%*Involved Forces %*Pulley Pull (Ts & Tt) Ts=(0.1*Pengine)/(9*Theta1dot*Rpulley) Tt=10*Ts %*Maximum gas pressure Pgas=P3*((pi/4)*Dbore^2) %*Generator induced Torque Tgen=Rpulley*Tt-Rc*Pgas-Rpulley*Ts
%*Static Analysis : First Case Study - Detailed in Report
Raz=-((d+c)*Pgas+(d+c+b+a)*(Ts+Tt))/(d+c+b)
Rbz=-(Raz+Pgas+Tt+Ts)
Ray=-((d+c)*Finertia)/(d+c+b)
Rby=-Ray
%*Equivalent Reactions
Ra=(Raz^(2)+(Ray))^(1/2)
Rb=(Rbz^(2)+(Rby))^(1/2)
%*FIRST CASE STUDY & SIZING syms x %*Shear Equations & Diagrams (OK-Need Check) %*Plane X-Y x=0:.001:a+b+c+d Fxy(0<x & x<=a)=0 Fxy(a<x & x<=a+b)=Ray Fxy(a+b<x & x<=a+b+c)=Ray+Finertia Fxy(a+b+c<x & x<=a+b+c+d)=Ray+Finertia-Finertia %*X-Y Shear Stress Diagram figure(1) plot(x,Fxy) grid title('Shear Stress Diagram : X-Y Plane'); xlabel('X (m)') ylabel('Shear Stress(N)') hold on %*Plane X-Z x=0:.001:a+b+c+d Fxz(0<x & x<=a)=Tt+Ts Fxz(a<x & x<=a+b)=Tt+Ts+Raz Fxz(a+b<x & x<=a+b+c+d)=Tt+Ts+Raz+Pgas %*X-Z Shear Stress Diagram figure(2) plot(x,Fxz) grid title('Shear Stress Diagram : X-Z Plane'); xlabel('X (m)') ylabel('Shear Stress(N)') hold on
%*Moment Equations & Diagrams (OK-Need Check)
%*Plane X-Y
x=0:.001:a+b+c+d
Mxy(0<=x & x<=a)=0
Mxy(a<=x & x<=a+b)=Ray.*(x(a<=x & x<=a+b)-(a))
Mxy(a+b<=x & x<=a+b+c)=Ray.*((a+b)-(a+b))+Finertia.*(x(a+b<=x & x<=a+b+c)-(a+b))
Mxy(a+b+c<=x & x<=a+b+c+d)=Ray.*((a+b)-(a+b))+Finertia.*((a+b+c)-(a+b))-Finertia.*(x(a+b+c<=x & x<=a+b+c+d)-(a+b+c))
%*X-Y Moment Plot
figure(3)
plot(x,Mxy)
grid
title('Bending Moment Diagram : X-Y Plane');
xlabel('X (m)')
ylabel('Bending Moment(Nm)')
hold on
%*Plane X-Z
x=0:.001:a+b+c+d
Mxz(0<=x & x<=a)=(Tt+Ts).*x(0<=x & x<=a)
Mxz(a<=x & x<=a+b)=(Tt+Ts).*a+Raz.*(x(a<=x & x<=a+b)-a)
Mxz(a+b<=x & x<=a+b+c+d)=(Tt+Ts).*a+Raz.*((a+b)-a)+Pgas.*(x(a+b<=x & x<=a+b+c+d)-(a+b))
%*X-z Moment Plot
figure(4)
plot(x,Mxz)
grid
title('Bending Moment Diagram : X-Z Plane');
xlabel('X (m)')
ylabel('Bending Moment(Nm)')
hold on
%*Torque equations & Diagram (OK-Need Check)
x=0:.001:a+b+c+d
Tx(0<=x & x<=a)=-(Tt-Ts).*Rpulley
Tx(a<=x & x<=a+b+c+d+e)=-(Tt-Ts).*Rpulley+Pgas.*Rc
%*Torque Along X Diagram
figure(5)
plot(x,Tx)
grid
title('Torque along X-Axis')
xlabel('X (m)')
ylabel('Bending Moment (Nm)')
hold on
%*Location of the Critical sections wrt X-Axis %*Shaft & Pins C1=0.5*Lp C2=a-0.5*Lmb-Lch C3=a C4=a+0.5*Lmb C5=a+b-0.5*Lcr C6=a+b C7=a+b+0.5*Lcr C8=a+b+c-Lcr C9=a+b+c C10=a+b+c+Lcr C11=a+b+c+d-0.5*Lmb C12=a+b+c+d C13=a+b+c+d+0.5*Lmb+Lch C14=a+b+c+d+e %*Webs W1=a+0.5*b W2=a+b+0.5*c W3=a+b+c+0.5*d
%*Fatigue design sress computation %*Material Properties (AISI 1040 annealed steel) Su=518.8 %MPa Sy=353.4 %MPa Hb=149 %*Pin Elements %%%%%%%%%%TEMPORARY VALUES RFcs=99.9 %Reliability Factor %*Adjusting Factors CLs=1 %For shaft elements, bending is the only alternating stress if Dpin<10 %mm CGs=1 end if 10<Dpin<50 %mm CGs=0.9 end if 50<Dpin<100 CGs=0.8 end CSs=0.9 %Journals are polished CTs=1 %Lower temperature than the oil degradation temp.
if RFcs==99.9
CRs=0.753
end
if RFcs==99
CRs=0.814
end
if RFcs==95
CRs=0.868
end
if RFcs==90
CRs=0.897
end
if RFcs==50
CRs=1
end
%*Fatigue Stress Computing
Snprimes=0.5*Su
Sns=Snprimes*CLs*CGs*CSs*CTs*CRs
%*Stress Concentration factor %*Notch Factor & Notch Sensitivity %%%%%%%%%TEMPORARY VALUE r=0.0005%m Notch if r==0.5 q=0.6 qs=0.65 end if r==1 q=0.67 qs=0.72 end if r==1.5 q=0.72 qs=0.78 end if r==2 q=0.75 qs=0.8 end %*Stress Concentration factor Kt Dond=1.2 %As 20% steps are common practice, the ratio %D/d is set to 1.2, where d=Djournal. rdRatio=r/Dpin %Taking Larges Kt of each interval rdRatio r Dpin if 0.2<=rdRatio<0.3 Kt=1.3 Kts=1.2 end if 0.1<=rdRatio<0.2 Kt=1.55 Kts=1.3 end if 0.0<=rdRatio<0.1 Kt=2.7 Kts=2.4 end Kf=1+(Kt-1)*q Kfs=1+(Kts-1)*qs %*Determining Journal Critical Section %Assuming square journals (not uncommon) Dpin=Lcr= x=.2 Mxz(x)

Answers (2)

Mike Hosea
Mike Hosea on 20 Oct 2014
Edited: Mike Hosea on 20 Oct 2014
You have not defined Mxz as a function, rather as an array, so Mxz(0.2) asks for the 0.2th element of the array. The syntax you are using is called logical indexing. You have defined a set of discrete (x,y) coordinates, and you could use interp1 to interpolate at an arbitrary x value, or you could write a proper function in a MATLAB file. Actually, the function would look much like what you have there, but you would make x an input argument, and you would not have the line defining x. If the function name us Mxz then you would need a different output variable name for the logical indexing.
However, after thinking about this some, I suspect your simplest course of action is to use interp1. You must keep the x array that you used to define the Mxz array (don't overwrite it with 0.2 or anything else). Then individual values can be computed via interpolation:
y = interp1(x,Mxz,0.2)
will extract the value for x == 0.2 if 0.2 is in the x array, otherwise interpolate for it.
  1 Comment
Dave
Dave on 21 Oct 2014
Thank you. I should have wrote this earlier, but your comment got me unstuck :).
Matlab is interesting & useful. Just started learning it... by necessity. For honor's project, we must design an engine from A to Z in three months. No sleep, only work :P.

Sign in to comment.


Community Treasure Hunt

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

Start Hunting!