Warning: Matrix is singular to working precision. how to solve?

Hi, I have a problem with some matrixcomputations. I have to analyse an 8-barlinkage and have to find the angle, his velocity and acceleration. Calculating position and velocity doesn't give any problems, but for the acceleration I get a warning : "Warning: Matrix is singular to working precision. how to solve?". How do I get the results. The matrix computation for the velocity has the same form.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Kinematica en werkuigendynamica.
%
% Voorbeeldanalyse van een vierstangenmechanisme.
%
% Bram Demeulenaere <[e-mail]>
% Maarten De Munck <[e-mail]>
% Johan Rutgeerts <[e-mail]>
% Wim Meeussen <[e-mail]>
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [phi2, phi3, phi5, phi6gj, phi6fg, phi7, phi8, dphi2, dphi3, dphi5, dphi6gj, dphi6fg, dphi7, dphi8,...
    ddphi2, ddphi3, ddphi5, ddphi6gj, ddphi6fg, ddphi7, ddphi8] = ...
    kinematics_8bar(r2, r3bc, r3ce, r4, r5, r6fg, r6gj, r7, r8, ad, dg, gi,phifgj, phiad, phigd, phigi, phi4, dphi4, ddphi4,...
    phi2_init,phi3_init,phi5_init, phi6fg_init, phi7_init, phi8_init, t,fig_kin_8bar)

