Solving Cubic Equations With Python
Introduction
In algebra the student learns to find the roots of cubic equations using the Rational Root Theorem. Presented here is python code implementing each of the steps the student is taught. What is NOT done is to go directly to the numpy.roots() function to find the roots since this does not teach the underlying algebraic operations of the Rational Root Theorem.
Algebraic Operations
Definition of a cubic polynomial:
One finds the first root by testing all factors, plus and minus, of the ratio until one solution is found:
Then x minus this root in divided into f(x) using long division, which results in a quadratic equation.
The remaining two roots are found by use of the formula for the roots of a quadratic equation.
The Complete Code
This version works with deprecated python2 as well as with python3 and importing numpy and cmath libraries in needed. The standard math library does not handle complex numbers but cmath does.
#tested with python 3.5.2 and 2.7.12
import numpy as np
import cmath
import sys
def solve_quadratic(a, b, c):
# Check if 'a' is zero, as it would not be a quadratic equation
if a == 0:
raise ValueError("Coefficient 'a' cannot be zero for a quadratic equation.")
# Calculate the discriminant
discriminant = (b**2) - (4 * a * c)
# Calculate the two roots using the quadratic formula
root1 = (-b - cmath.sqrt(discriminant)) / (2 * a)
root2 = (-b + cmath.sqrt(discriminant)) / (2 * a)
return root1, root2
def get_factors(n):
"""
Returns a list of all factors of a given integer n.
"""
factors = []
for i in range(1, n + 1):
if n % i == 0:
factors.append(i)
return factors
def get_factor_ratios(num1, num2):
"""
Calculates all possible ratios of factors of num1 to factors of num2.
"""
factors1 = get_factors(abs(num1))
factors2 = get_factors(abs(num2))
ratios = set() # Using a set to automatically handle unique ratios
for f1 in factors1:
for f2 in factors2:
# Avoid division by zero
if f2 != 0:
ratios.add(f1 / f2)
return sorted(list(ratios)) # Convert to list and sort for ordered output
def cubic_polynomial(x, a, b, c, d):
"""Evaluates a cubic polynomial for a given x and coefficients."""
return a * (x ** 3) + b * (x ** 2) + c * x + d
#get the coefficents as arguments
num = len(sys.argv)
if num != 5:
print ("must provide cubic coefficients as arguments, e.g. A0 A1 A2 A3")
print (" where f(x) = A3x^3 + A2x^2 + A1x + A0")
sys.exit()
A0 = int(sys.argv[1])
A1 = int(sys.argv[2])
A2 = int(sys.argv[3])
A3 = int(sys.argv[4])
if A0 != 0 :
#develop all ratios of A0/A3
all_ratios = get_factor_ratios(A0, A3)
print("All ratios of factors of " + str(A0) + " to factors of " + str(A3) + ":")
print(all_ratios)
#interate over all the ratios, plus and minus, to find
#the first instance of f=0, thus finding a root
for i in range(0, len(all_ratios)):
x = all_ratios[i]
f = cubic_polynomial(x, A3, A2, A1, A0)
if f == 0:
break
if i+1 == len(all_ratios):
#do all negative ratios
for i in range(0, len(all_ratios)):
x = -all_ratios[i]
f = cubic_polynomial(x, A3, A2, A1, A0)
if f == 0:
break
else:
print("continue")
print("function returns: " + str(f) + " for ratio: " + str(x) )
print("This is then a root of the cubic polynomial")
R = x
# Define dividend and divisor polynomials (coefficients from highest to lowest degree)
dividend_coeffs = np.array([A3, A2, A1, A0]) # Represents x^3 - 7x + 6
divisor_coeffs = np.array([1, -R]) # Represents x - 1
# Perform polynomial division
quotient_coeffs, remainder_coeffs = np.polydiv(dividend_coeffs, divisor_coeffs)
print("Quotient coefficients:", quotient_coeffs)
print("Remainder coefficients:", remainder_coeffs)
a = quotient_coeffs[0]
b = quotient_coeffs[1]
c = quotient_coeffs[2]
else:
a = A3
b = A2
c = A1
x = 0
try:
roots = solve_quadratic(a, b, c)
# Print the roots, handling the case of purely real roots
if roots[0].imag == 0 and roots[1].imag == 0:
# print(f"The roots are: {roots[0].real} and {roots[1].real} and {x}")
print("The roots are: " + str(roots[0].real) + " and " + str(roots[1].real) + " and " + str(x))
else:
print("The roots are: " + str(roots[0]) + " and " + str(roots[1]) + " and " + str(x))
except ValueError as e:
print("Error: " + str(e))
except Exception as e:
print("An unexpected error occurred: " + str(e))
"""
examples:
A0 A1 A2 A3
6 5 -12 4
24 2 1 1
-5 9 -5 1
1 0 -7 6
-26 37 -20 4
"""
Code Description
The first part of the code are functions that implement the steps. The first solves the quadratic equation, which is the final step in the implementation of the Rational Root Theorem.
def solve_quadratic(a, b, c):
The next and second function generates all of the integer factors for each of the two coefficients,
def get_factors(n):
and the third returns a sorted list of the ratios.
def get_factor_ratios(num1, num2):
The last function evaluates the polynomial numerically so that a ratio can be tested if it is a root.
def cubic_polynomial(x, a, b, c, d):
The main code body begins with receiving the coefficents a0 through a3 as arguments and converting them from strings to integers. Then the value of a0 (=A0) is tested. If it is non-zero, then the general steps of solution are followed. If it is zero, then one root is x=0 and only the resulting quadratic equation needs to be solved.
#develop all ratios of A0/A3
A sorted list of all integer factors of the ratio is obtained.
#interate over all the ratios, plus and minus, to find
#the first instance of f=0, thus finding a root
Each ratio is tested in the cubic equation to determine if it is a root. First all positive values followed by all negative values. The time a root is found the code progresses to the next step.
# Define dividend and divisor polynomials (coefficients from highest to lowest degree)
Long division is done using the numpy library. To achieve this the code creates two numpy arrays. The first holds the coefficients of the cubic polynomial, and the second holds the root that had just been found.
# Perform polynomial division
The numpy polydiv function is used and it returns a list holding the coefficients of the resulting quadratic equation.
Finally, the solve_quadratic() function is invoked using the coefficients. It is here where the cmath square root function is used with which both real roots and complex number roots can be returned.
Testing
At the end of the code are some sets of cubic coefficients that can be used to test the code. Here is an example.
nvidia@tegra-ubuntu:~/polynomials$ python x_full_cubic_solve.py -26 37 -20 4
Output:
All ratios of factors of -26 to factors of 4:
[0, 1, 2, 3, 6, 13, 26]
continue
function returns: 0 for ratio: 2
This is then a root of the cubic polynomial
('Quotient coefficients:', array([ 4., -12., 13.]))
('Remainder coefficients:', array([ 0.]))
The roots are: (1.5-1j) and (1.5+1j) and 2In this example one can see that there is one real root and two complex roots.
Recent Comments