MATLAB Answers

0

I need help with my code using gaussian elimination.

S 님이 질문을 제출함. 23 Oct 2019
최근 활동 John D'Errico 님이 편집함. 23 Oct 2019
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

  댓글 수: 0

로그인 to comment.

답변 수: 1

Answer by John D'Errico
on 23 Oct 2019
Edited by John D'Errico
on 23 Oct 2019

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.

  댓글 수: 0

로그인 to comment.



Translated by