Question

A cubic polynomial is of the form: f(x) = Ax^3 + Bx^2 + Cx + D   

A root of the polynomial is a value, x, such that f(x)=0.  

Write a program that takes in the coefficients of a cubic polynomial: A, B, C, and D. The program finds and reports all three roots of the polynomial.

Hint: First use bisection method to determine a single root of the polynomial. Then divide the polynomial by its factor to come up with a quadratic equation. You may use the procedure described here: https://www.purplemath.com/modules/synthdiv.htm. Subsequently, solve the quadratic equation.

The code should be written in python.

Answer 1

NOTE: PLEASE REFER IMAGE FOR INDENTATION

import math

import numpy as np

def solve(a, b, c, d):

    if (a == 0 and b == 0):                     # Case for handling Liner Equation

        return np.array([(-d * 1.0) / c])                 # Returning linear root as numpy array.

    elif (a == 0):                              # Case for handling Quadratic Equations

        D = c * c - 4.0 * b * d                       # Helper Temporary Variable

        if D >= 0:

            D = math.sqrt(D)

            x1 = (-c + D) / (2.0 * b)

            x2 = (-c - D) / (2.0 * b)

        else:

            D = math.sqrt(-D)

            x1 = (-c + D * 1j) / (2.0 * b)

            x2 = (-c - D * 1j) / (2.0 * b)

           

        return np.array([x1, x2])               # Returning Quadratic Roots as numpy array.

    f = findF(a, b, c)                          # Helper Temporary Variable

    g = findG(a, b, c, d)                       # Helper Temporary Variable

    h = findH(g, f)                             # Helper Temporary Variable

    if f == 0 and g == 0 and h == 0:            # All 3 Roots are Real and Equal

        if (d / a) >= 0:

            x = (d / (1.0 * a)) ** (1 / 3.0) * -1

        else:

            x = (-d / (1.0 * a)) ** (1 / 3.0)

        return np.array([x, x, x])              # Returning Equal Roots as numpy array.

    elif h <= 0:                                # All 3 roots are Real

        i = math.sqrt(((g ** 2.0) / 4.0) - h)   # Helper Temporary Variable

        j = i ** (1 / 3.0)                      # Helper Temporary Variable

        k = math.acos(-(g / (2 * i)))           # Helper Temporary Variable

        L = j * -1                              # Helper Temporary Variable

        M = math.cos(k / 3.0)                  # Helper Temporary Variable

        N = math.sqrt(3) * math.sin(k / 3.0)    # Helper Temporary Variable

        P = (b / (3.0 * a)) * -1                # Helper Temporary Variable

        x1 = 2 * j * math.cos(k / 3.0) - (b / (3.0 * a))

        x2 = L * (M + N) + P

        x3 = L * (M - N) + P

        return np.array([x1, x2, x3])           # Returning Real Roots as numpy array.

    elif h > 0:                                 # One Real Root and two Complex Roots

        R = -(g / 2.0) + math.sqrt(h)           # Helper Temporary Variable

        if R >= 0:

            S = R ** (1 / 3.0)                  # Helper Temporary Variable

        else:

            S = (-R) ** (1 / 3.0) * -1          # Helper Temporary Variable

        T = -(g / 2.0) - math.sqrt(h)

        if T >= 0:

            U = (T ** (1 / 3.0))                # Helper Temporary Variable

        else:

            U = ((-T) ** (1 / 3.0)) * -1        # Helper Temporary Variable

        x1 = (S + U) - (b / (3.0 * a))

        x2 = -(S + U) / 2 - (b / (3.0 * a)) + (S - U) * math.sqrt(3) * 0.5j

        x3 = -(S + U) / 2 - (b / (3.0 * a)) - (S - U) * math.sqrt(3) * 0.5j

        return np.array([x1, x2, x3])           # Returning One Real Root and two Complex Roots as numpy array.

# Helper function to return float value of f.

def findF(a, b, c):

    return ((3.0 * c / a) - ((b ** 2.0) / (a ** 2.0))) / 3.0

# Helper function to return float value of g.

def findG(a, b, c, d):

    return (((2.0 * (b ** 3.0)) / (a ** 3.0)) - ((9.0 * b * c) / (a **2.0)) + (27.0 * d / a)) /27.0

# Helper function to return float value of h.

def findH(g, f):

    return ((g ** 2.0) / 4.0 + (f ** 3.0) / 27.0)

x=solve(1,2,3,4)

print(x)