How to Solve this Big System of Equations

Dear Community,

I have made a stiffness matrix of 443x443 dimension and I have an associated load vector of 443x1. Since the matrix is diagonal, how would you solve this system? Attached the script:

The system was created as follows: [K]{x} = {L} and I need the {x} vector

Thanks!

%% 1st Order Cylindrical Segment

ABD_cyl_13 =  [157142.857142857,47142.8571428571,    0,               0,               0,0;
               47142.8571428571,157142.857142857,    0,               0,               0,0;
                              0,               0,55000,               0,               0,0;
                              0,               0,    0,22130.9523809524,6639.28571428571,0;
                              0,               0,    0,6639.28571428571,22130.9523809524,0;
                              0,               0,    0,               0,               0,7745.83333333333];
r_cyl = 608.9/2;

% det_11 = - ABD_cyl_13(1,1)* ABD_cyl_13(4,4) + ABD_cyl_13(1,4)^2;
% a11_12 = - ABD_cyl_13(4,4)*ABD_cyl_13(1,2)+ABD_cyl_13(1,4)*ABD_cyl_13(1,5);
% a11_32 = ABD_cyl_13(1,1) * ABD_cyl_13(1,5) - ABD_cyl_13(1,4) * ABD_cyl_13(1,2);
% a11_54 = ABD_cyl_13(1,4) * ABD_cyl_13(2,4) - ABD_cyl_13(4,4) * ABD_cyl_13(2,1);
% a11_56 = ABD_cyl_13(2,1) * ABD_cyl_13(1,4) - ABD_cyl_13(2,4) * ABD_cyl_13(1,1);
% a11_52 = -ABD_cyl_13(2,4) * ABD_cyl_13(4,4) *ABD_cyl_13(1,2) + ABD_cyl_13(2,4) * ABD_cyl_13(1,4) * ABD_cyl_13(1,5) - ABD_cyl_13(5,4) * ABD_cyl_13(1,1) *ABD_cyl_13(1,5) + ABD_cyl_13(5,4) * ABD_cyl_13(1,4) * ABD_cyl_13(1,2);

A_cyl_131 = [      0,                        1/r_cyl - ABD_cyl_13(1,2)/(ABD_cyl_13(1,1)*r_cyl),0,            1/(r_cyl*ABD_cyl_13(1,1)),      0,                        0,                  0;
             -1/r_cyl,                                                                        0,1,                                    0,      0,                        0,                  0;
                   0,                                                                        0,0,                                    0,      0,1/(r_cyl*ABD_cyl_13(4,4)),                  0;
                   0,                                                                        0,0,                                    0,1/r_cyl,                        0,                  0;
                   0,1/r_cyl*(ABD_cyl_13(2,2)-ABD_cyl_13(1,2)*ABD_cyl_13(2,1)/ABD_cyl_13(1,1)),0,- (1/r_cyl) + ABD_cyl_13(2,1)/(r_cyl),      0,                        0,-1.892*(1-0.5*0.33);%-1.892
                   0,                                                                        0,0,                                    0,     -1,                        0,                  0;
                   0,                                                                        0,0,                                    0,      0,                        0,                  0];
T1{1,:} = eye(7,7);
for i = 1:60-1
    if i==1
    T1{2,:} = A_cyl_131*T1{1,:}  ;
    elseif i == 2
    T1{i+1,:} = A_cyl_131*T1{2,:} + A_cyl_131*T1{1,:};
    elseif i == 3
    T1{i+1,:} = A_cyl_131*T1{i,:} + A_cyl_131*T1{i-1,:} + A_cyl_131*T1{i-2,:};
    elseif i > 3
    T1{i+1,:} = A_cyl_131*T1{i,:} + A_cyl_131*T1{i-1,:} + A_cyl_131*T1{i-2,:} + A_cyl_131*T1{i-3,:} ;
    end
end 

