Index in position 3 exceeds array bounds. (Must not exceed 1).

I'm trying to run this code but every time I get the error: Index in position 3 exceeds array bounds (must not exceed 1).

Error in Design_3303 (line 28)

if AOA>theta1(i1,i2,i3)

How can I fix this?

Code:

clear;
clc;
close all
chord=1;
AOA=8*ones(1,101);
M_infinity=2.5*ones(1,101);
gamma=1.4*ones(1,101);

x_u=0:0.01:1;
x_l=0:0.01:1;
t_u=0:0.001:0.1;
t_l=0.1-t_u;

for  i1=1:length(x_u)
    for i2=1:length(x_l)
        for i3=1:length(t_u)
            M1=M_infinity;                        
            tmax=chord*(t_u+t_l);
            theta1=abs(atand(t_u./x_u));
            theta2=abs(atand(t_u./(1-x_u)));
            theta3=abs(atand(t_l./x_l));
            theta4=abs(atand(t_l./(1-x_l)));
            Hyp1=sqrt(t_u.^2+x_u.^2);
            Hyp2=sqrt(t_u.^2+(1-x_u).^2);
            Hyp3=sqrt(t_l.^2+x_l.^2);
            Hyp4=sqrt(t_l.^2+(1-x_l).^2);
            
            if AOA>theta1(i1,i2,i3)
                %expansion waves==>expansion wave relations
                %Prandtl-Meyer
                %Upper surface 1 relations:
                %M1=M_infinity;
                %Total Pressure to freestream pressure
                p0_pin=(1+((gamma-1)/2).*M1.^2).^(gamma/(gamma-1));
                theta=AOA-theta1;
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M1.^2-1)));
                vf2=atand(sqrt(M1.^2-1));
                v_infinity=vf*vf1-vf2;
                v_u1=theta+v_infinity;
                M_u1=arrayfun(@pm, M1, theta, gamma);
                p0u1_pu1=(1+((gamma-1)/2).*M_u1.^2).^(gamma/(gamma-1));
                pu1_pin=(1./p0u1_pu1).*p0_pin;
                %Upper surface 2 relations
                theta_umax=theta1+theta2;
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M1.^2-1)));
                vf2=atand(sqrt(M1.^2-1));
                v_infinity=vf*vf1-vf2;
                v_u1=theta+v_infinity;
                v_u2=v_u1+theta_umax;   
                M_u2=arrayfun(@pm, M_u1, theta_umax, gamma);%%%%%%%%%
                p0u2_pu2=(1+((gamma-1)./2).*M_u2.^2).^(gamma/(gamma-1));
                pu2_pin=(1./p0u2_pu2).*(p0_pin);
                %Lower surface 1 relations
                theta_L1=AOA+theta3;
                beta_L1=arrayfun(@oswbeta, M1, theta_L1, gamma);
                Mn_L1=M1.*sind(beta_L1);
                Mn_L2=sqrt(((gamma-1).*Mn_L1.^2+2)./(2.*gamma.*Mn_L1.^2-(gamma-1)));
                pL1_pin=1+((2*gamma)./(gamma+1)).*(Mn_L1.^2-1);
                M_L1=Mn_L2./(sind(beta_L1-theta_L1));
                p0L1_pL1=(1+((gamma-1)./2).*M_L1.^2).^(gamma/(gamma-1));
                %Lower surface 2 relations
                theta_umax2=theta3+theta4; 
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1).*(M_L1.^2-1)));
                vf2=atand(sqrt(M_L1.^2-1));
                v_L1=vf*vf1-vf2;
                v_L2=v_L1+theta_umax2;   
                M_L2=arrayfun(@pm, M_L1, theta_umax2, gamma);
                p0L2_pL2=(1+((gamma-1)/2).*M_L2.^2).^(gamma/(gamma-1));
                pL2_pin=(1./p0L2_pL2).*p0L1_pL1.*pL1_pin;
            elseif AOA==theta1(i1,i2,i3)
                %Upper surface 1 Relations:
                M1=M_infinity;
                p0_pin=(1+((gamma-1)/2)*M1^2)^(gamma/(gamma-1));
                theta=AOA-theta1;
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M1^2-1)));
                vf2=atand(sqrt(M1^2-1));
                v_infinity=vf*vf1-vf2;
                v_u1=theta+v_infinity;
                M_u1=pm(M1, theta, gamma);
                p0u1_pu1=(1+((gamma-1)/2)*M_u1^2)^(gamma/(gamma-1));
                pu1_pin=(1/p0u1_pu1)*p0_pin;
                %Upper Surface 2
                theta_umax=theta1+theta2;
                v_u2=v_u1+theta_umax;
                M_u2=pm(M_u1, theta_umax, gamma);
                p0u2_pu2=(1+((gamma-1)/2)*M_u2^2)^(gamma/(gamma-1));
                pu2_pin=(1/p0u2_pu2)*(p0_pin);
                %Lower surface 1 relations
                theta_L1=AOA+theta3;
                beta_L1=oswbeta(M1, theta_L1, gamma);
                Mn_L1=M1*sind(beta_L1);
                Mn_L2=sqrt(((gamma-1)*Mn_L1^2+2)/(2*gamma*Mn_L1^2-(gamma-1)));
                pL1_pin=1+((2*gamma)/(gamma+1))*(Mn_L1^2-1);
                M_L1=Mn_L2/(sind(beta_L1-theta_L1));
                %Lower surface 2 relations
                theta_umax2=theta3+theta4; 
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M_L1^2-1)));
                vf2=atand(sqrt(M_L1^2-1));
                v_L1=vf*vf1-vf2;
                v_L2=v_L1+theta_umax2;   
                M_L2=pm(M_L1, theta_umax2, gamma);
                p0L2_pL2=(1+((gamma-1)/2)*M_L2^2)^(gamma/(gamma-1));
                pL2_pin=(1/p0L2_pL2)*p0L1_pL1*pL1_pin;
            elseif AOA<theta1(i1,i2,i3) %==>oblique shockwave
                %Upper Surface 1 Relations
                M1=M_infinity;
                theta=AOA;
                beta=oswbeta(M1, theta, gamma);
                M_n1=M1*sind(beta);
                M_n2=sqrt(((gamma-1)*M_n1^2+2)/(2*gamma*M_n1^2-(gamma-1)));
                pu1_pin=1+((2*gamma)/(gamma+1))*(M_n1^2-1);
                M_u1=M_n2/(sind(beta-theta));
                %Upper Surface 2 Relations
                theta_umax=theta1+theta2;
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M1^2-1)));
                vf2=atand(sqrt(M1^2-1));
                v_infinity=vf*vf1-vf2;
                v_u1=theta+v_infinity;
                v_u2=v_u1+theta_umax;
                M_u2=pm(M_u1, theta_umax, gamma);
                p0u2_pu2=(1+((gamma-1)/2)*M_u2^2)^(gamma/(gamma-1));
                pu2_pin=(1/p0u2_pu2)*(p0_pin);
                
                %Lower surface 1 relations
                theta_L1=AOA+theta3;
                beta_L1=oswbeta(M1, theta_L1, gamma);
                Mn_L1=M1*sind(beta_L1);
                Mn_L2=sqrt(((gamma-1)*Mn_L1^2+2)/(2*gamma*Mn_L1^2-(gamma-1)));
                pL1_pin=1+((2*gamma)/(gamma+1))*(Mn_L1^2-1);
                M_L1=Mn_L2/(sind(beta_L1-theta_L1));
                
                cp_u1=(pu1_pin-1)/(0.5*gamma*M_infinity^2);
                cp_u2=(pu2_pin-1)/(0.5*gamma*M_infinity^2);
                cp_L1=(pL1_pin-1)/(0.5*gamma*M_infinity^2);
                cp_L2=(pL2_pin-1)/(0.5*gamma*M_infinity^2);
                
                %Normal and axial force coefficients
                c_n=0.5*((cp_L1+cp_L2)-(cp_u1+cp_u2)); %Normal
                c_a=0.5*(cp_u1*tand(theta1)+cp_u2*tand(-theta2)-cp_L1*tand(-theta3)-cp_L2*tand(theta4));%Axial
                
                %Coefficient of lift
                c_l=c_n*cosd(AOA)-c_a*sind(AOA);
                
                %Drag coefficient
                c_d=c_n*sind(AOA)+c_a*cosd(AOA);
                %cm_LE=cp_u1*((0.5-0)/2)*(0.5-0)+cp_u2*((1+0.5)/2)*(1-0.5)+cp_L1*tan(theta3)*(
            end
        end
    end
end

[Walter Roberson]

Dknow, you marked this as Duplicate. However it appears to be your only question posted. What is it a duplicate of?

[Stephen23]

https://web.archive.org/web/20211120185229/https://www.mathworks.com/matlabcentral/answers/1590124-index-in-position-3-exceeds-array-bounds-must-not-exceed-1

Answer 1

[Cris LaPierre]

Your variable theta1 is apparently an mxnx1 array, and your code is telling MATLAB to index the 3rd dimension with a value >1.

theta1 = rand(3,2,1);
% Works
theta1(1,2,1)

ans = 0.7707

% Doesn't work
theta1(1,2,3)

Index in position 3 exceeds array bounds. Index must not exceed 1.

In your code, theta1 is a 1x101 vector, so indexing the third dimension does not make sense.

[Dknow]

well, I need the three for loops to have matlab iterate the equations each time values for x_u, x_l, t_u and t_l change and also, if the conditions in the if elseif statements are met. I had to index those variables so as to get the resulting solutions in an array. How am I supposed to change the code to achieve these?

[Cris LaPierre]

How do the 3 loops relate to theta1? Currently you are expecting theta1 to be a 3D array, but it is not. Should it be?

[Dknow]

theta1 is supposed to be dependent on x_u (loop 1) and t_u (loop 3)

[Dknow]

I'm not sure if I have them indexed correctly based on my x_u, t_u, t_l and x_l arrays above

[Dknow]

theta1=abs(atand(t_u./x_u));

[Dknow]

Is making it a 3D array the only way the code will work? If so, how do I do that? Is there a way I can run the code without making Theta1 3D?

[Cris LaPierre]

I am not familiar with what you are trying to do, but I suspect it can be done without needing a 3D array. You maybe could use a refresher on indexing. Consider going through Ch 5 of MATLAB Onramp.

Answer 2

[Walter Roberson]

You have a lot of confusion about whether you are dealing with arrays or scalars. Your triple loop structure is the kind of coding you would use for working with scalars all over the place, but you incorrectly work with entire vectors instead, and it gets really confusion as to which element of the vector you should be using instead of the entire vector.

You also hard-code the number 101 for various sizes, which just happens to be the size of several of your vectors. But which 101 corresponds to which vector ? If t_u were to change to be 1001 elements instead of 101 then how would that affect your code? (Oh look, t_u is 1001 long already...)

In the below, I have adjusted your code so that theta1 is consistently indexed, as you wanted to do that in your code. It is not clear why you wanted to index theta1, as you do not use it after the loop, and you only ever appear to intend to use one element at a time.

I did not index theta2, theta3, or theta4, as you never indexed them yourselves.

In the code you will find a few lines marked % ? !! which are lines that you need to edit.

AOA is 101 elements? Is that intended to match the length of x_u, the length of x_l, or the length of t_u ? I store the vector lengths in variables named num_* and on the first three of the % ? !! lines you need to fill in the appropriate name length.

Just after the for i3 line, you will see another three lines marked % ? !! in which AOA, gamma, and M_infinity are assigned to by selecting one particular value as the "current" value out of the vector of possible values. But I cannot tell from your code whether you want to select the AOA associated with the "current" x_u, or the AOA associated with the "current" x_l (both of which have 101) elements, or even the AOA associated with the "current" t_u (which has 1001 elements but you might have overlooked that when you hard-coded 101)

Notice that I create vectors of different names holding all of the allowed values for AOA, M_infinity, and gamma, and in the code pull out the "current" values of each. This avoids confusing using the names as scalars and using the names as the entire vectors.

We know that you were confusing them because you test AOA>theta1(i1,i2,i3) but theta(i1,i2,i3) is a scalar and if you really intended to be comparing the entire vector of values for AOA there, you would be comparing a vector to a scalar in an if statement. Which is permitted, but means the same thing as if all(AOA>theta1(i1,i2,i3)) -- only true if all of the AOA values are greater than theta1(i1,i2,i3) . Which is confusing, and is also a waste of time (you could just pre-compute the minimum of AOA and compare that; if the minimum is greater than all the other values will be greater as well.)

Now... could we talk about why you are creating vectors of values for AOA, M_infinity, and gamma, but the vectors are all the same value? Why bother to create them as vectors, why not just scalars? Or did I miss something in your code where you are using algebraic matrix multiplication against another vector you created as a column vector, effectively calculating a dot product? But if you were intending to do that, then since the elements are all the same, the dot product would be the scalar value times the sum of the other vector, and no vector of identical constants would be needed...

x_u=0:0.01:1;
x_l=0:0.01:1;
t_u=0:0.001:0.1;

num_xu = length(x_u);
num_xl = length(x_l);
num_tu = length(t_u);

t_l=0.1-t_u;

chord=1;
AOA_vals = 8*ones(1,num_??); % ? !!
M_infinity_vals = 2.5*ones(1,num_??); % ? !!
gamma_vals = 1.4*ones(1,num_??); % ? !!

theta1 = zeros(num_xu, num_xl, num_tu);  %

for  i1=1:num_xu
    for i2=1:num_xl
        for i3=1:num_tu
            AOA = AOA_vals(i?);      % ? !!
            gamma = gamma_vals(i?);  % ? !!
            M_infinity = M_infinity_vals(i?); % ? !!
            M1=M_infinity(i?); % ? !!                
            tmax=chord*(t_u(i3)+t_l(i2));  %
            theta1(i1,i2,i3) = abs(atand(t_u(i3)./x_u(i1))); %
            theta2=abs(atand(t_u(i3)./(1-x_u(i1)))); %
            theta3=abs(atand(t_l(i3)./x_l(i2))); %
            theta4=abs(atand(t_l(i3)./(1-x_l(i2)))); %
            Hyp1=sqrt(t_u(i3).^2+x_u(i1).^2); %
            Hyp2=sqrt(t_u(i3).^2+(1-x_u(i1)).^2); %
            Hyp3=sqrt(t_l(i3).^2+x_l(i2).^2); %
            Hyp4=sqrt(t_l(i3).^2+(1-x_l(i2)).^2); %
            
            if AOA > theta1(i1,i2,i3)              %
                %expansion waves==>expansion wave relations
                %Prandtl-Meyer
                %Upper surface 1 relations:
                %M1=M_infinity;
                %Total Pressure to freestream pressure
                p0_pin=(1+((gamma-1)/2).*M1.^2).^(gamma/(gamma-1));
                theta = AOA - theta1(i1,i2,i3); %
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M1.^2-1)));
                vf2=atand(sqrt(M1.^2-1));
                v_infinity=vf*vf1-vf2;
                v_u1=theta+v_infinity;
                M_u1=arrayfun(@pm, M1, theta, gamma);
                p0u1_pu1=(1+((gamma-1)/2).*M_u1.^2).^(gamma/(gamma-1));
                pu1_pin=(1./p0u1_pu1).*p0_pin;
                %Upper surface 2 relations
                theta_umax = theta1(i1,i2,i3) + theta2; %
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M1.^2-1)));
                vf2=atand(sqrt(M1.^2-1));
                v_infinity=vf*vf1-vf2;
                v_u1=theta+v_infinity;
                v_u2=v_u1+theta_umax;   
                M_u2=arrayfun(@pm, M_u1, theta_umax, gamma);%%%%%%%%%
                p0u2_pu2=(1+((gamma-1)./2).*M_u2.^2).^(gamma/(gamma-1));
                pu2_pin=(1./p0u2_pu2).*(p0_pin);
                %Lower surface 1 relations
                theta_L1 = AOA + theta3;
                beta_L1=arrayfun(@oswbeta, M1, theta_L1, gamma);
                Mn_L1=M1.*sind(beta_L1);
                Mn_L2=sqrt(((gamma-1).*Mn_L1.^2+2)./(2.*gamma.*Mn_L1.^2-(gamma-1)));
                pL1_pin=1+((2*gamma)./(gamma+1)).*(Mn_L1.^2-1);
                M_L1=Mn_L2./(sind(beta_L1-theta_L1));
                p0L1_pL1=(1+((gamma-1)./2).*M_L1.^2).^(gamma/(gamma-1));
                %Lower surface 2 relations
                theta_umax2=theta3+theta4; 
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1).*(M_L1.^2-1)));
                vf2=atand(sqrt(M_L1.^2-1));
                v_L1=vf*vf1-vf2;
                v_L2=v_L1+theta_umax2;   
                M_L2=arrayfun(@pm, M_L1, theta_umax2, gamma);
                p0L2_pL2=(1+((gamma-1)/2).*M_L2.^2).^(gamma/(gamma-1));
                pL2_pin=(1./p0L2_pL2).*p0L1_pL1.*pL1_pin;
            elseif AOA == theta1(i1,i2,i3) %
                %Upper surface 1 Relations:
                M1 = M_infinity;
                p0_pin=(1+((gamma-1)/2)*M1^2)^(gamma/(gamma-1));
                theta = AOA - theta1(i1,i2,i3); %
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M1^2-1)));
                vf2=atand(sqrt(M1^2-1));
                v_infinity=vf*vf1-vf2;
                v_u1=theta+v_infinity;
                M_u1=pm(M1, theta, gamma);
                p0u1_pu1=(1+((gamma-1)/2)*M_u1^2)^(gamma/(gamma-1));
                pu1_pin=(1/p0u1_pu1)*p0_pin;
                %Upper Surface 2
                theta_umax = theta1(i1,i2,i3) + theta2;
                v_u2=v_u1+theta_umax;
                M_u2=pm(M_u1, theta_umax, gamma);
                p0u2_pu2=(1+((gamma-1)/2)*M_u2^2)^(gamma/(gamma-1));
                pu2_pin=(1/p0u2_pu2)*(p0_pin);
                %Lower surface 1 relations
                theta_L1 = AOA + theta3;
                beta_L1=oswbeta(M1, theta_L1, gamma);
                Mn_L1=M1*sind(beta_L1);
                Mn_L2=sqrt(((gamma-1)*Mn_L1^2+2)/(2*gamma*Mn_L1^2-(gamma-1)));
                pL1_pin=1+((2*gamma)/(gamma+1))*(Mn_L1^2-1);
                M_L1=Mn_L2/(sind(beta_L1-theta_L1));
                %Lower surface 2 relations
                theta_umax2=theta3+theta4; 
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M_L1^2-1)));
                vf2=atand(sqrt(M_L1^2-1));
                v_L1=vf*vf1-vf2;
                v_L2=v_L1+theta_umax2;   
                M_L2=pm(M_L1, theta_umax2, gamma);
                p0L2_pL2=(1+((gamma-1)/2)*M_L2^2)^(gamma/(gamma-1));
                pL2_pin=(1/p0L2_pL2)*p0L1_pL1*pL1_pin;
            elseif AOA < theta1(i1,i2,i3) %==>oblique shockwave  %
                %Upper Surface 1 Relations
                M1 = M_infinity;
                theta = AOA;
                beta=oswbeta(M1, theta, gamma);
                M_n1=M1*sind(beta);
                M_n2=sqrt(((gamma-1)*M_n1^2+2)/(2*gamma*M_n1^2-(gamma-1)));
                pu1_pin=1+((2*gamma)/(gamma+1))*(M_n1^2-1);
                M_u1=M_n2/(sind(beta-theta));
                %Upper Surface 2 Relations
                theta_umax = theta1(i1,i2,i3) + theta2; %
                vf=sqrt((gamma+1)/(gamma-1));
                vf1=atand(sqrt((gamma-1)/(gamma+1)*(M1^2-1)));
                vf2=atand(sqrt(M1^2-1));
                v_infinity=vf*vf1-vf2;
                v_u1=theta+v_infinity;
                v_u2=v_u1+theta_umax;
                M_u2=pm(M_u1, theta_umax, gamma);
                p0u2_pu2=(1+((gamma-1)/2)*M_u2^2)^(gamma/(gamma-1));
                pu2_pin=(1/p0u2_pu2)*(p0_pin);
                
                %Lower surface 1 relations
                theta_L1 = AOA + theta3;
                beta_L1=oswbeta(M1, theta_L1, gamma);
                Mn_L1=M1*sind(beta_L1);
                Mn_L2=sqrt(((gamma-1)*Mn_L1^2+2)/(2*gamma*Mn_L1^2-(gamma-1)));
                pL1_pin=1+((2*gamma)/(gamma+1))*(Mn_L1^2-1);
                M_L1=Mn_L2/(sind(beta_L1-theta_L1));
                
                cp_u1=(pu1_pin-1)/(0.5*gamma*M_infinity^2);
                cp_u2=(pu2_pin-1)/(0.5*gamma*M_infinity^2);
                cp_L1=(pL1_pin-1)/(0.5*gamma*M_infinity^2);
                cp_L2=(pL2_pin-1)/(0.5*gamma*M_infinity^2);
                
                %Normal and axial force coefficients
                c_n=0.5*((cp_L1+cp_L2)-(cp_u1+cp_u2)); %Normal
                c_a = 0.5*(cp_u1*tand(theta1(i1,i2,i3))+cp_u2*tand(-theta2)-cp_L1*tand(-theta3)-cp_L2*tand(theta4));%Axial
                
                %Coefficient of lift
                c_l = c_n*cosd(AOA) - c_a*sind(AOA);
                
                %Drag coefficient
                c_d = c_n*sind(AOA) + c_a*cosd(AOA);
                %cm_LE=cp_u1*((0.5-0)/2)*(0.5-0)+cp_u2*((1+0.5)/2)*(1-0.5)+cp_L1*tan(theta3)*(
            end
        end
    end
end

[Walter Roberson]

Also, at various places you use arrayfun() to calculate between values that look to me as if they should be scalars.

If you are going to vectorize, then vectorize -- but in such a case you should have coded with at least one : as an index, such as theta1(i1, i2, :) ... and if you do that then you should probably be using logical indexing to select branches instead of if statements.