1500字范文,内容丰富有趣,写作好帮手!
1500字范文 > 【math】利用Cardano方法对一元三次方程求解及python实现

【math】利用Cardano方法对一元三次方程求解及python实现

时间:2023-03-18 04:40:04

相关推荐

【math】利用Cardano方法对一元三次方程求解及python实现

文章目录

【参考】【问题描述】求解一元三次方程 【代码实现】现成的包 cardano_method根据公式编写求解代码 【总结】

【参考】

用Cardano方法求解三次方程介绍cardano方法求解下载cardano方法包x^3+1=0求解问题、三次方程反函数问题Micorsoft-Math-solver 微软数学工具WolframAlpha: inverse of a function/反函数百度百科-一元三次方程求根公式

【问题描述】

求解一元三次方程

求解一元三次方程:

a x 3 + b x 2 + c x + d = 0 ax^3+bx^2+cx+d=0 ax3+bx2+cx+d=0

求解有多种方法,其中Cardano方法得到的解为:

一元三次多项式: y = f ( x ) = a x 3 + b x 2 + c x + d y=f(x)=ax^3+bx^2+cx+d y=f(x)=ax3+bx2+cx+d => a x 3 + b x 2 + c x + d − y = 0 ax^3+bx^2+cx+d-y=0 ax3+bx2+cx+d−y=0

【代码实现】

现成的包 cardano_method

可以使用python包:安装cardano_method【参考这里】

pip install cardano_method

该包的使用:CubicEquation函数对应第一个参数列表是:[a, b, c, d],即求解方程的各系数。

from cardano_method.cubic import CubicEquationa = CubicEquation([1, 3, 4, 4])print(a.answers) # j表示虚部后缀# [(-2+0j), (-0.5+1.322875j), (-0.5-1.322875j)]print(a.answers[0].real) # 获取第一个解的实部print(a.answers[0].imag) # 获取第一个解的虚部

但是发现一个问题:解方程 x 3 + 1 = 0 x^3+1=0 x3+1=0时,居然报错分母是0:ZeroDivisionError

a = CubicEquation([1, 0, 0, 1])

但实际上, x 3 + 1 = 0 x^3+1=0 x3+1=0的解是对应: x 1 = − 1 x_1=-1 x1​=−1, x 2 = 1 2 + 3 2 i x_2=\frac{1}{2}+\frac{\sqrt{3}}{2}i x2​=21​+23 ​​i, x 3 = 1 2 − 3 2 i x_3=\frac{1}{2}-\frac{\sqrt{3}}{2}i x3​=21​−23 ​​i(这里 i i i表示虚数)。需要详细看下包中bug如何解决?

根据公式编写求解代码

可以根据对应解的公式写出求解函数:【注意:关于浮点型计算有问题?】

def cardano_solution_v0(a, b, c, d):ab = -b/float(3*a)q = (3*a*c-(b**2)) / (9*(a**2))r = (9*a*b*c-27*(a**2)*d-2*(b**3)) / (54*(a**3))delta_sqrt = (q**3+r**2)**(1.0/2)s = (r+delta_sqrt)**(1.0/3)t = (r-delta_sqrt)**(1.0/3)imag = complex(0, (s-t)*(3**(1.0/3))/2)x1 = s+t+abx2 = -(s+t)/2+ab+imagx3 = -(s+t)/2+ab-imagreturn x1, x2, x3

测试修改保留浮点,还是有问题,跟cardano_method包的结果对不上???【注:cardano_solution_v0cardano_solution_v1都有问题,正确的解见下方函数:cardano_solution

def round_ri(xo, n=4):xr, xi = round(xo.real, n), round(xo.imag, n)if xi == 0:return xrelse:return complex(xr, xi)def cardano_solution_v1(a, b, c, d):ab = -b/float(3*a)q = (3*a*c-(b**2)) / (9*(a**2))r = (9*a*b*c-27*(a**2)*d-2*(b**3)) / (54*(a**3))delta_sqrt = round((q**3+r**2)**(1.0/2), 4)# print("r, delta_sqrt:", r, delta_sqrt, round(r-delta_sqrt, 4))s = round((r+delta_sqrt)**(1.0/3), 8)t = round(r-delta_sqrt, 4)**(1.0/3)print("st,ab:", s, t, ab)imag = complex(0, (s-t)*(3**(1.0/3))/2)x1 = s+t+abx2 = -(s+t)/2+ab+imagx3 = -(s+t)/2+ab-imagreturn round_ri(x1), round_ri(x2), round_ri(x3)

是公式本身写的有误?还是浮点/开根号问题?需要再检查。。。

测试其他公式: 百度百科-一元三次方程求根公式

这里是说,一元三次方程的系数是复数时,用Cardano公式有问题(什么问题?),旧使用如下通用的求根公式:

def cardano_solution(a, b, c, d):#u = round((9*a*b*c-27*(a**2)*d-2*(b**3)) / (54*(a**3)), 4)#v = round(3*(4*a*c**3 - b**2*c**2-18*a*b*c*d+27*a**2*d**2+4*b**3*d) / (18**2*a**4), 4) ** (1.0/2)u = (9*a*b*c-27*(a**2)*d-2*(b**3)) / (54*(a**3))v = (3*(4*a*c**3 - b**2*c**2-18*a*b*c*d+27*a**2*d**2+4*b**3*d) / (18**2*a**4)) ** (1.0/2)if abs(u+v) >= abs(u-v):m = (u+v) ** (1.0/3)else:m = (u-v) ** (1.0/3)if m == 0: n == 0else:n = (b**2-3*a*c) / (9*a**2*m)# w = complex(0, -0.5+(3/4)**(1.0/2))# w2 = complex(0, -0.5-(3/4)**(1.0/2))w = -0.5+(-3/4)**(1.0/2)w2 = -0.5-(-3/4)**(1.0/2)ab = -b/float(3*a)x1 = m+n+abx2 = w*m+w2*n+abx3 = w2*m+w*n+ab# return x1, x2, x3return round_ri(x1), round_ri(x2), round_ri(x3)

测试结果cardano_solutioncardano_method包可以对上,且求解 x 3 + 1 = 0 x^3+1=0 x3+1=0 没有问题:

print(cardano_solution(1,3,4,4))print(cardano_solution(1,0,0,-1))print(cardano_solution(1,0,0,1))# ((-0.5+1.3229j), -2.0, (-0.5-1.3229j))# (1.0, (-0.5+0.866j), (-0.5-0.866j))# ((0.5+0.866j), -1.0, (0.5-0.866j))print(CubicEquation([1,3,4,4]).answers)print(CubicEquation([1,0,0,-1]).answers)print(CubicEquation([1,0,0,1]).answers) # 报错# [(-2+0j), (-0.5+1.322875j), (-0.5-1.322875j)]# [(1+0j), (-0.5+0.866025j), (-0.5-0.866025j)]# ZeroDivisionError ...

【总结】

根据公式自己编写的函数cardano_solution可求解一元三次方程。下载的包cardano_method和直接使用Cardano公式写的函数有问题。原因是?

附:python 获取复数的实部和虚部。使用j后缀表示虚数,比如a+bj中,a是实部,b是虚部,python中用complex(a,b)生成一个复数。

x = 2+1.5jprint(x.real) # 打印实部:2print(x.imag) # 打印虚部:2x1 = complex(2,1.5) # 使用`complex`生成复数 2+1.5j

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。