%% Geometrical Aspects - 1.3mm Segment
Step1 = 300/59; % Definition of the step - 1,1 mm thickness skirt 100
x2 = 0:Step1:300; % Extension of the 1,1 mm thickness skirt
Segment_11 = 300;
%Step11 = 300/5900;
L1 = length(x2);
B_13 = cell(1,59);                                % initialzation
B_13{:,1} = eye(7,7);
for tal = 1:60
    B1{tal,:} = eye(7,7) + A_cyl_131 *(300/59) + (A_cyl_131^2*(300/59)^2)/4 + (A_cyl_131^3*(300/59)^3)/6;
    
end
TMM_cyl_1 = cell2mat(T1);
% TMM_cyl_1 = cell2mat(B1);
TMM_cyl_1(7:7:end,:) = [];
Load_cyl = TMM_cyl_1(:,7);
TMM_cyl_1(:,7) = [];
Load_cyl_1 = mat2cell(Load_cyl,6*ones(60,1),ones(1,1));
B_matrix_cyl_1 = mat2cell(TMM_cyl_1,6*ones(60,1),6*ones(1,1));
for jj=1:60
    Z{jj,:} = pinv(B_matrix_cyl_1{jj,:})*Load_cyl_1{jj,:};
    
end
%% 1st Order Dome Segment

ABD_dome_065 = [80237.9082033442,26478.5097071036,0,26077.3201660869,8605.51565480867,0;26478.5097071036,80237.9082033442,0,8605.51565480867,26077.3201660869,0;0,0,26879.6992481203,0,0,8735.90225563910;26077.3201660869,8605.51565480867,0,11300.1720719710,3729.05678375042,0;8605.51565480867,26077.3201660869,0,3729.05678375042,11300.1720719710,0;0,0,8735.90225563910,0,0,3785.55764411028];
%% Spherical Dome - 2mm thickness
%ABD_dome_2 = zeros(7,7);
% a = (ABD_dome_2(1,1)*ABD_dome_2(1,5) - ABD_dome_2(1,4)* ABD_dome_2(1,2));
% a2 = (ABD_dome_2(1,1)* ABD_dome_2(1,5) - ABD_dome_2(1,4)*ABD_dome_2(1,2));
% a4 = (- ABD_dome_2(2,1)* ABD_dome_2(4,4)* ABD_dome_2(1,2)+ ABD_dome_2(2,1)* ABD_dome_2(1,4)* ABD_dome_2(1,5)- ABD_dome_2(5,4)* ABD_dome_2(1,1)* ABD_dome_2(1,5)+ ABD_dome_2(2,4)*ABD_dome_2(1,4)*ABD_dome_2(1,2));
a11= - ABD_dome_065(1,2) * ABD_dome_065(4,4) + ABD_dome_065(1,5) * ABD_dome_065(1,4);
% a1 = - ABD_dome_2(1,2) * ABD_dome_2(4,4) + ABD_dome_2(2,4) * ABD_dome_2(1,4);
a12 = - ABD_dome_065(1,5)* ABD_dome_065(4,4) + ABD_dome_065(1,4) * ABD_dome_065(4,5);
a31 = - ABD_dome_065(1,2)* ABD_dome_065(1,4) + ABD_dome_065(1,5) * ABD_dome_065(1,1);
a33 = - ABD_dome_065(1,4)* ABD_dome_065(1,5) + ABD_dome_065(4,5) * ABD_dome_065(1,1);
det_dome = - ABD_dome_065(1,1)* ABD_dome_065(4,4) + ABD_dome_065(1,4) * ABD_dome_065(1,4);
a41 = - ABD_dome_065(2,1)*ABD_dome_065(4,4)*ABD_dome_065(1,2) + ABD_dome_065(2,1)*ABD_dome_065(1,4)*ABD_dome_065(1,5) + ABD_dome_065(2,4)*ABD_dome_065(1,4)*ABD_dome_065(1,2) - ABD_dome_065(2,4)*ABD_dome_065(1,5)*ABD_dome_065(1,1);
%a43 = - ABD_dome_2(2,1)*ABD_dome_2(4,4)*ABD_dome_2(1,5) + ABD_dome_2(2,1)*ABD_dome_2(1,4)*ABD_dome_2(1,5) + ABD_dome_2(2,4)*ABD_dome_2(1,4)*ABD_dome_2(1,2) - ABD_dome_2(2,4)*ABD_dome_2(1,5)*ABD_dome_2(1,1);
a46 = - ABD_dome_065(2,1)* ABD_dome_065(1,4) + ABD_dome_065(2,4)* ABD_dome_065(1,1);
a44 = - ABD_dome_065(2,1)* ABD_dome_065(4,4) + ABD_dome_065(2,4)* ABD_dome_065(1,1);
a61 = - ABD_dome_065(2,4) * ABD_dome_065(4,4) * ABD_dome_065(1,2) + ABD_dome_065(2,4) * ABD_dome_065(1,4) *  ABD_dome_065(1,5) + ABD_dome_065(5,4) * ABD_dome_065(1,4) * ABD_dome_065(1,2) - ABD_dome_065(5,4) * ABD_dome_065(1,1) * ABD_dome_065(1,5);
% a62 = - ABD_dome_2(2,4) * ABD_dome_2(4,4) * ABD_dome_2(1,2) + ABD_dome_2(2,4) * ABD_dome_2(1,4) *  ABD_dome_2(1,5) + ABD_dome_2(5,4) * ABD_dome_2(1,4) * ABD_dome_2(1,2) - ABD_dome_2(5,4) * ABD_dome_2(1,1) * ABD_dome_2(1,5);
a66 =  ABD_dome_065(2,4)* ABD_dome_065(1,4) - ABD_dome_065(5,4)* ABD_dome_065(1,1);
a64 = - ABD_dome_065(2,4)* ABD_dome_065(4,4) + ABD_dome_065(5,4)* ABD_dome_065(1,4);
a63 = - ABD_dome_065(2,4)*ABD_dome_065(4,4)*ABD_dome_065(1,5) + ABD_dome_065(4,5)*ABD_dome_065(2,4)*ABD_dome_065(1,4) + ABD_dome_065(5,4)*ABD_dome_065(1,4)*ABD_dome_065(1,5) - ABD_dome_065(5,4)*ABD_dome_065(1,1)*ABD_dome_065(4,5);

