MATLAB: Error finding the natural frequency of a cantilever beam

beamcantilevereiginvMATLABmode shapenatural frequency

Hi, I currently have an assignment where I have to find the 6 first natural frequencies of a cantilever beam. I've writen the code but I get the error shown bellow and would like to ask for some hints on the matter. I started by finding the stifness and mass matrix and then turn them to the global matrixes. I don't know however how to get the reduced matrixes and I think that's the problem with code although I'm not sure if the eig of the inverted matrix will give me what I need. The cross section is a hollow square and the thickness reduces along the lenght Thank you
Warning: Matrix is singular to working precision.
> In trab2 (line 228) %T=inv(mred)*kred;
Error using eig
Input to EIG must not contain NaN or Inf.
Error in trab2 (line 230) [V,D]=eig(T);
E = 70e+6;                    %[Pa] Young module
L = 5;                        %[m] length
h= 0.12;                      %[m] cross section width
b = 0.12;                     %[m] cross section hight
nel = 32;                     %nr of elements 
t_0 = 0.005;                  %inital thickness
t_1 = 0.0003;                 %final thickness
rho=2700;                     %density
%(distance between nodes)
dt = L/nel; 
% Cálculo de numero de nós (nr of nodes)
n_node = nel+1;
% position of the nodes
for i=1:n_node;
    xe_coord(i) = (i-1)*dt;
end
% Cálculo da espessura da viga (beam thickness)
for i=1:nel+1;
    t(i) = t_0 - t_1*xe_coord(i);
end
%(element thickness)
for i=1:nel;
    xe_t(i)=(t(i)+t(i+1))/2;
end
% (Inercia)
for i=1:nel;
    I_el(i)=((b*h^3)/12.0)-((b-2*xe_t(i))*(h-2*xe_t(i))^3)/12.0;
end
for i=1:nel;
    A(i)=b*h-(b-2*t(i))*(h-2*t(i));         %Area of cross section
end
% (Stifness matrix)
for i = 1:nel;
      k_el(1,1,i) = I_el(i)*6.0;
      k_el(1,2,i) = I_el(i)*(-3.0*dt);
      k_el(1,3,i) = I_el(i)*(-6.0);
      k_el(1,4,i) = I_el(i)*(-3.0*dt);
      k_el(2,1,i) = I_el(i)*(-3.0*dt);
      k_el(2,2,i) = I_el(i)*(2.0*dt^2);
      k_el(2,3,i) = I_el(i)*(3.0*dt);
      k_el(2,4,i) = I_el(i)*(dt^2);
      k_el(3,1,i) = I_el(i)*(-6.0);
      k_el(3,2,i) = I_el(i)*(3.0*dt);
      k_el(3,3,i) = I_el(i)*6.0;
      k_el(3,4,i) = I_el(i)*(3.0*dt);
      k_el(4,1,i) = I_el(i)*(-3.0*dt);
      k_el(4,2,i) = I_el(i)*(dt^2);
      k_el(4,3,i) = I_el(i)*(3.0*dt);
      k_el(4,4,i) = I_el(i)*(2.0*dt^2);
