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

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)

Answer 1

[Mike Hosea]

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.

[Dave]

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.

Answer 2

[Stalin Samuel]

http://www.mathworks.com/matlabcentral/answers/102145-why-do-i-get-the-error-subscript-indices-must-either-be-real-positive-integers-or-logicals