CODEWITHSUNDEEP
pythonnumpymathsympy
rambler
rambler

Reputation: 25

python solve Cubic equation without using sympy

Is it possible to solve Cubic equation without using sympy?

Example:

import sympy as sp
xp = 30
num = xp + 4.44
sp.var('x, a, b, c, d')
Sol3 = sp.solve(0.0509 * x ** 3 + 0.0192 * x ** 2 + 3.68 * x - num, x)

The result is:

[6.07118098358257, -3.2241955998463 - 10.0524891203436*I, -3.2241955998463 + 10.0524891203436*I]

But I want to find a way to do it with numpy or without 3 part lib at all

I tried with numpy:

import numpy as np
coeff = [0.0509, 0.0192, 3.68, --4.44]
print(np.roots(coeff))

But the result is :

[ 0.40668245+8.54994773j  0.40668245-8.54994773j -1.19057511+0.j]

Upvotes: 2

Views: 2580

Answers (2)

Alain T.
Alain T.

Reputation: 42143

You could implement the cubic formula

enter image description here

this Youtube video from mathologer could help understand it.

Based on that, the cubic function for ax^3 + bx^2 + cx + d = 0 can be written like this:

def cubic(a,b,c,d):
    n = -b**3/27/a**3 + b*c/6/a**2 - d/2/a
    s = (n**2 + (c/3/a - b**2/9/a**2)**3)**0.5
    r0 = (n-s)**(1/3)+(n+s)**(1/3) - b/3/a
    r1 = (n+s)**(1/3)+(n+s)**(1/3) - b/3/a
    r2 = (n-s)**(1/3)+(n-s)**(1/3) - b/3/a
    return (r0,r1,r2)

The simplified version of the formula only needs to get c and d as parameters (aka p and q) and can be implemented like this:

def cubic(p,q):
    n = -q/2
    s = (q*q/4+p**3/27)**0.5
    r0 = (n-s)**(1/3)+(n+s)**(1/3)
    r1 = (n+s)**(1/3)+(n+s)**(1/3)
    r2 = (n-s)**(1/3)+(n-s)**(1/3)
    return (r0,r1,r2)

print(cubic(-15,-126))
(5.999999999999999, 9.999999999999998, 2.0)

I'll let you mix in complex number operations to properly get all 3 roots

Upvotes: 1

scientist
scientist

Reputation: 147

In your numpy method you are making two slight mistakes with the final coefficient.

In the SymPy example your last coefficient is - num, this is, according to your code: -num = - (xp + 4.44) = -(30 + 4.44) = -34.44

In your NumPy example yout last coefficient is --4.44, which is 4.44 and does not equal -34.33.

If you edit the NumPy code you will get:

import numpy as np
coeff = [0.0509, 0.0192, 3.68, -34.44]
print(np.roots(coeff))
[-3.2241956 +10.05248912j -3.2241956 -10.05248912j
  6.07118098 +0.j        ]

The answer are thus the same (note that NumPy uses j to indicate a complex number. SymPy used I)

Upvotes: 4

Related Questions