end
%Matriz Rigidez de um elemento finito unidimensional arbitrário  (
k_el=2*E/dt^3*k_el;
% CALCULO DA MATRIZ DE MASSA DOS ELEMENTOS (Mass Matrix)
for i = 1:nel;
      m_el(1,1,i) = A(i)*156;
      m_el(1,2,i) = A(i)*(-22)*dt;
      m_el(1,3,i) = A(i)*54;
      m_el(1,4,i) = A(i)*13*dt;
      m_el(2,1,i) = A(i)*(-22*dt);
      m_el(2,2,i) = A(i)*4*dt^2;
      m_el(2,3,i) = A(i)*(-13)*dt;
      m_el(2,4,i) = A(i)*(-3)*dt^2;
      m_el(3,1,i) = A(i)*54;
      m_el(3,2,i) = A(i)*(-13)*dt;
      m_el(3,3,i) = A(i)*156;
      m_el(3,4,i) = A(i)*4*dt^2;
      m_el(4,1,i) = A(i)*13*dt;
      m_el(4,2,i) = A(i)*(-3)*dt^2;
      m_el(4,3,i) = A(i)*22*dt;
      m_el(4,4,i) = A(i)*4*dt^2;
end
m_el=rho*dt/420*m_el;
%Matrizes Globais (Global Matrixes)
k_glob=zeros(2*nel+2,2*nel+2);
for i=1:nel;
      k_glob(2*i-1,2*i-1) = k_glob(2*i-1,2*i-1)+k_el(1,1,i);
      k_glob(2*i-1,2*i) = k_glob(2*i-1,2*i)+k_el(1,2,i);
      k_glob(2*i-1,2*i+1) = k_glob(2*i-1,2*i+1)+k_el(1,3,i);
      k_glob(2*i-1,2*i+2) = k_glob(2*i-1,2*i+2)+k_el(1,4,i);
      k_glob(2*i,2*i-1) = k_glob(2*i,2*i-1)+k_el(2,1,i);
      k_glob(2*i,2*i) = k_glob(2*i,2*i)+k_el(2,2,i);
      k_glob(2*i,2*i+1) = k_glob(2*i,2*i+1)+k_el(2,3,i);
      k_glob(2*i-2*i+2) = k_glob(2*i,2*i+2)+k_el(2,4,i);
      k_glob(2*i+1,2*i-1) = k_glob(2*i+1,2*i-1)+k_el(3,1,i);
      k_glob(2*i+1,2*i) = k_glob(2*i+1,2*i)+k_el(3,2,i);
      k_glob(2*i+1,2*i+1) = k_glob(2*i+1,2*i+1)+k_el(3,3,i);
      k_glob(2*i+1,2*i+2) = k_glob(2*i+1,2*i+2)+k_el(3,4,i);
      k_glob(2*i+2,2*i-1) = k_glob(2*i+2,2*i-1)+k_el(4,1,i);
      k_glob(2*i+2,2*i) = k_glob(2*i+2,2*i)+k_el(4,2,i);
      k_glob(2*i+2,2*i+1) = k_glob(2*i+2,2*i+1)+k_el(4,3,i);
      k_glob(2*i+2,2*i+2) = k_glob(2*i+2,2*i+2)+k_el(4,4,i);
end
m_glob=zeros(2*nel+2,2*nel+2);
for i=1:nel;
      m_glob(2*i-1,2*i-1) = m_glob(2*i-1,2*i-1)+m_el(1,1,i);
      m_glob(2*i-1,2*i) = m_glob(2*i-1,2*i)+m_el(1,2,i);
      m_glob(2*i-1,2*i+1) = m_glob(2*i-1,2*i+1)+m_el(1,3,i);
      m_glob(2*i-1,2*i+2) = m_glob(2*i-1,2*i+2)+m_el(1,4,i);
      m_glob(2*i,2*i-1) = m_glob(2*i,2*i-1)+m_el(2,1,i);
      m_glob(2*i,2*i) = m_glob(2*i,2*i)+m_el(2,2,i);
      m_glob(2*i,2*i+1) = m_glob(2*i,2*i+1)+m_el(2,3,i);
      m_glob(2*i-2*i+2) = m_glob(2*i,2*i+2)+m_el(2,4,i);
      m_glob(2*i+1,2*i-1) = m_glob(2*i+1,2*i-1)+m_el(3,1,i);
      m_glob(2*i+1,2*i) = m_glob(2*i+1,2*i)+m_el(3,2,i);
      m_glob(2*i+1,2*i+1) = m_glob(2*i+1,2*i+1)+m_el(3,3,i);
      m_glob(2*i+1,2*i+2) = m_glob(2*i+1,2*i+2)+m_el(3,4,i);
      m_glob(2*i+2,2*i-1) = m_glob(2*i+2,2*i-1)+m_el(4,1,i);
      m_glob(2*i+2,2*i) = m_glob(2*i+2,2*i)+m_el(4,2,i);
      m_glob(2*i+2,2*i+1) = m_glob(2*i+2,2*i+1)+m_el(4,3,i);
      m_glob(2*i+2,2*i+2) = m_glob(2*i+2,2*i+2)+m_el(4,4,i);
end
%reduced matrixes
kred=zeros(2*nel,2*nel);
for i=1:2*nel;
    for j=1,2*nel;
        kred(i,j) = k_glob(i+2,j+2);
    end
end
mred=zeros(2*nel,2*nel);
for i=1:2*nel;
    for j=1,2*nel;
        mred(i,j) = m_glob(i+2,j+2);
    end
end
T=inv(mred)*kred;
[V,D]=eig(T);
Related Question