% allocation of the result vectors (this results in better performance because we don't have to reallocate and
% copy the vector each time we add an element.
phi2 = zeros(size(t));
phi3 = zeros(size(t));
phi5 = zeros(size(t));
phi6gj = zeros(size(t));
phi6fg = zeros(size(t));
phi7 = zeros(size(t));
phi8 = zeros(size(t));
dphi2 = zeros(size(t));
dphi3 = zeros(size(t));
dphi5 = zeros(size(t));
dphi6gj = zeros(size(t));
dphi6fg = zeros(size(t));
dphi7 = zeros(size(t));
dphi8 = zeros(size(t));
ddphi2 = zeros(size(t));
ddphi3 = zeros(size(t));
ddphi5 = zeros(size(t));
ddphi6gj = zeros(size(t));
ddphi6fg = zeros(size(t));
ddphi7 = zeros(size(t));
ddphi8 = zeros(size(t));

% fsolve options (help fsolve, help optimset)
optim_options = optimset('Display','off');

% *** loop over positions ***
Ts = t(2) - t(1);      % timestep
t_size = size(t,1);    % number of simulation steps
for k=1:t_size

      % *** position analysis ***

      % fsolve solves the non-linear set of equations
      % loop closure equations: see loop_closure_eqs.m
      % argument loop_closure_eqs: file containing closure equations
      % argument [..]': initial values of unknown angles phi3 and phi4
      % argument optim options: parameters for fsolve
      % argument phi2(k): input angle phi2 for which we want to calculate the unknown angles phi3 and phi4
      % argument a1 ... phi1: constants
      % return value x: solution for the unknown angles phi3 and phi4
      % return exitflag: indicates convergence of algorithm
      [x, fval, exitflag]=fsolve('loop_closure_eqs',[phi2_init,phi3_init,phi5_init, phi6fg_init, phi7_init, phi8_init]',...
          optim_options,phi4(k),r2, r3bc, r3ce, r4, r5, r6fg, r6gj, r7, r8, ad, dg, gi,phifgj, phiad, phigd, phigi);
      if (exitflag ~= 1)
          disp 'The fsolve exit flag was not 1, probably no convergence!'
          exitflag;
      end

      % save results of fsolve
      phi2(k)=x(1);
      phi3(k)=x(2);
      phi5(k)=x(3);
      phi6fg(k)=x(4);
      phi7(k)=x(5); 
      phi8(k)=x(6);
      phi6gj(k)= phi6fg(k)-phifgj;

      % *** velocity analysis ***

      A = [-r2*sin(phi2(k)), -r3bc*sin(phi3(k)), 0,                0,                    0,                0;
           r2*cos(phi2(k)),  r3bc*cos(phi3(k)),  0,                0,                    0,                0;
           0,                -r3ce*sin(phi3(k)), r5*sin(phi5(k)),  r6fg*sin(phi6fg(k)),  0,                0;
           0,                r3ce*cos(phi3(k)),  -r5*cos(phi5(k)), -r6fg*cos(phi6fg(k)), 0,                0;
           0,                0,                  0 ,               -r6gj*sin(phi6gj(k)), -r7*sin(phi7(k)), -r8*sin(phi8(k));
           0,                0,                  0 ,               r6gj*cos(phi6gj(k)),  r7*cos(phi7(k)),  r8*cos(phi8(k))];   
      B = [-r4*sin(phi4(k))*dphi4(k);
           r4*cos(phi4(k))*dphi4(k);
           r4*sin(phi4(k))*dphi4(k);
           -r4*cos(phi4(k))*dphi4(k); 
           0;
           0];

      x = A\B;

      % save results
      dphi2(k)=x(1);
      dphi3(k)=x(2);
      dphi5(k)=x(3);
      dphi6fg(k)=x(4);
      dphi7(k)=x(5);
      dphi8(k)=x(6);

      dphi6gj(k)=x(4);

      % *** acceleration analysis ***

      A = [-r2*sin(phi2(k)), -r3bc*sin(phi3(k)), 0,                0,                    0,                0;
           r2*cos(phi2(k)),  r3bc*cos(phi3(k)),  0,                0,                    0,                0;
           0,                -r3ce*sin(phi3(k)), r5*sin(phi5(k)),  r6fg*sin(phi6fg(k)),  0,                0;
           0,                r3ce*cos(phi3(k)),  -r5*cos(phi5(k)), -r6fg*sin(phi6fg(k)), 0,                0;
           0,                0,                  0,                -r6gj*sin(phi6gj(k)), -r7*sin(phi7(k)), -r8*sin(phi8(k));
           0,                0,                  0,                r6gj*sin(phi6gj(k)),  r7*sin(phi7(k)),  r8*sin(phi8(k))];

      B = [r2*cos(phi2(k))*dphi2(k)^2 + r3bc*cos(phi3(k))*dphi3(k)^2 - r4*sin(phi4(k))*ddphi4(k) - r4*cos(phi4(k))*dphi4(k)^2;
           r2*sin(phi2(k))*dphi2(k)^2 + r3bc*sin(phi3(k))*dphi3(k)^2 + r4*cos(phi4(k))*ddphi4(k) - r4*sin(phi4(k))*dphi4(k)^2; 
           r4*sin(phi4(k))*ddphi4(k) + r4*cos(phi4(k))*dphi4(k)^2 + r3ce*cos(phi3(k))*dphi3(k)^2 - r5*cos(phi5(k))*dphi5(k)^2 - r6fg*cos(phi6fg(k))*dphi6fg(k)^2;
           -r4*cos(phi4(k))*ddphi4(k) + r4*sin(phi4(k))*dphi4(k)^2 + r3ce*sin(phi3(k))*dphi3(k)^2 - r5*sin(phi5(k))*dphi5(k)^2 - r6fg*sin(phi6fg(k))*dphi6fg(k)^2;
           r6gj*cos(phi6gj(k))*dphi6gj(k)^2 + r7*cos(phi7(k))*dphi7(k)^2 + r8*cos(phi8(k))*dphi8(k)^2;
           r6gj*sin(phi6gj(k))*dphi6gj(k)^2 + r7*sin(phi7(k))*dphi7(k)^2 + r8*sin(phi8(k))*dphi8(k)^2];

      x = A\B;
      % save results
      ddphi2(k)=x(1);
      ddphi3(k)=x(2);
      ddphi5(k)=x(3);
      ddphi6fg(k)=x(4);
      ddphi7(k)=x(5);
      ddphi8(k)=x(6);

      ddphi6gj(k)=x(4);

      % *** calculate initial values for next iteration step ***
      phi2_init = phi2(k)+Ts*dphi2(k);
      phi3_init = phi3(k)+Ts*dphi3(k);
      phi5_init = phi5(k)+Ts*dphi5(k);
      phi6fg_init = phi6fg(k)+Ts*dphi6fg(k);
      phi7_init = phi7(k)+Ts*dphi7(k);
      phi8_init = phi8(k)+Ts*dphi8(k);

end % loop over positions

% *** create movie ***

% point D = fixed
D = 0;
% point S = fixed
A = ad * exp(j*phiad);
G = -dg * exp(1i*phigd);
I= G + gi * exp(j*phigi);

% define which positions we want as frames in our movie
frames = 200;    % number of frames in movie
delta = floor(t_size/frames); % time between frames
index_vec = [1:delta:t_size]';

% Create a window large enough for the whole mechanisme in all positions, to prevent scrolling.
% This is done by plotting a diagonal from (x_left, y_bottom) to (x_right, y_top), setting the
% axes equal and saving the axes into "movie_axes", so that "movie_axes" can be used for further
% plots.
x_left = -20;
y_bottom = -10;
x_right = 20;
y_top = 30;

figure(1)
hold on
plot([x_left, x_right], [y_bottom, y_top]);
axis equal;
movie_axes = axis;   %save current axes into movie_axes

% draw and save movie frame
for m=1:length(index_vec);
    index = index_vec(m);
    B = A + r2 * exp(j*phi2(index));
    C = D + r4 * exp(j*phi4(index));
    E = C + r3ce * exp(j*phi3(index));
    F = E - r5 * exp(j*phi5(index));
    J = G + r6gj * exp(j*phi6gj(index));
    H = I - r8 * exp(j*phi8(index)); 

      loop1 = [A B C D];
      loop2 = [D C E F G];
      loop3 = [G J H I];
      figure(1)
      clf
      hold on
      plot(real(loop1),imag(loop1),'-o',real(loop2),imag(loop2),'-o',real(loop3),imag(loop3),'-o')

      axis(movie_axes);     % set axes as in movie_axes
      Movie(m) = getframe;  % save frame to a variable Film
  end

% save movie
save fourbar_movie Movie
close(1)

% *** plot figures ***

if fig_kin_8bar

      %plot assembly at a certain timestep 
      index = 1; %select 1st timestep
      D = 0;
      A = ad * exp(j*phiad);
      G = -dg * exp(j*phigd);
      I= G + gi * exp(j*phigi);
      B = A + r2 * exp(j*phi2(index));
      C = D + r4 * exp(j*phi4(index));
      E = C + r3ce * exp(j*phi3(index));
      F = E - r5 * exp(j*phi5(index));
      J = G + r6gj * exp(j*phi6gj(index));
      H = I - r8 * exp(j*phi8(index)); 

      figure
      assembly1 = [A B C D];
      assembly2 = [D C E F G];
      assembly3 = [G J H I];
      plot(real(assembly1), imag(assembly1), 'ro-', real(assembly2), imag(assembly2), 'ro-',...
           real(assembly3), imag(assembly3), 'ro-')
      xlabel('[cm]')
      ylabel('[cm]')
      title('achtangenmechanisme')
      axis equal

      figure('Name','Hoeken','NumberTitle','off')
      subplot(421)
      plot(t,phi2)
      ylabel('\theta_2 [rad]')    
      subplot(422)
      plot(t,phi3)
      ylabel('\theta_3 [rad]')
      subplot(423)
      plot(t,phi4)
      ylabel('\theta_4 [rad]')
      subplot(424)
      plot(t,phi5)
      ylabel('\theta_5 [rad]')
      subplot(425)
      plot(t,phi6fg)
      ylabel('\theta_6_fg [rad]')
      subplot(426)
      plot(t,phi6gj)
      ylabel('\theta_6_gj [rad]')
      subplot(427)
      plot(t,phi7)
      ylabel('\theta_7 [rad]')
      xlabel('t [s]')
      subplot(428)
      plot(t,phi8)
      ylabel('\theta_8 [rad]')    
      xlabel('t [s]')

      figure('Name','Hoeksnelheid','NumberTitle','off')
      subplot(421)
      plot(t,dphi2)
      ylabel('\omega_2 [rad/s]')
      subplot(422)
      plot(t,dphi3)
      ylabel('\omega_3 [rad/s]')
      subplot(423)
      plot(t,dphi4)
      ylabel('\omega_4 [rad/s]')
      subplot(424)
      plot(t,dphi5)
      ylabel('\omega_5 [rad/s]')
      subplot(425)
      plot(t,dphi6fg)
      ylabel('\omega_6fg [rad/s]')
      subplot(426)
      plot(t,dphi6gj)
      ylabel('\omega_6gj [rad/s]')
      subplot(427)
      plot(t,dphi7)
      ylabel('\omega_7 [rad/s]')
      subplot(428)
      plot(t,dphi8)
      ylabel('\omega_8 [rad/s]')
      xlabel('t [s]')

      figure('Name','Hoekversnelling','NumberTitle','off')
      subplot(421)
      plot(t,ddphi2)
      ylabel('\alpha_2 [rad/^s]')
      subplot(422)
      plot(t,ddphi3)
      ylabel('\alpha_3 [rad/^s]')
      subplot(423)
      plot(t,ddphi4)
      ylabel('\alpha_4 [rad/^s]')
      subplot(424)
      plot(t,ddphi5)
      ylabel('\alpha_5 [rad/^s]')
      subplot(425)
      plot(t,ddphi6fg)
      ylabel('\alpha_6fg [rad/^s]')
      subplot(426)
      plot(t,ddphi6gj)
      ylabel('\alpha_6gj [rad/^s]')
      subplot(427)
      plot(t,ddphi7)
      ylabel('\alpha_7 [rad/^s]')
      subplot(428)
      plot(t,ddphi8)
      ylabel('\alpha_8 [rad/^s]')
      xlabel('t [s]')

end

Answer 1

[njj1]

There are ways of getting around a singular matrix, but they are a little ad hoc. have you tried QR decomposition ([Q,R] = qr(A))? You can also try adding small values along the diagonal.