I need help with my code using gaussian elimination.

조회 수: 22 (최근 30일)
S
S 2019년 10월 23일
편집: John D'Errico 2019년 10월 23일
I keep getting NaN as my answer when I use gaussian elimination but according to f = A\b, there's only a couple of them that are NaN and the rest are actual answers. I would like to know how to fix my code to get it using the elimination method. This is my code:
% load
F1 = 10; %tons
F2 = 10;
F3 = 15;
F4 = 10;
F5 = 10;
c = cosd(45); % 45 degree angle
s = sind(45);
% coefficient matrix
A = [ 1 -c 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 -s 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 0 1 c -c 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 s s 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 -c 0 0 0 0 1 -1 c 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 -s 0 0 0 -1 0 0 -s 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 c 0 0 0 0 0 0 -c 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
0 s 0 0 0 -1 0 0 -s 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 1 0 0 0 -c 0 -1 c 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 s 0 0 s 1 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 c 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -s 0 1;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -c 0 0 0 -c 1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -s 0 0 -1 c 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 -1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0];
% load vector
b = [ 0
0
0
F1
0
F2
0
F3
0
0
0
0
0
0
0
F4
0
F5
0
0
0
0
0
0 ];
% solve linear system
% f = A\b
% plotTruss(f);
% function x = GuassNaive(A,b)
[m,n]=size(A);
if m~=n, error('Matrix A must be square'); end
nb = n+1;
Aug = [A,b];
for k = 1:n-1
% [big,i]=max(abs(Aug(k:n,k)));
% ipr=i+k-1;
% if ipr~=k
% Aug([k,ipr],:)=Aug([ipr,k],:);
% end
for i = k+1:n
factor = Aug(i,k)/Aug(k,k);
Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);
end
end
x = zeros(n,1);
x(n)=Aug(n,nb)/Aug(n,n);
for i = n-1:-1:1
x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
end
% end
x

답변 (1개)

John D'Errico
John D'Errico 2019년 10월 23일
편집: John D'Errico 2019년 10월 23일
Well, we COULD try to help you solve it using Gaussian elimination, but actually we can't. In fact, nobody can.
rank(A)
ans =
23
size(A)
ans =
24 24
Your matrix is a singular matrix. Those NaNs you got were an indication that the matrix is singular, but the test above is proof positive. The rank of A is less than the smaller of the row and column dimensions of A.
You cannot use Gaussian elimination to solve the problem. In fact, no exact solution exists at all.
rank([A,b])
ans =
24
This additional test verifies that the vector b CANNOT be written in the form of a linear combination of the columns of A. Therefore, NO solution exits. At best, you could find an approximate solution.
xhat = pinv(A)*b;
xhat
xhat =
-10.644
9.1002
6.4348
-17.079
-17.079
9.9871
1.9418e-15
19.147
-5.0054
9.8118e-15
14.956
1.3986e-14
10.552
9.9924
6.5194
-9.9612
-0.018303
-14.595
15.71
0.06249
10.075
-0.23924
-10.52
-0.16917
[b,A*xhat]
ans =
0 7.1054e-15
0 -7.1054e-15
0 -3.5527e-15
10 9.9871
0 2.9619e-15
10 10
0 -7.8701e-15
15 14.956
0 -0.031245
0 0.031245
0 0.031245
0 9.8118e-15
0 0.012942
0 -0.012942
0 0.044187
10 10
0 0.044187
10 10.075
0 -0.10668
0 1.7764e-14
0 -0.031245
0 0.075433
0 -0.031245
0 0.06249
So you can see the solution is not terrible. The first column is your target vector, b. The second column is our approximate solution.
You can never do any better than that, for this matrix A and right hand side vector b.
It may be possible that your matrix A was created incorrectly, but I can never know how or where your error came from, because only you know what A means and how it was derived.
If I had to make a wild guess as to the cause, I would need to first guess what you were trying to solve. The angles there, as well as the comments about tons of load for the F values, suggest this is a truss problem of some sort. In turn, that suggests you have a problem when you defined the truss. Most of the time, that means you forgot to specify sufficient constraints on one of the nodes of the truss, so the truss is effectively free to translate in space without bound. Or, possibly, you forgot to prevent free rotational motion.

카테고리

Help CenterFile Exchange에서 Structural Analysis에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by