R = 304.5;
A_sp = cell(65,1);
A_sp{1,:} = zeros(7,7);

%% In case Bij not 0
for phi = 1:89
    
    phi_rad = (phi-1) * pi/180; %angle in radians 5*phi/304.5;
    r = 304.5 / cos(phi_rad);
    A_sp{phi,:} = [- sin(phi_rad)* a11/(det_dome*r),1/304.5 - (cos(phi_rad) * a11)/(det_dome*r),(sin(phi_rad) * a12)/(det_dome*r),-ABD_dome_065(4,4)/(det_dome*r),0,-ABD_dome_065(1,4)/(det_dome*r),0;
        -1/R,0,1,0,0,0,0;
        -(sin(phi_rad) * a31)/(r * det_dome),-(cos(phi_rad) * a31)/(r * det_dome),(sin(phi_rad)*a33)/(r*det_dome),-ABD_dome_065(1,4)/(det_dome*r),0,-ABD_dome_065(1,1)/(det_dome*r),0;
        (sin(phi_rad)^2/r)*(ABD_dome_065(2,2) - (a41)/det_dome) , (cos(phi_rad)*sin(phi_rad)/r)*(ABD_dome_065(2,2) - (a41)/det_dome) ,(sin(phi_rad)^2/r) *  (ABD_dome_065(2,5) - a41/det_dome), ((sin(phi_rad)*ABD_dome_065(2,1))/(ABD_dome_065(1,1)*r)) , 1/R ,0,0;
        (cos(phi_rad)*sin(phi_rad)/r)*(ABD_dome_065(2,2) - (a41)/det_dome),(cos(phi_rad)^2/r)*(ABD_dome_065(2,2) - (a41)/det_dome), (sin(phi_rad)*cos(phi_rad)/r) *  (ABD_dome_065(2,5) - a41/det_dome) ,- 1/R + (cos(phi_rad) * a44)/(det_dome*r) ,0,(-cos(phi_rad) * a46)/(det_dome*r),-1.892*(0.5-0.5*0.33);%-1.892
        (sin(phi_rad)^2/r)*(-ABD_dome_065(2,5) + a61 /(det_dome)),    (-sin(phi_rad)*cos(phi_rad)/r)*(ABD_dome_065(2,5)/r - a61/(det_dome)),   (sin(phi_rad)^2/R) * (ABD_dome_065(5,5) -  a63/det_dome)  ,  (-sin(phi_rad) * a64)/(det_dome*r),-1,(sin(phi_rad) * a66 )/(det_dome*r),0;
        0,0,0,0,0,0,0];
    
