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

3 vues (au cours des 30 derniers jours)
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)

Réponses (2)

Mike Hosea
Mike Hosea le 20 Oct 2014
Modifié(e) : Mike Hosea le 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 commentaire
Dave
Dave le 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.

Connectez-vous pour commenter.


Stalin Samuel
Stalin Samuel le 20 Oct 2014

Catégories

En savoir plus sur Stress and Strain dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by