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

17 次查看（过去 30 天）
louis snyders
louis snyders 2018-3-16
回答： njj1 2018-3-16
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 <bram.demeulenaere@mech.kuleuven.be>
% Maarten De Munck <maarten.demunck@mech.kuleuven.be>
% Johan Rutgeerts <johan.rutgeerts@mech.kuleuven.be>
% Wim Meeussen <wim.meeussen@mech.kuleuven.be>
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

请先登录，再回答此问题。

回答（1 个）

njj1
njj1 2018-3-16
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.

类别

在 Help Center 和 File Exchange 中查找有关 Data Distribution Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by