Gaussian elimination
Problem: Solve Ax=b using Gaussian elimination then backwards substitution. A being an n by n matrix. Also, x and b are n by 1 vectors. To improve accuracy, please use partial pivoting and scaling.
J
%. , J's matrix divide verb, directly solves systems of determined and of over-determined linear equations directly. This example J session builds a noisy sine curve on the half circle, fits quintic and quadratic equations, and displays the results of evaluating these polynomials.
<lang J>
f=: 6j2&": NB. formatting verb
sin=: 1&o. NB. verb to evaluate circle function 1, the sine
add_noise=: ] + -:@:[ - (* ?@:(#&0@:#)) NB. AMPLITUDE add_noise SIGNAL
f RADIANS=: o.@:(%~ i.@:>:)5 NB. monadic circle function is pi times 0.00 0.63 1.26 1.88 2.51 3.14
f SINES=: sin RADIANS 0.00 0.59 0.95 0.95 0.59 0.00
f NOISY_SINES=: 0.1 add_noise SINES _0.01 0.61 0.91 0.99 0.60 0.02
A=: (^/ i.@:#) RADIANS NB. A is the quintic coefficient matrix
NB. display the equation to solve (f A) ; 'x' ; '=' ; f@:,. NOISY_SINES
┌────────────────────────────────────┬─┬─┬──────┐ │ 1.00 0.00 0.00 0.00 0.00 0.00│x│=│ _0.01│ │ 1.00 0.63 0.39 0.25 0.16 0.10│ │ │ 0.61│ │ 1.00 1.26 1.58 1.98 2.49 3.13│ │ │ 0.91│ │ 1.00 1.88 3.55 6.70 12.62 23.80│ │ │ 0.99│ │ 1.00 2.51 6.32 15.88 39.90100.28│ │ │ 0.60│ │ 1.00 3.14 9.87 31.01 97.41306.02│ │ │ 0.02│ └────────────────────────────────────┴─┴─┴──────┘
f QUINTIC_COEFFICIENTS=: NOISY_SINES %. A NB. %. solves the linear system _0.01 1.71 _1.88 1.48 _0.58 0.08
quintic=: QUINTIC_COEFFICIENTS&p. NB. verb to evaluate the polynomial
NB. %. also solves the least squares fit for overdetermined system quadratic=: (NOISY_SINES %. (^/ i.@:3:) RADIANS)&p. NB. verb to evaluate quadratic. quadratic
_0.0200630695393961729 1.26066877804926536 _0.398275112136019516&p.
NB. The quintic is agrees with the noisy data, as it should f@:(NOISY_SINES ,. sin ,. quadratic ,. quintic) RADIANS _0.01 0.00 _0.02 _0.01 0.61 0.59 0.61 0.61 0.91 0.95 0.94 0.91 0.99 0.95 0.94 0.99 0.60 0.59 0.63 0.60 0.02 0.00 0.01 0.02
f MID_POINTS=: (+ -:@:(-/@:(2&{.)))RADIANS
_0.31 0.31 0.94 1.57 2.20 2.83
f@:(sin ,. quadratic ,. quintic) MID_POINTS _0.31 _0.46 _0.79 0.31 0.34 0.38 0.81 0.81 0.77 1.00 0.98 1.00 0.81 0.83 0.86 0.31 0.36 0.27
</lang>
Mathematica
<lang Mathematica>GaussianElimination[A_?MatrixQ, b_?VectorQ] := Last /@ RowReduce[Flatten /@ Transpose[{A, b}]]</lang>
MATLAB
<lang MATLAB> function [ x ] = GaussElim( A, b)
% Ensures A is n by n sz = size(A); if sz(1)~=sz(2)
fprintf('A is not n by n\n');
clear x;
return;
end
n = sz(1);
% Ensures b is n x 1. if n~=sz(1)
fprintf('b is not 1 by n.\n');
return
end
x = zeros(n,1); aug = [A b]; tempmatrix = aug;
for i=2:sz(1)
% Find maximum of row and divide by the maximum
tempmatrix(1,:) = tempmatrix(1,:)/max(tempmatrix(1,:));
% Finds the maximum in column
temp = find(abs(tempmatrix) - max(abs(tempmatrix(:,1))));
if length(temp)>2
for j=1:length(temp)-1
if j~=temp(j)
maxi = j; %maxi = column number of maximum
break;
end
end
else % length(temp)==2
maxi=1;
end
% Row swap if maxi is not 1
if maxi~=1
temp = tempmatrix(maxi,:);
tempmatrix(maxi,:) = tempmatrix(1,:);
tempmatrix(1,:) = temp;
end
% Row reducing
for j=2:length(tempmatrix)-1
tempmatrix(j,:) = tempmatrix(j,:)-tempmatrix(j,1)/tempmatrix(1,1)*tempmatrix(1,:);
if tempmatrix(j,j)==0 || isnan(tempmatrix(j,j)) || abs(tempmatrix(j,j))==Inf
fprintf('Error: Matrix is singular.\n');
clear x;
return
end
end
aug(i-1:end,i-1:end) = tempmatrix;
% Decrease matrix size
tempmatrix = tempmatrix(2:end,2:end);
end
% Backwards Substitution x(end) = aug(end,end)/aug(end,end-1); for i=n-1:-1:1
x(i) = (aug(i,end)-dot(aug(i,1:end-1),x))/aug(i,i);
end
end </lang>