HomogenousEquationSolver

HomogenousEquationSolver.lua

function find_pivot_positions(matrix)
	local m, n = #matrix, #matrix[1]
	local pivot_positions = {}
	for i = 1, m do
		for j = 1, n do
			if matrix[i][j] == 1 then
				table.insert(pivot_positions, {i, j})
				break
			end
		end
	end
	return pivot_positions
end
-- Function to interchange rows to bring the matrix into upper triangular form
function interchangeRows(matrix)
	local non_zero_indices = {}
	for row_index = 1, #matrix do
		local row_has_non_zero_element = true
		for element_index = 1, #matrix[row_index] do
			if matrix[row_index][element_index] ~= 0 then
				if non_zero_indices[element_index] then
					table.insert(non_zero_indices[element_index], row_index)
				else
					non_zero_indices[element_index] = {row_index}
				end
				row_has_non_zero_element = false
				break
			end
		end
		if row_has_non_zero_element then
			local last_element_index = #matrix[row_index]
			if non_zero_indices[last_element_index] then
				table.insert(non_zero_indices[last_element_index], row_index)
			else
				non_zero_indices[last_element_index] = {row_index}
			end
		end
	end

	local sorted_non_zero_indices = {}
	for key in pairs(non_zero_indices) do
		table.insert(sorted_non_zero_indices, key)
	end
	table.sort(sorted_non_zero_indices)

	local sorted_matrix = {}
	for _, column_index in ipairs(sorted_non_zero_indices) do
		for _, row_index in ipairs(non_zero_indices[column_index]) do
			table.insert(sorted_matrix, matrix[row_index])
		end
	end

	return sorted_matrix
end

-- Function to scale the pivot row to have a leading entry of 1
function scalePivotRow(matrix, current_row, pivot_index)
	local pivot_value = matrix[current_row][pivot_index]

	for element_index = 1, #matrix[current_row] do
		if pivot_value ~= 0 then
			matrix[current_row][element_index] = matrix[current_row][element_index] / pivot_value
		end
	end

	return matrix
end

-- Function to perform row operations to eliminate non-zero entries above and below the pivot
function performRowOperations(matrix, current_row, pivot_index)
	for row_index = 1, #matrix do
		local pivot_multiplier = matrix[row_index][pivot_index]
		if current_row ~= row_index then
			for element_index = 1, #matrix[row_index] do
				matrix[row_index][element_index] = matrix[row_index][element_index] - matrix[current_row][element_index] * pivot_multiplier
				if math.abs(matrix[row_index][element_index]) < 0.00000000000001 then
					matrix[row_index][element_index] = 0
				end
			end
		end
	end

	return matrix
end

-- Function to compute the Reduced Row Echelon Form (RREF) of a matrix
function rref(matrix)
	for current_row = 1, (#matrix-1) do
		matrix = interchangeRows(matrix)
		local pivot_index = 0
		for element_index = 1, #matrix[current_row] do
			if matrix[current_row][element_index] ~= 0 then
				pivot_index = element_index
				break
			end
		end
		matrix = scalePivotRow(matrix, current_row, pivot_index)
		matrix = performRowOperations(matrix, current_row, pivot_index)
	end

	return matrix
end



local function deepCopy(original)
	local copy = {}
	for k, v in pairs(original) do
		if type(v) == "table" then
			v = deepCopy(v)
		end
		copy[k] = v
	end
	return copy
end


function HomogeneousSolver(Matrix)
	local mat = deepCopy(Matrix)
	local RREF = rref(mat)

	local pivots = find_pivot_positions(RREF)
	local n_pivots = #pivots
	local no_of_cols = #RREF[1]
	local n_free = no_of_cols - n_pivots
	if n_free == 0 then
		error("Solution doesn't exist")
	else
		local ans_mat = {}
		for i = 1, n_free + 1 do
			local row = {}
			for j = 1, no_of_cols do
				table.insert(row, 0)
			end
			table.insert(ans_mat, row)
		end

		local non_pivot_cols = {}
		for i = 1, no_of_cols do
			local found = false
			for _, pivot in ipairs(pivots) do
				
				if pivot[2] == i then
					found = true
					break
				end
			end
			if not found then
				table.insert(non_pivot_cols, i)
			end
		end

		for i = 1, #mat do
			local pivot = -1
			for j = 1, no_of_cols do
				if RREF[i][j] == 1 then
					pivot = j
				end
				if pivot ~= -1 and pivot ~= j and RREF[i][j] ~= 0 then
					local ind = table.find(non_pivot_cols, j) + 1
					ans_mat[ind][pivot] = ans_mat[ind][pivot] - RREF[i][j]
				end
			end
		end

		for i = 1, #non_pivot_cols do
			local ind = table.find(non_pivot_cols, non_pivot_cols[i]) + 1
			ans_mat[ind][non_pivot_cols[i]] = 1
		end
		return ans_mat
	end
end

return HomogeneousSolve