一元三次、四次方程求解

一、一元四次方程求解

一元四次方程求根公式，百科：        https://baike.baidu.com/item/%E4%B8%80%E5%85%83%E5%9B%9B%E6%AC%A1%E6%96%B9%E7%A8%8B%E6%B1%82%E6%A0%B9%E5%85%AC%E5%BC%8F/10721996?fr=aladdin

1.网上搜到的求解代码

python求解代码，见：https://github.com/Larissa1990/Solve-cubic-and-quartic-equations-with-one-unknown/blob/master/Equations.py   经测试， 有bug，求解不正确

说明，见：https://www.cnblogs.com/larissa-0464/p/11706131.html

 

2. 一元四次方程在线求解工具：https://www.osgeo.cn/app/s2085   可用来验证求解程序是否正确；

 

3. 一元四次方程， 沈天珩简化求根公式,  具体公式 详见百科链接；

import math
import cmath
import numpy as np


def cal_quartic_ik(args_list):
    a, b, c, d, e = args_list

    D = 3*pow(b,2) - 8*a*c
    E = -pow(b, 3) + 4*a*b*c - 8*pow(a, 2)*d
    F = 3*pow(b, 4) + 16*pow(a, 2)*pow(c, 2) - 16*a*pow(b, 2)*c + 16*pow(a, 2)*b*d - 64*pow(a, 3)*e

    A = D**2 - 3*F
    B = D*F - 9*pow(E, 2)
    C = F**2 - 3*D*pow(E, 2)

    delta = B**2 - 4*A*C  # 总判别式

    if (D == 0) & (E == 0) & (F == 0):
        """ 四重实根"""
        x = -b/(4*a)
        return 1, [x]
    if (A == 0) & (B == 0) & (C == 0) & (D*E*F != 0):
        """ 两个实根，其中一个三重实根"""
        x1 = (-b*D + 9*E)/(4*a*D)
        x234 = (-b * D - 3 * E) / (4 * a * D)
        return 2, [x1, x234]
    if (E == 0) & (F == 0) & (D != 0):
        """ 一对二重根"""
        if D>0:  # 根为实数
            x13 = (-b + math.sqrt(D))/(4*a)
            x24 = (-b - math.sqrt(D)) / (4 * a)
            return 2, [x13, x24]

        if D<0:  # 根为虚数
            # x13 = (-b + cmath.sqrt(D))/(4*a)
            # x24 = (-b - cmath.sqrt(D)) / (4 * a)
            return 0, 0
    if (A*B*C != 0) & (delta == 0):
        """ 一对二重实根 """
        x3 = (-b - np.sign(A*B*E)*math.sqrt( D - B/A))/(4*a)
        x4 = (-b - np.sign(A*B*E)*math.sqrt( D - B/A))/(4*a)
        if A*B>0 :  # 其余两根为不等实根
            x1 = (-b + np.sign(A*B*E)*math.sqrt( D - B/A) + math.sqrt( 2*B/A) )/(4*a)
            x2 = (-b + np.sign(A * B * E) * math.sqrt(D - B / A) - math.sqrt(2 * B / A)) / (4 * a)
            return 4, [x1, x2, x3, x4]
        if A*B < 0:  # 其余两根为共轭虚根
            # x1 = (-b + np.sign(A * B * E) * math.sqrt(D - B / A) + cmath.sqrt(2 * B / A)) / (4 * a)
            # x2 = (-b + np.sign(A * B * E) * math.sqrt(D - B / A) - cmath.sqrt(2 * B / A)) / (4 * a)
            return 2,  [x3, x4]
    if delta > 0:
        """" 两个不等实根和一对共轭虚根"""
        z1 = A*D + 3*(( -B + math.sqrt(delta))/2.0)
        z2 = A * D + 3 * ((-B - math.sqrt(delta)) / 2.0)

        # print """ z1 =  """, z1
        # print """ z2 =  """, z2
        # print """ abs(z1) =  """, abs(z1)
        # print """ abs(z2) =  """, abs(z2)

        z = D**2 - D*(np.sign(z1)*pow(abs(z1), 1.0/3.0) + np.sign(z2)*pow(abs(z2), 1.0/3.0)) + \
            (np.sign(z1)*pow(abs(z1), 1.0/3.0) + np.sign(z2)*pow(abs(z2), 1.0/3.0))**2 - 3*A

        x1 = (-b + np.sign(E)*math.sqrt((D + np.sign(z1)*pow(abs(z1), 1.0/3.0) + np.sign(z2)*pow(abs(z2), 1.0/3.0))/3.0)
              + math.sqrt((2*D - np.sign(z1)*pow(abs(z1), 1.0/3.0) - np.sign(z2)*pow(abs(z2), 1.0/3.0)
                 + 2*math.sqrt(z))/3.0))/(4*a)
        x2 = (-b + np.sign(E)*math.sqrt((D + np.sign(z1)*pow(abs(z1), 1.0/3.0) + np.sign(z2)*pow(abs(z2), 1.0/3.0))/3.0)
              - math.sqrt((2*D - np.sign(z1)*pow(abs(z1), 1.0/3.0) - np.sign(z2)*pow(abs(z2), 1.0/3.0)
                 + 2*math.sqrt(z))/3.0))/(4*a)

        # 虚根忽略
        return 2, [x1, x2]
    if delta < 0:
        if E == 0:
            if (D>0) & (F>0) :
                """ 四个不等实根 """
                x1 = (-b + math.sqrt(D + 2*math.sqrt(F))) / (4 * a)
                x2 = (-b - math.sqrt(D + 2 * math.sqrt(F))) / (4 * a)
                x3 = (-b + math.sqrt(D - 2 * math.sqrt(F))) / (4 * a)
                x4 = (-b - math.sqrt(D - 2 * math.sqrt(F))) / (4 * a)
                return 4, [x1, x2, x3, x4]
            else:
                """ 两对不等共轭虚根 """
                # 虚根忽略
                print " 两对不等共轭虚根 "
                return 0, 0
        else:
            if (D > 0) & (F > 0):
                """ 四个不等实根 """
                theta = math.acos((3*B-2*A*D)/(2*A*math.sqrt(A)))
                y1 = (D - 2*math.sqrt(A)*math.cos(theta/3.0))/3.0
                y2 = (D + math.sqrt(A)*(math.cos(theta/3.0) + math.sqrt(3)*math.sin(theta/3.0)))/3.0
                y3 = (D + math.sqrt(A) * (math.cos(theta / 3.0) - math.sqrt(3) * math.sin(theta / 3.0))) / 3.0

                x1 = (-b + np.sign(E) * math.sqrt(y1) + ( math.sqrt(y2) + math.sqrt(y3)))/(4*a)
                x2 = (-b + np.sign(E) * math.sqrt(y1) - (math.sqrt(y2) + math.sqrt(y3))) / (4 * a)
                x3 = (-b - np.sign(E) * math.sqrt(y1) + (math.sqrt(y2) - math.sqrt(y3))) / (4 * a)
                x4 = (-b - np.sign(E) * math.sqrt(y1) - (math.sqrt(y2) - math.sqrt(y3))) / (4 * a)

                return 4, [x1, x2, x3, x4]
            else:
                """ 两对不等共轭虚根 """
                # 虚根忽略
                print " 两对不等共轭虚根 "
                return 0, 0

        经测试，求解正确；

         程序中，忽略了复数解，返回的解，均是实数解； 可作为一基础求解工具；

 