end
% A_sp = A_sp(25:end,:);

%% Dome - 0.65mm Thickness
Step2 = 304.5*pi/180;
theta = 1:89;
for tah = 1:1:89
    B_sp1{tah,:} = (eye(7,7) + A_sp{tah,:} *(Step2) + (A_sp{tah,:}^2*(Step2)^2)/4 + (A_sp{tah,:}^3*(Step2)^3)/6 + (A_sp{tah,:}^4*(Step2)^4)/24);
%     B_sp1{tah,:} = (eye(7,7) + A_sp{tah,:}*Step2);
end

T{1,:} = eye(7,7);

for i = 1:89-1
    if i==1
    T{2,:} = A_sp{1,:}*T{1,:}  ;
    elseif i == 2
    T{i+1,:} = A_sp{1,:}*T{2,:} + A_sp{2,:}*T{1,:};
    elseif i == 3
    T{i+1,:} = A_sp{1,:}*T{i,:} + A_sp{2,:}*T{i-1,:} + A_sp{3,:}*T{i-2,:};
    elseif i > 3
    T{i+1,:} = A_sp{1,:}*T{i,:} + A_sp{2,:}*T{i-1,:} + A_sp{3,:}*T{i-2,:} + A_sp{4,:}*T{i-3,:} ;
    end
end 
TMM_dome_1 = cell2mat(T);
% TMM_dome_1 = cell2mat(B_sp1);
TMM_dome_1(7:7:end,:) = [];
Load_dome = TMM_dome_1(:,7);
TMM_dome_1(:,7) = [];
Load_dome_1 = mat2cell(Load_dome,6*ones(89,1),ones(1,1));
B_matrix_dome_1 = mat2cell(TMM_dome_1,6*ones(89,1),6*ones(1,1));
for kK = 1:89
    theta_transformation = theta(kK)* pi / 180;
    Transformation_Matrix{kK,:} = [cos(theta_transformation),-sin(theta_transformation),0,0,0,0;
                                   sin(theta_transformation),cos(theta_transformation),0,0,0,0;
                                                           0,                        0,1,0,0,0;
                                   0,0,0,cos(theta_transformation),-sin(theta_transformation),0;
                                   0,0,0,sin(theta_transformation),cos(theta_transformation),0;
                                                                                   0,0,0,0,0,1];
                               
                               %% Transformation of the Transfer Matrix
                               
                               TM_transformed{kK,:} = Transformation_Matrix{kK,:}*B_matrix_dome_1{kK,:}*transpose(Transformation_Matrix{kK,:});
end
 B_global = [B_matrix_cyl_1;B_matrix_dome_1];
Load_general_vector = [Load_cyl_1;Load_dome_1];
Load_System = mat2cell(cell2mat(Load_general_vector),3*ones(298,1),ones(1,1));
Gesamtsystem_13 = mat2cell(cell2mat(B_global), 3 * ones(1,298), 3 * ones(1,2));
for nn=1:2:297
    KS{nn,1} = pinv(Gesamtsystem_13{nn,2}) * Gesamtsystem_13{nn,1};
    KS{nn,2} = pinv(-Gesamtsystem_13{nn,2});
