MATLAB: The final B matrix should be a zero matrix but i’m getting very small values in it, despite the code being correct; please help.

incorrect solutionswriong answers

%Properties = [E1 E2 G12 v12 v23]
Properties = [155 12.1 4.4 0.248 0.458];
X = [0 90 90 0];
h = 150*10^-6;
z = length(X);
H = h*z;
Z =[];
ex = 10^-3;
ey = 0;
exy = 0;
kx = 0;
ky = 0;
kxy = 0;
for i = -z/2:z/2
Z = [Z,i*h];
end
Z2 = Z.^2;
E1 = Properties(1)*10^9;
E2 = Properties(2)*10^9;
v12 = Properties(4);
G12 = Properties(3)*10^9;
%On-Axis Compliance and Stiffness matrix
onS11 = 1/E1;
onS12 = -v12/E1;
onS22 = 1/E2;
onS66 = 1/G12;
onS = [onS11 onS12 0;onS12 onS22 0;0 0 onS66];
disp('On-Axis Compliance Matrix (S)(1/Gpa)= \n')
disp(onS)
qden = onS11*onS22 – onS12*onS12;
onQ11 = onS22/qden;
onQ12 = -onS12/qden;
onQ11 = onS22/qden;
onQ22 = onS11/qden;
onQ66 = 1/onS66;
onQ = [onQ11 onQ12 0;onQ12 onQ22 0;0 0 onQ66];
disp('On-Axis Stiffness Matrix (Q) (Gpa): \n')
disp(onQ)
S= [];
Q= [];
%Compliance and Stiffness matrix of each layer
for i =1:length(X)
m=cosd(X(i));
n=sind(X(i));
S11 = (onS11*m^4)+(2*onS12+onS66)*(n^2)*(m^2)+onS22*n^4;
S12 = (onS11+onS22-onS66)*(n^2)*(m^2)+ onS12*((n^4)+(m^4));
S16 = (2*onS11-2*onS12-onS66)*n*m^3 – (2*onS22-2*onS12-onS66)*m*n^3;
S22 = (onS11*n^4)+(2*onS12+onS66)*(n^2)*(m^2)+onS22*m^4;
S26 = (2*onS11-2*onS12-onS66)*m*n^3 – (2*onS22-2*onS12-onS66)*n*m^3;
S66 = 2*(2*onS11+2*onS22-4*onS12-onS66)*(n^2)*(m^2)+onS66*((n^4)+(m^4));
S(:,:,i)= [S11 S12 S16;S12 S22 S26;S16,S26,S66];
Q11 = (onQ11*m^4)+2*(onQ12+2*onQ66)*(n^2)*(m^2)+onQ22*n^4;
Q12 = (onQ11+onQ22-4*onQ66)*(n^2)*(m^2)+ onQ12*((n^4)+(m^4));
Q16 = (onQ11-onQ12-2*onQ66)*n*m^3 + (onQ12-onQ22+2*onQ66)*m*n^3;
Q22 = (onQ11*n^4)+2*(onQ12+2*onQ66)*(n^2)*(m^2)+onQ22*m^4;
Q26 = (onQ11-onQ12-2*onQ66)*m*n^3 + (onQ12-onQ22+2*onQ66)*n*m^3;
Q66 = (onQ11+onQ22-2*onQ12-2*onQ66)*(n^2)*(m^2)+onQ66*((n^4)+(m^4));
Q(:,:,i)= [Q11,Q12,Q16;Q12,Q22,Q26;Q16,Q26,Q66];
end
%A,B,D Matrices
A= [0 0 0;0 0 0;0 0 0];
B= [0 0 0;0 0 0;0 0 0];
D= [0 0 0;0 0 0;0 0 0];
p=0;
q=0;
r=0;
for i=1:3
for j=1:3
for k=1:z
p=Q(i,j,k)*(Z(k+1)-Z(k));
q=Q(i,j,k)*(Z2(k+1)-Z2(k));
r=Q(i,j,k)*((Z(k+1))^3-(Z(k))^3);
A(i,j)= A(i,j)+p;
B(i,j)= B(i,j)+q;
D(i,j)= D(i,j)+r;
end
end
end
B = B/2;
D = D/3;MathWorks

Best Answer

You are seeing the cumulative results of floating-point approximation error.
See the documentation section on Accuracy of Floating-Point Data (link) for a relevant discussion.