一元三次方程求根公式及其Fortran代码

类别:    标签: 数理   阅读次数:   版权: (CC) BY-NC-SA

2012-10-30 08:34:12

一元三次方程的求解比一元二次方程困难, 求根公式看起来也很复杂, 中文维基百科词条 上原先的三角函数解公式有误, 我已经修改过了. 但是我觉得还有必要将这些公式简化一下, 给出具体的算法, 以便编程使用.

首先需要明确, 一元三次方程 $f(x)=ax^3+bx^2+cx+d=0,(a\neq0)$ 至少有一个实根. 这是由于 $f(-\infty)=-\infty, f(\infty)=\infty$, 而 $f(x)$ 在整个区间上连续, 所以由微积分的介值定理, 必存在一点 $x_0$ 满足 $f(x_0)=0$ . 直观上也很容易理解, 曲线 $f(x)$ 取值范围为负无穷到正无穷, 那么它与x轴至少有一个交点, 这个交点就是方程的实根. 这个结论可以推广到所有奇数次的多项式方程.

对一元三次方程 $f(x)=ax^3+bx^2+cx+d=0, (a\neq0)$, 首先化为 $f(x)=x^3+{3b \over 3a} x^2+{6c \over 6a} x+{2d \over 2a}=x^3+3b’x^2+6c’x+2d’=0$ 的形式, 为简便我们仍记作 $f(x)=x^3+3bx^2+6cx+2d=0$ , 令

\[\begin{align} \alpha &=-b^3+3bc-d \\ \beta &=b^2-2c \\ \Delta &=\alpha^2-\beta^3 \\ R_1 &=\sqrt[3]{\alpha+\sqrt{\Delta}} \\ R_2 &=\sqrt[3]{\alpha-\sqrt{\Delta}}={\beta \over R_1} \end{align}\]

则方程的三个根可写作:

\[\begin{align} x_1 &=-b+R_1+R_2 \\ x_2 &=-b-{(R_1+R_2) \over 2} +{\sqrt{3}(R_1-R_2) \over 2} i \\ x_3 &=-b-{(R_1+R_2) \over 2} -{\sqrt{3}(R_1-R_2) \over 2} i \end{align}\]

其中 $\Delta$ 即为三次方程根的判别式, 由其容易看出:

很显然, $\Delta \lt 0$ 的情况最复杂, 也很让人迷惑, 因为实根必须借助于复数才能求得, 这也是最初要引入复数的原因. 由于我没有学过复变函数, 下面的说法可能有误, 说出来只是为了帮助大家理解.

当 $\Delta \lt 0$ 时, 记 $z=\alpha+\sqrt{|\Delta|}i, \bar z=\alpha-\sqrt{|\Delta|}i, R_1=\sqrt[3]{z}, R_2=\sqrt[3]{\bar z}$, 若我们认为 $\sqrt[3]{\bar z}=\overline {\sqrt[3]{z} }$ , 则 $R_2=\bar {R_1}$, 记 $R_1=R, R_2=\bar R, \omega= { {-1+\sqrt{3} i} \over 2}, \bar \omega={ {-1-\sqrt{3} i} \over 2}$ , 则三个根可写作:

\[\begin{align} x_1 &=-b+R_1+R_2=-b+R+\bar R \\ x_2 &=-b-{(R_1+R_2) \over 2} +{\sqrt{3}(R_1-R_2) \over 2} i =-b+\omega R+\bar \omega \bar R \\ x_3 &=-b-{(R_1+R_2) \over 2} -{\sqrt{3}(R_1-R_2) \over 2} i =-b+\bar \omega R+ \omega \bar R \end{align}\]

很显然

\[\begin{align} x_1 &=-b+2\Re(R) \\ \Re(R) &=\Re \sqrt[3]{z} = \sqrt[3]{|z|} \cos({ {\theta+2k\pi} \over 3}), k=0, 1, -1 \\ |z| &=\sqrt{\alpha^2+|\Delta|}=\sqrt{\alpha^2-\alpha^2+\beta^3} = \beta^{3/2} \\ \theta &= \arccos({\alpha \over |z|}) = \arccos({\alpha \over \beta^{3/2}}) \end{align}\]

由于复变函数的多值性, 我们一下就得到了所有的三个根, 正如上面所给出的. 另外由于 $\omega$ 对应的正是旋转 $2\pi \over 3$ 的复数, 所以 $x_2, x_3$ 并不能给出新的根.

代码

# Language: fortran
Subroutine getCubicRoot(P, X)
real*8, parameter:: TwoPi = 8.D0*atan(1.D0)
real*8 P(*), a, b, c, d, Alph, Beta, Delt, R1, R2, tht, X(3)

X = 0.D0
a = P(1)
b = P(2)/(3.D0*a)
c = P(3)/(6.D0*a)
d = P(4)/(2.D0*a)

Alph = -b*b*b + 3.D0*b*c - d
Beta =  b*b - 2.D0*c
Delt = Alph*Alph-Beta*Beta*Beta

if(Delt>0.D0) then
	tht = Alph+sqrt(Delt); R1 = sign(abs(tht)**(1.D0/3.D0), tht)
	tht = Alph-sqrt(Delt); R2 = sign(abs(tht)**(1.D0/3.D0), tht)
	X(1) = -b+R1+R2
else if(Delt==0.D0) then
	R1 = sign(abs(Alph)**(1.D0/3.D0), Alph)
	if(R1==0.D0) then
		X(1) = -b
	else
		X(1) = -b+2.D0*R1
		X(2) = -b-R1
	end if
else if(Delt<0.D0) then
	tht = acos(Alph/(sqrt(Beta)*Beta))
	X(1)  = -b+2.D0*sqrt(Beta)*cos(tht/3.D0)
	X(2)  = -b+2.D0*sqrt(Beta)*cos((tht+TwoPi)/3.D0)
	X(3)  = -b+2.D0*sqrt(Beta)*cos((tht-TwoPi)/3.D0)
end if
End Subroutine getCubicRoot
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 【转】Fortran的奇淫技巧
后一篇: GnuPlot:简单数据统计处理及取整函数int/floor的问题

访问人次(2015年7月 9日起): | 最后更新: 2024-01-20 10:40:28 UTC | 版权所有 © 2008 - 2024 Jerkwin