一、一元三次方程求解

 

盛金公式： https://baike.baidu.com/item/%E7%9B%9B%E9%87%91%E5%85%AC%E5%BC%8F#9

def cal_cubic_ik(args_list):
    a, b, c, d = args_list

    A = b**2 - 3*a*c
    B = b*c - 9*a*d
    C = c**2 - 3*b*d

    delta = B**2 - 4*A*C  # 总判别式

    if (A == 0) & (B == 0):
        """ 一个三重实根"""
        x = -b/(3*a)
        return 1, [x]

    if delta > 0:
        """" 一个实根和一对共轭复根"""

        y1 = A*b + 3*a*(( -B + math.sqrt(delta))/2.0)
        y2 = A*b + 3*a*(( -B - math.sqrt(delta))/2.0)
        print """ y1 =  """, y1
        print """ y2 =  """, y2
        print """ abs(y2) =  """, abs(y2)
        # 虚根忽略
        x1 = (-b - (np.sign(y1)*pow(abs(y1), 1.0 / 3.0) + np.sign(y2)*pow(abs(y2), 1.0 / 3.0)) )/ (3 * a)  # 负数直接开立方根会报错
        return 1, [x1]

    if delta == 0:
        """三个实根，其中有一个二重根 """
        K = B/float(A)
        x1 = -b/float(a) + K
        x23 = -K/2.0
        return 2,  [x1, x23]

    if delta < 0:
        """ 三个不等实根 """
        T = (2*A*b - 3*a*B)/(2*math.sqrt(pow(A, 3)))
        theta = math.acos(T)

        x1 = (-b - 2 * math.sqrt(A) * math.cos(theta / 3.0)) / (3 * a)
        x2 = (-b + math.sqrt(A)*(math.cos(theta/3.0)) + math.sqrt(3)*math.sin(theta/3.0))/(3*a)
        x3 = (-b + math.sqrt(A) * (math.cos(theta / 3.0)) - math.sqrt(3) * math.sin(theta / 3.0)) / (3 * a)

        return 3, [x1, x2, x3]

 程序中，忽略了复数解，返回的解，均是实数解；