How to Solve this Big System of Equations

2 次查看（过去 30 天）
Marcelo Boldt
Marcelo Boldt 2021-1-25
评论： Walter Roberson 2021-1-26
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,:) = [];

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

采纳的回答

Bruno Luong
Bruno Luong 2021-1-26
编辑：Bruno Luong 2021-1-26
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.

更多回答（2 个）

Anna Lischke
Anna Lischke 2021-1-26
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
  7 个评论
显示 5更早的评论
Bruno Luong
Bruno Luong 2021-1-26
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
Walter Roberson 2021-1-26
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 .

请先登录，再进行评论。

Bjorn Gustavsson
Bjorn Gustavsson 2021-1-26
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

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by