Post by dj54 on Dec 19, 2022 10:27:16 GMT -5
The Liberty Basic code provides a calling routine that illustrates the use of the system solver in the included library. It also prints the original system of equations along with the solution, error in the solution, and run time. You can specify the size of the matrices involved. It has been tested to solve a system of 1000 linear equations (not recommended due to time required).
' This is the calling program to illustrate the
' use of the subroutine LSSGaussianElim
'
' The printed output produces this sample output, headings auto adjust for matrix size "n":
'
'Gaussian Elimination Linear System Solver
'-----------------------------------------
'Matrix size n: 10
'Elapsed time in milliseconds: 7
'Error message: No singularities encountered
'
'Linear System AX = b where X is the solution vector
'
'---------------------------- A ---------------------------- --- X --- - b -
' 9.0 69.0 21.0 1.0 79.0 22.0 56.0 36.0 55.0 3.0 1.29860826 | 60.0
' 81.0 69.0 29.0 92.0 69.0 44.0 4.0 98.0 99.0 26.0 0.63962691 | 100.0
' 38.0 68.0 61.0 89.0 28.0 63.0 42.0 30.0 89.0 11.0 -2.04725213 | 81.0
' 75.0 34.0 8.0 25.0 14.0 10.0 82.0 53.0 61.0 70.0 1.01993463 | 92.0
' 62.0 45.0 85.0 51.0 71.0 73.0 53.0 42.0 38.0 39.0 0.04491942 | 63.0
' 87.0 74.0 91.0 29.0 14.0 86.0 10.0 89.0 68.0 7.0 1.27646213 | 44.0
' 35.0 8.0 31.0 72.0 47.0 98.0 52.0 69.0 91.0 88.0 1.47122975 | 33.0
' 46.0 18.0 6.0 93.0 10.0 46.0 24.0 31.0 97.0 55.0 0.08882901 | 54.0
' 77.0 52.0 66.0 86.0 86.0 24.0 30.0 53.0 84.0 39.0 -1.21343471 | 11.0
' 3.0 23.0 2.0 60.0 14.0 62.0 25.0 87.0 60.0 23.0 -1.44265765 | 94.0
'
'Error in solution vector X:
'---------------------------
'Maximum error = 0
'Minimum error = -0.20605739e-12
'Avg error = 0.0 where Avg error=SQR(Sum(Error^2))
'Below you should replace n=10 for the number of linear equations to be solved.
n=10
dim A(n,n+1), X(n), AA(n,n+1), XError(n)
global n, A, X, ErrorMsg$
'******************************************************************
' Assign value to matrix A and save original values in matrix AA
'******************************************************************
for i= 1 to n
for j = 1 to n+1
A(i,j) = int(rnd(1)*100) + 1 : REM 1 <= Aij <= 100
AA(i,j) = A(i,j)
next j
next i
StartTime = time$("ms") : REM milliseconds since midnight
call LSSGaussianElim
StopTime = time$("ms") : REM milliseconds since midnight
'******************************************************
' Print the Linear System and the Solution Vector X(k)
'******************************************************
print "Gaussian Elimination Linear System Solver"
print "-----------------------------------------"
print "Matrix size n: ";n
print "Elapsed time in milliseconds: "; StopTime-StartTime
print "Error message: ";ErrorMsg$
print
print "Linear System AX = b where X is the solution vector"
print
heading$ = RptStr$("-",3*n-2) + " A " + RptStr$("-",3*n-2) + " --- X ---" + " - b -"
print heading$
'****************************
' Print Solution Vector X(k)
'****************************
template$ = "###.#"
template2$ = "####.##############"
for i = 1 to n
for j = 1 to n
print using(template$,AA(i,j)); " ";
next j
print using(template2$,X(i));" | "; using(template$,AA(i,n+1))
next i
print
'***********************************************************************
' Using the original system, calculate the error in solution vector X()
'***********************************************************************
XErrorMax = 0
XErrorMin = 0
XErrorSum = 0
for i = 1 to n
XError(i) = 0
for j= 1 to n
XError(i) = XError(i) + AA(i,j)*X(j)
next j
XError(i) = AA(i,n+1) - XError(i)
if XError(i) < XErrorMin then XErrorMin = XError(i)
if XError(i) > XErrorMax then XErrorMax = XError(i)
XErrorSum = XErrorSum + XError(i)^2
next i
print "Error in solution vector X:"
print "---------------------------"
print "Maximum error = ";XErrorMax
print "Minimum error = ";XErrorMin
print "Avg error = ";sqr(XErrorSum);" where Avg error=SQR(Sum(Error^2))"
end
'*******************************************************************************************
'*******************************************************************************************
'** Library of Functions and Subroutines **
'*******************************************************************************************
'*******************************************************************************************
function round(nnum,n)
' nnum=number to be rounded, n=number of decimal places
round = int(100*nnum)/(10^n)
end function
'------------------------------------------------------------------------------------------
function RptStr$(char$,length)
'Builds a string of length "length" using "char$"
'char$=character(s) to be repeated
'length=number of times char$ is to be repeated
RptStr$ = ""
for i = 1 to length
RptStr$ = RptStr$ + char$
next i
end function
'------------------------------------------------------------------------------------------
function StrFormat$(string$,strlen,just$,trail)
' string$ - string to be formatted
' strlen - output length of string which should be the greater of column title or longest anticipated string
' just$ - justification L or R
' trail - number of blanks to append to the end of the string$ - value of 1 means 1 space between columns
if len(trim$(string$)) < strlen then
blanks = strlen - len(trim$(string$))
if just$ = "L" then
StrFormat$ = trim$(string$) + space$(blanks) + space$(trail)
else
StrFormat$ = space$(blanks) + trim$(string$) + space$(trail)
end if
else
StrFormat$ = left$(trim$(string$),strlen) + space$(trail)
end if
end function
'------------------------------------------------------------------------------------------
sub LSSGaussianElim
'Gaussian Elimination Algorithm for Linear Systems
'
'USAGE: in your calling program, include these 3 statements:
' n=10 change to 10 to the size of your A matrix
' dim A(n,n+1), X(n)
' global n, A, X, ErrorMsg$
'
' Solve the linear system AX = b where
'
' A is the n x n matrix of coefficients to the unknowns X(1), ..., X(n) and
' b is the vector of constants b(1), ..., b(n)n
'
' In the algorithm below, you must augment the matrix A with column n+1 to house
' the values of the vector b. Thus b(i) = A(i,n+1) for i = 1, ..., n
'
'1. Start
'2. Declare the variables and read the order of the matrix n.
'3. Take the coefficients of the linear equation as:
' Do for k=1 to n
' Do for j=1 to n+1
' Read Akj
' End for j
' End for k
'
'Do for k=1 to n-1
' Do for i=k+1 to n
' Do for j=k+1 to n+1
' Aij = Aij - (Aik / Akk) * Akj
' End for j
' End for i
'End for k
'Xn = Ann+1 / Ann
'Do for k=n-1 to 1 step -1
' sum = 0
' Do for j=k+1 to n
' sum = sum + Akj * Xj
' End for j
' Xk = (1/Akk) * (Akn+1 - sum)
'End for k
'
'7. Display the result Xk
'8. Stop
'
'ERROR HANDLING (in the upper triangular matrix calculation)
'A check is performed to avoid dividing by zero. This is called a singularity.
'If a singularity occurs along the diagonal (when Ajj=0), then ErrorMsg$
'is set to "STOP: Matrix is singular, upper triangular matrix k,i,j=" and
'values supplied for indexes k,i,j. If you see this message, the rest of the
'printed solution vector X() is NOT correct. The zero value for Ajj is
'artificially set to 1. If there are NO zero values, thus no singularities,
'then the ErrorMsg$ is set to "No singularities encountered".
'
'Runtimes to Solve Matrix
'on HP Pavilion 360 made in 2019
'
'N= Milliseconds
'--- ------------
'5 1
'10 3
'20 11
'35 78
'50 250
'75 878
'90 1412
'100 1884 = 1.884 seconds
'200 14243 = 14.243 seconds
'400 111795 = 1.86 minutes
'600 380505 = 6.34 minutes
'
'Using Excel to draw a scatter plot and curve fit the growth in
'milliseconds as function of N yields this nearly perfect fit
'Milliseconds = 0.0018N^3 - 0.0435N^2 + 6.9428N - 128.01
'
'As you get past N=100, the growth in solve time may suggest
'using a faster algorithm, depending on the time sensitivity
'of the real-world application. If your application requires
'solving multi-thousands of N=1000+ systems, such as might
'be encountered with partial differential equations, then
'definitely consider a faster algorithm.
'*****************************************
' Calculate the upper triangular matrix
'*****************************************
for k=1 to n-1
for i=k+1 to n
for j=k+1 to n+1
if A(k,k) = 0 then
ErrorMsg$ = "STOP: Matrix is singular, upper triangular matrix k,i,j= "+str$(k)+" "+str$(i)+" "+str$(j)
EXIT FOR
else
A(i,j) = A(i,j) - (A(i,k) / A(k,k)) * A(k,j)
ErrorMsg$ = "No singularities encountered"
end if
next j
if left$(ErrorMsg$,4) = "STOP" then EXIT FOR
next i
if left$(ErrorMsg$,4) = "STOP" then EXIT FOR
next k
'*****************************************
' Backsolve to calculate values for X()
'*****************************************
X(n) = A(n,n+1) / A(n,n)
for k=n-1 to 1 step -1
sum = 0
for j=k+1 to n
sum = sum + A(k,j)*X(j)
next j
X(k) = ( A(k,n+1) - sum ) / A(k,k)
next k
end sub
'------------------------------------------------------------------------------------------