end
for nn=2:2:298
%     KS{nn,1} = Gesamtsystem_13{nn,1} - Gesamtsystem_13{nn,2}*pinv(Gesamtsystem_13{nn-1,2})*Gesamtsystem_13{nn-1,1};
    KS{nn,1} = transpose(pinv(-Gesamtsystem_13{nn-1,2}));%transpose
    KS{nn,2} = (Gesamtsystem_13{nn,2}) * pinv(Gesamtsystem_13{nn-1,2});
end
for mmm = 1:2:297
    
    LOAD{mmm,1} = pinv(- Gesamtsystem_13{mmm,2}) * Load_System {mmm,:}; % General Local Load vectors
end

for mmm = 2:2:298
    
    LOAD{mmm,1} = - Load_System {mmm,1} + Gesamtsystem_13{mmm,2} * pinv(Gesamtsystem_13{mmm-1,2}) * Load_System {mmm-1,1}; % General Local Load vectors
end


for i = 1:149
    for j = 1:149
        
        if (i==1) && (i==j)
            
            Structure_Global_Stiffness1{i,i} = KS{2*i,2};
            Structure_Global_Stiffness1{i,i+1} = zeros(3);
            
        elseif  (i==j) && (i>1) && (i<149)
            Structure_Global_Stiffness1{i,i-1} =  KS{2*i,1};
            Structure_Global_Stiffness1{i,i+1} =  KS{2*i+1,2};
         
        elseif (i==j) && (i==149)
            Structure_Global_Stiffness1{i,i-1} =  KS{2*i,1};
            Structure_Global_Stiffness1{i,i} =  KS{2*i-2,2};
        elseif (j < i-1)
            Structure_Global_Stiffness1{i,j} = zeros(3); 
        elseif (i <149) && (j > i+1)
            Structure_Global_Stiffness1{i,j} = zeros(3);
        elseif (i==149) && (j > i+2)
            Structure_Global_Stiffness1{i,j} = zeros(3);
        end
    end
end

for i = 2:148
    for j = 2:148
        
        if (i==j) && (i == 1)
         Structure_Global_Stiffness1{i,j} = KS{2,2} + KS{3,1};
        elseif (i == j) && (i > 1)
         Structure_Global_Stiffness1{i,i} = KS{2*i,2} + KS{2*i+1,1}; 
        end
    end
end 


Stiffness_Matrix = cell2mat(Structure_Global_Stiffness1);

   for k = 1:149
    if (k==1)
        Load_System_a1{k,:} = LOAD{2*k,:};
    elseif (k > 1) && (k < 149)
        Load_System_a1{k,:} = LOAD{2*k,:} + LOAD{2*k+1,:};
    elseif (k == 149)
        Load_System_a1{k,:} = LOAD{2*k-1,:};
    end
   end
LOAD_VECTOR_SYSTEM_a1 = cell2mat(Load_System_a1);
   

% Stiffness_Matrix(:,447) = []; % Moment
% Stiffness_Matrix(:,446) = []; % Shear
Stiffness_Matrix(:,445) = []; % Axial
% Stiffness_Matrix(:,444) = []; %Rotation
% Stiffness_Matrix(:,443) = []; % w displacement
Stiffness_Matrix(:,442) = [];   % horizontal displacement
% Stiffness_Matrix(:,6) = []; % Momentum
% Stiffness_Matrix(:,5) = []; % Shear
Stiffness_Matrix(:,4) = []; % Axial
% Stiffness_Matrix(:,3) = []; % Rotation
% Stiffness_Matrix(:,2) = []; % w displacement
Stiffness_Matrix(:,1) = []; % u displacement

% Stiffness_Matrix(447,:) = [];
% Stiffness_Matrix(446,:) = [];
Stiffness_Matrix(445,:) = [];
% Stiffness_Matrix(444,:) = [];
% Stiffness_Matrix(443,:) = [];
Stiffness_Matrix(442,:) = [];
% Stiffness_Matrix(6,:) = [];
% Stiffness_Matrix(5,:) = [];
Stiffness_Matrix(4,:) = [];
% Stiffness_Matrix(3,:) = [];
% Stiffness_Matrix(2,:) = [];
Stiffness_Matrix(1,:) = [];

