Wrote this file to compute ABD matrix, running into index/array issue

The code as commented is to calculate the ABD matrix of a material under plane-stress. As far as I can tell there should be no issues running this, yet when run, I receive these errors:

Index in position 1 exceeds array bounds (must not exceed 2).

Error in ESCI325HarnettAssignment63>ReducedStiffnessMatrix (line 67)

Qbar(1,1) = Q(1,1)*(cosd(theta)^4) + 2*(Q(1,2) + 2*Q(3,3))*(sind(theta)^2)*(cosd(theta)^2) + Q(2,2)*(sind(theta)^4);

Error in ESCI325HarnettAssignment63 (line 34)

[Qbar] = ReducedStiffnessMatrix(E1, E2, G12, v12, v21, theta(i));

Any help in a direction to fix this would be appreciated.

% The purpose of this program is to calculate the ABD Matrix for an orthotropic
% material under plane-stress.

close all;

E1   = 5.6*10^6; % psi
E2   = 1.2*10^6; % psi
G12  = .6*10^6; % psi
v12 = .26;
v21 = .0557;
den = (1 - (v12*v21));

Nlayers = 4; % the number of layers
n = Nlayers + 1; % the number of edges

% angle of each layer
theta(1) = 0; 
theta(2) = 0;
theta(3) = 0;
theta(4) = 0;

% edge thickness of each layer in inches
z(1) = .4;
z(2) = .2;
z(3) = .0;
z(4) = -.2;
z(5) = -.4;

A = zeros(3,3);
B = zeros(3,3);
D = zeros(3,3);

for i = 1:Nlayers
  [Qbar] = ReducedStiffnessMatrix(E1, E2, G12, v12, v21, theta(i));
  [A] = ExtensionalMatrix(A, Qbar, z(i), z(i+1));
  [B] = BendTwistCouplingMatrix(B, Qbar, z(i), z(i+1));
  [D] = BendingMatrix(D, Qbar, z(i), z(i+1));
end

[ABD] = [A, B; 

B, D];

function [A] = ExtensionalMatrix(Aprevious, Qbar, ztop, zbottom)
A = zeros(3,3);
A = Aprevious + Qbar*(ztop - zbottom);
end

function [B] = BendTwistCouplingMatrix(Bprevious, Qbar, ztop, zbottom)
B = zeros(3,3);
B = Bprevious + (1/2)*Qbar*(ztop^2 - zbottom^2);
end

function [D] = BendingMatrix(Dprevious, Qbar, ztop, zbottom)
D = zeros(3,3);
D = Dprevious + (1/3)*Qbar*(ztop^3 - zbottom^3);
end

function [Qbar] = ReducedStiffnessMatrix(E1, E2, G12, v12, den, theta)

Q(1,1) = E1/(den);
Q(1,2) = v12*E2/(den);
Q(1,2) = v12*E2/(den);

Q(2,1) = Q(1,2);
Q(2,2) = E2/(den);
Q(2,3) = G12;

Qbar(1,1) = Q(1,1)*(cosd(theta)^4) + 2*(Q(1,2) + 2*Q(3,3))*(sind(theta)^2)*(cosd(theta)^2) + Q(2,2)*(sind(theta)^4);
Qbar(1,2) = (Q(1,1) + Q(2,2) - 4*Q(3,3))*(sind(theta)^2)*(cosd(theta)^2) + Q(1,2)*((sind(theta)^4) + (cosd(theta)^4));
Qbar(1,3) = (Q(1,1) - Q(1,2) - 2*Q(3,3))*(sind(theta))*(cosd(theta)^3) + (Q(1,2) - Q(2,2) + 2*Q(3,3))*(sind(theta)^3)*(cosd(theta));

Qbar(2,1) = Qbar(1,2);
Qbar(2,2) = Q(1,1)*(sind(theta)^4) + 2*(Q(1,2) + 2*Q(3,3))*(sind(theta)^2)*(cosd(theta)^2) + Q(2,2)*(cosd(theta)^4); 
Qbar(2,3) = (Q(1,1) - Q(1,2) - 2*Q(3,3))*(sind(theta)^3)*(cosd(theta)) + (Q(1,2) - Q(2,2) + 2*Q(3,3))*(sind(theta))*(cosd(theta)^3);

Qbar(3,1) = Qbar(1,3);
Qbar(3,2) = Qbar(2,3);
Qbar(3,3) = (Q(1,1) + Q(2,2) - 2*Q(1,2) - 2*Q(3,3))*(sind(theta)^2)*(cosd(theta)^2) + Q(3,3)*((sind(theta)^4) + (cosd(theta)^4));
end

[KALYAN ACHARJYA]

Such errors are very easy to detect. On the error line, some variables that you are trying to access its elements, there may be not having enough lengths.

Use whos to see details

Answer 1

[Dinesh Yadav]

Hi Kevin,

The error here is occurs in ReduceStiffnessMatrix function in the 7th line where you are accessing Q(3,3) but at that point of time the size of Q is . Therefore error occurs. The max index at position 1 of Q you can access is 2 as you only have 2 rows, and you are trying to access element of 3rd row of Q which does not exist.

[Kevin Harnett]

I see the mistake. Thank you.