% gauss_el.m % Programmed by ¼°æ¿ø
function x = gauss_el(A,b); n = length(A); R = [A b];
for j=1:n-1 for i=j+1:n if R(j,j)==0 max_id = find( max( abs(R(j+1:n,j)) ) ) + j; % Find pivot temp = R(j,:); % Change two rows R(j,:) = R(max_id,:); R(max_id,:) = temp; end a = R(i,j)/R(j,j); R(i,:) = R(i,:) - a*R(j,:); end end A = R(:,1:n); b = R(:,n+1); x = zeros(n,1); for j=n:-1:1 x(j) = ( b(j) - A(j,j+1:n)*x(j+1:n) )/A(j,j); end norm( b-A*x )
%%%%
>> A = myrand(8,8,1,10); b = myrand(8,1,10,20)
b =
17 13 14 19 17 12 18 16
>> x = gauss_el(A,b)
ans =
1.9860e-015
x =
1.8490 0.2148 -0.4661 1.3613 -0.0846 -0.0117 0.8911 -1.3097
>>
>> A = myrand(10,10,1,10); b = myrand(10,1,10,20);
>> gauss_el(A,b)
ans =
6.3632e-015
ans =
-7.2008 3.4116 5.7833 -1.6668 2.8751 -0.2326 -2.5189 3.5883 1.6561 -1.0576
>>
|
LAST UPDATE: 2013.10.10 - 15:00 |
|