Error finding the natural frequency of a cantilever beam

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);

Answer 1

[KSSV]

YOu may refer this link:

http://in.mathworks.com/matlabcentral/fileexchange/30970-natural-frequencies---buckling-loads-of-columns-using-finite-element-method

[Pedro Pires]

Thank you for the help. I will try to adapt my code based on what you made.