BBBB = eye(443);

%% Boundary Conditions in the Load Vector

% LOAD_VECTOR_SYSTEM_a1(447,:) = [];
% LOAD_VECTOR_SYSTEM_a1(446,:) = [];
LOAD_VECTOR_SYSTEM_a1(445,:) = [];
% LOAD_VECTOR_SYSTEM_a1(444,:) = [];
% LOAD_VECTOR_SYSTEM_a1(443,:) = [];
LOAD_VECTOR_SYSTEM_a1(442,:) = [];
% LOAD_VECTOR_SYSTEM_a1(6,:) = [];
% LOAD_VECTOR_SYSTEM_a1(5,:) = [];
LOAD_VECTOR_SYSTEM_a1(4,:) = [];
% LOAD_VECTOR_SYSTEM_a1(3,:) = [];
% LOAD_VECTOR_SYSTEM_a1(2,:) = [];
LOAD_VECTOR_SYSTEM_a1(1,:) = [];
   

[Marcelo Boldt]

Would you do {x} = pinv([K])*{L} ?

Answer 1

[Bruno Luong]

Possible causes:

• your matrix is wrongly constructed
• you forget to include proper boundary conditions
• you forget to add conditions that make your system well posed (such as free solid displacement in 3D mechanics)
• the EDP you solve is not well posed (unlikely in forward mechanical problem), you solve a physical ill-posed physical system.

None of that we can help you since it is not direcly MATLAB related.

Answer 2

[Anna Lischke]

It looks as though your Stiffness_Matrix contains several rows that are all zero. In this case, the matrix is not invertible, and attempting to invert this matrix in MATLAB will result in a warning telling you that your matrix is 'singular to working precision'.

I suggest you review your code that assembles your matrix for correctness.

You can use the nnz function to count how many nonzero values occur in an array. For example, you can see that the first row of Stiffness_Matrix is all zeros by running the following command:

nnz(Stiffness_Matrix(1,:))

The result will be zero, meaning there are no nonzero values in the first row of Stiffness_Matrix.

Once you have an invertible matrix K, you can solve the equation Kx = L for the vector x using mldivide (i.e., 'backslash'):

x = K\L

[Marcelo Boldt]

Hello Anna,

Thank you for your reply, I have just use the nnz command and now I have an invertible matrix but not matlab is sending me the following warning:

Warning: Matrix is singular to working precision.

What are the reasons for that?

[Walter Roberson]

how does rank of the matrix compare to its size?

[Marcelo Boldt]

The rank of the matrix is 354 < to the matrix's dimension. How do I proceed?

[Walter Roberson]

You have significant problems in how you construct your matrix. 354 is not even close to 443.

• none of your rows may be all 0
• none of your columns may be all zero
• none of your rows may be identical to other rows
• none of your columns may be identical to other columns.

[Marcelo Boldt]

Now I improved the code and I am getting way closer rank = 442 < 443. Thanks for the guidance

[Bruno Luong]

Then the null space is dimension 1.

You might "plot" null(K) and see the form of this "blind" solution. It migh give a hint what condition is still missing (such as integral solution == 0) in your equation system.

[Walter Roberson]

One of the ways this sort of problem can arise is if people accidentally miss one degree of freedom in a solution, such as carefully constraining everything on the XY plane but not noticing that at one location they did not prevent twisting into Z space. The equivalent of forgetting to put a nut on a bolt .

Answer 3

[Bjorn Gustavsson]

Then you proceed by reading up on linear inverse problem theory, you can start with looking at Tikhonov regularization. In physics and engineering the forward/theory matrix, your K, are inherently underdetermined. There are a well devloped set of methods and techniques to properly handle this. There are a good matlab toolbox on the file exchange for handling this: regtools, it also contains an introduction to these types of problems.

HTH

[Marcelo Boldt]

Thank you very much for the guidance!