0

我试图在 R 中解析三次方程的解，而不是数值解。

我在互联网上查找并获得了三次根的公式并编写了以下代码：链接是： http: //www.math.vanderbilt.edu/~schectex/courses/cubic/

cub <- function(a,b,c,d) {
  p <- -b/3/a
  q <- p^3 + (b*c-3*a*(d))/(6*a^2)
  r <- c/3/a
  x <- (q+(q^2+(r-p^2)^3)^0.5)^(1/3)+(q-(q^2+(r-p^2)^3)^0.5)^(1/3)+p
  x
}

然而这个函数在大多数情况下不起作用，我猜这是因为公式中负数的力量，例如我注意到 R 不能得到 (-8)^(1/3) 的实根，即 -2 . 但我不确定如何修复我的代码，以便它可以用来解决一般的精确三次解决方案。

谢谢。

4

4 回答 4

4

我会用polyroot(). 有关更多详细信息，请参见此处。

polyroot(z = c(8,0,0,1))
# [1]  1+1.732051i -2+0.000000i  1-1.732051i
于 2014-03-24T03:06:34.937 回答
3

试试这个：

# calcaulate -8 as a complex number
z <- as.complex(-8)  # or z <- -8 + 0i

# find all three cube roots
zroot3 <- z^(1/3) * exp(2*c(0:2)*1i*pi/3)
zroot3
## [1]  1+1.732051i -2+0.000000i  1-1.732051i


# check that all three cube roots cube to original
zroot3^3
## [1] -8+0i -8+0i -8-0i
于 2014-03-24T00:58:49.027 回答
1

如果您只想要真正的根目录，那么这里是另一种选择：

> x <- c( -8,8 )
> sign(x) * abs(x)^(1/3)
[1] -2  2

或者您可能对 Ryacas 包或 polynom 包感兴趣以获取其他选项。

于 2014-03-24T02:15:47.750 回答
0

这是一个计算所有解析解的函数： 'cubsol' 。任何意见都将受到欢迎。一个问题 - 目前，代码搜索效率相当低，真正的解决方案是... s2 = cuberoot(q-s0^0.5); 产生的三个复杂解决方案之一；xtemp[1:3] <- s1+ s2 +p; 有没有一种更有效的方法可以在计算之前知道它是哪一个？

# -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -
# Return all  the complex cube roots of a number
cuberoot <- function(x){
  return( as.complex(x)^(1/3)*exp(c(0,2,4)*1i*pi/3) );
}
# -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -   -
# cubsol solves analytically the cubic equation and
# returns a list whose first element is the real roots and the
# second element the complex roots.
# test with :
#a = -1; b=-10; c=0; d=50; x=0.01*(-1000:1500); plot(x,a*x^3+b*x^2+c*x+d,t='l'); abline(h=0)
# coefs = c(a,b,c,d)
cubsol <- function(coeffs) {
  if (!(length(coeffs) == 4)){
    stop('Please provide cubsol with a 4-vector of coefficients')
  }
  a = coeffs[1]; b=coeffs[2]; c=coeffs[3]; d=coeffs[4];
  rts = list();  

  p <- -b/3/a
  q <- p^3 + (b*c-3*a*(d))/(6*a^2)
  r <- c/3/a

  s0 = q^2+(r-p^2)^3; 
  xtemp = as.complex(rep(0,9));    
  if (s0 >= 0){ nReRts=1; } else {nReRts=3; }
  # Now find all the roots in complex space:
  s0 = as.complex(s0);
  s1 = cuberoot(q+s0^0.5)
  s2 = cuberoot(q-s0^0.5);
  xtemp[1:3] <- s1+ s2 +p;  # I think this is meant to always contain
                            # the sure real soln.
  # Second and third solution;
  iSqr3 = sqrt(3)*1i; 
  xtemp[4:6] = p - 0.5*(s1+s2 + iSqr3*(s1-s2));
  xtemp[7:9] = p - 0.5*(s1+s2 - iSqr3*(s1-s2));
  ind1 = which.min(abs(a*xtemp[1:3]^3 + b*xtemp[1:3]^2 +c*xtemp[1:3] +d))
  ind2 = 3+which.min(abs(a*xtemp[4:6]^3 + b*xtemp[4:6]^2 +c*xtemp[4:6] +d))
  ind3 = 6+which.min(abs(a*xtemp[7:9]^3 + b*xtemp[7:9]^2 +c*xtemp[7:9] +d))

  if (nReRts == 1){  
    rts[[1]] = c(Re(xtemp[ind1])); 
    rts[[2]] = xtemp[c(ind2,ind3)]
  } else {  # three real roots
    rts[[1]] = Re(xtemp[c(ind1,ind2,ind3)]); 
    rts[[2]] = numeric();
  }
  return(rts)

} # end of function cubsol
于 2018-06-17T00:28:09.367 回答