Question:
Solve the following system of linear equations using Gaussian Elimination using Fortran.
x1
+ x2 - x3 + x4 - x5 = 2
2x1+
2 x2 + x3 - x4 + x5 = 4
3x1
+ x2-3 x3 -2x4+ 3 x5 = 8
4x1
+ x2 - x3 +4 x4 -5x5 =16
16x1
-x2 + x3 - x4 - x5 =32
Solution:
program main
implicit none
integer, parameter :: n=5
double precision:: a(n,n), b(n), x(n)
integer:: i,j
! matrix A
data (a(1,i), i=1,5) / 1.0, 1.0, -1.0, 1.0, -1.0 /
data (a(2,i), i=1,5) / 2.0, 2.0, 1.0, -1.0, 1.0 /
data (a(3,i), i=1,5) / 3.0, 1.0, -3.0, -2.0, 3.0 /
data (a(4,i), i=1,5) / 4.0, 1.0, -1.0, 4.0, -5.0 /
data (a(5,i), i=1,5) / 16.0, -1.0, 1.0, -1.0, -1.0 /
! matrix b
data (b(i), i=1,5) / 2.0, 4.0, 8.0, 16.0, 32.0 /
! print a header and the original equations
write (*,200)
do i=1,n
write (*,201) (a(i,j),j=1,n), b(i)
end do
call gauss_2(a,b,x,n)
! print matrix A and vector b after the elimination
write (*,202)
do i = 1,n
write (*,201) (a(i,j),j=1,n), b(i)
end do
! print solutions
write (*,203)
write (*,201) (x(i),i=1,n)
200 format (' Gauss elimination with scaling and pivoting ' &
,/,/,' Matrix A and vector b')
201 format (6f12.6)
202 format (/,' Matrix A and vector b after elimination')
203 format (/,' Solutions x(n)')
end
subroutine gauss_2(a,b,x,n)
implicit none
integer:: n
double precision:: a(n,n), b(n), x(n)
double precision:: s(n)
double precision:: c, pivot, store
integer:: i, j, k, l
! step 1: begin forward elimination
do k=1, n-1
! step 2: "scaling"
! s(i) will have the largest element from row i
do i=k,n
s(i) = 0.0
do j=k,n
s(i) = max(s(i),abs(a(i,j)))
end do
end do
! "pivoting 1"
! row with the largest pivoting element
pivot = abs(a(k,k)/s(k))
l = k
do j=k+1,n
if(abs(a(j,k)/s(j)) > pivot) then
pivot = abs(a(j,k)/s(j))
l = j
end if
end do
! Check if the system has a sigular matrix
if(pivot == 0.0) then
write(*,*) ' The matrix is sigular '
return
end if
! "pivoting 2" interchange rows k and l (if needed)
if (l /= k) then
do j=k,n
store = a(k,j)
a(k,j) = a(l,j)
a(l,j) = store
end do
store = b(k)
b(k) = b(l)
b(l) = store
end if
! elimination (after scaling and pivoting)
do i=k+1,n
c=a(i,k)/a(k,k)
a(i,k) = 0.0
b(i)=b(i)- c*b(k)
do j=k+1,n
a(i,j) = a(i,j)-c*a(k,j)
end do
end do
end do
! backward substitution
x(n) = b(n)/a(n,n)
do i=n-1,1,-1
c=0.0
do j=i+1,n
c= c + a(i,j)*x(j)
end do
x(i) = (b(i)- c)/a(i,i)
end do
end subroutine gauss_2
