切割matlab的结果或精确的实数

Chopping the result of matlab or exact the reals

我有一个关于在MATLAB中求解四次方程的问题,请参阅

MATLAB函数如下所示:

MATLAB

function [res] = solveQuartic(a, b, c, d, e)
 p = (8*a*c - 3*b^2)/(8*a^2);
 q = (b^3 - 4*a*b*c + 8*a^2*d)/(8*a^3);

 delta0 = c^2-3*b*d + 12*a*e;
 delta1 = 2*c^3 - 9*b*c*d + 27*b^2*e + 27*a*d^2 - 72*a*c*e;
 Q = ((delta1 + sqrt(delta1^2 - 4*delta0^3))*0.5)^(1/3);
 S = 0.5*sqrt(-2/3*p + (Q + delta0/Q)/(3*a));

 x1 = -b/(4*a) - S - 0.5*sqrt(-4*S^2-2*p + q/S);
 x2 = -b/(4*a) - S + 0.5*sqrt(-4*S^2-2*p + q/S);
 x3 = -b/(4*a) + S - 0.5*sqrt(-4*S^2-2*p - q/S);
 x4 = -b/(4*a) + S + 0.5*sqrt(-4*S^2-2*p - q/S);
 res = [x1; x2; x3; x4];
end

测试

solveQuartic(-65, -312, -582, -488, -153)

(*
-1.400000000000000 - 0.627571632442189i
-1.000000000000000 + 0.000000000000000i
-1.400000000000000 + 0.627571632442189i
-1.000000000000000 + 0.000000000000000i
*)

显然,四次方程$$-65x^4-312x^3-582x^2-488x-153=0$$拥有两个real roots

现在我想求出四次方程的real roots,我的试验如下

function [res] = solveQuarticReals(a, b, c, d, e)
 p = (8*a*c - 3*b^2)/(8*a^2);
 q = (b^3 - 4*a*b*c + 8*a^2*d)/(8*a^3);

 delta0 = c^2-3*b*d + 12*a*e;
 delta1 = 2*c^3 - 9*b*c*d + 27*b^2*e + 27*a*d^2 - 72*a*c*e;
 Q = ((delta1 + sqrt(delta1^2 - 4*delta0^3))*0.5)^(1/3);
 S = 0.5*sqrt(-2/3*p + (Q + delta0/Q)/(3*a));

 x1 = -b/(4*a) - S - 0.5*sqrt(-4*S^2-2*p + q/S);
 x2 = -b/(4*a) - S + 0.5*sqrt(-4*S^2-2*p + q/S);
 x3 = -b/(4*a) + S - 0.5*sqrt(-4*S^2-2*p - q/S);
 x4 = -b/(4*a) + S + 0.5*sqrt(-4*S^2-2*p - q/S);
 res = [x1; x2; x3; x4];
 %exact the real roots
 j = 0;
  for i = 1 : length(res)
   if imag(res(i)) == 0
    j = j + 1;
    mid(j) = res(i); 
   end
  end
  res = mid;
end

然而,它失败了。

问题

a=-65;
b=-312;
c=-582;
d=-488;
e=-153;


p = (8*a*c - 3*b^2)/(8*a^2);
q = (b^3 - 4*a*b*c + 8*a^2*d)/(8*a^3);

delta0 = c^2-3*b*d + 12*a*e;
delta1 = 2*c^3 - 9*b*c*d + 27*b^2*e + 27*a*d^2 - 72*a*c*e;
Q = ((delta1 + sqrt(delta1^2 - 4*delta0^3))*0.5)^(1/3);
S = 0.5*sqrt(-2/3*p + (Q + delta0/Q)/(3*a));

x1 = -b/(4*a) - S - 0.5*sqrt(-4*S^2-2*p + q/S);
x2 = -b/(4*a) - S + 0.5*sqrt(-4*S^2-2*p + q/S);
x3 = -b/(4*a) + S - 0.5*sqrt(-4*S^2-2*p - q/S);
x4 = -b/(4*a) + S + 0.5*sqrt(-4*S^2-2*p - q/S);
res = [x1; x2; x3; x4];

%exact the real roots  
realNumber = real(res((abs(imag(res)) <= 1e-10)))

请尝试这些代码。 由于浮点数据类型的特性,imag(S)imag(0.5*sqrt(-4*S^2-2*p + q/S)) 之间可能存在一些差异,即使两个值在数学上是相同的。 为了检查错误,请在您的 matlab window 上键入 "format shortEng",然后键入 imag(S) - imag(0.5*sqrt(-4*S^2-2*p + q/S))。 答案是55.5112e-018i。 答案的虚数肯定不是0,而是55.5112e-018。 这种现象是由浮点数据类型的性质造成的。

PS。非常感谢 Stewie Griffin! :-)