La résolution d'une équation cubique
https://stackoverflow.com/questions/1829330
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
Dans le cadre d'un programme que je vous écris, je dois résoudre une équation cubique exactement (plutôt que d'utiliser un viseur numérique racine):
a*x**3 + b*x**2 + c*x + d = 0.
Je suis en train d'utiliser les équations de . Cependant, considérez le code suivant (ce qui est Python, mais il est un code assez générique):
a = 1.0
b = 0.0
c = 0.2 - 1.0
d = -0.7 * 0.2
q = (3*a*c - b**2) / (9 * a**2)
r = (9*a*b*c - 27*a**2*d - 2*b**3) / (54*a**3)
print "q = ",q
print "r = ",r
delta = q**3 + r**2
print "delta = ",delta
# here delta is less than zero so we use the second set of equations from the article:
rho = (-q**3)**0.5
# For x1 the imaginary part is unimportant since it cancels out
s_real = rho**(1./3.)
t_real = rho**(1./3.)
print "s [real] = ",s_real
print "t [real] = ",t_real
x1 = s_real + t_real - b / (3. * a)
print "x1 = ", x1
print "should be zero: ",a*x1**3+b*x1**2+c*x1+d
Mais la sortie est:
q = -0.266666666667
r = 0.07
delta = -0.014062962963
s [real] = 0.516397779494
t [real] = 0.516397779494
x1 = 1.03279555899
should be zero: 0.135412149064
la sortie est non nul, et donc x1 est pas vraiment une solution. Y at-il une erreur dans l'article de Wikipedia?
ps: Je sais que numpy.roots va résoudre ce genre d'équation, mais je dois le faire pour des millions d'équations et donc je dois mettre en œuvre ce pour travailler sur des tableaux de coefficients
.La solution
La (rho^(1/3), theta/3) de notation de Wikipédia ne signifie pas que rho^(1/3) est la partie réelle et theta/3 est la partie imaginaire. Au contraire, ceci est en coordonnées polaires. Ainsi, si vous voulez que la partie réelle, vous prendriez rho^(1/3) * cos(theta/3).
J'ai fait ces changements à votre code et cela a fonctionné pour moi:
theta = arccos(r/rho)
s_real = rho**(1./3.) * cos( theta/3)
t_real = rho**(1./3.) * cos(-theta/3)
(Bien sûr, s_real = t_real ici parce que cos est encore.)
Autres conseils
Je l'ai regardé l'article de Wikipedia et votre programme.
J'ai aussi résolu l'équation en utilisant Wolfram Alpha et les résultats ne correspondent pas là-bas ce que vous obtenez.
Je venais de passer par votre programme à chaque étape, utilise beaucoup de déclarations d'impression et obtenir chaque résultat intermédiaire. Ensuite, passer par une calculatrice et de le faire vous-même.
Je ne peux pas trouver ce qui se passe, mais où vos calculs à la main et le programme est un bon diverger endroit pour regarder.
Voici la solution de A. Rex en JavaScript:
a = 1.0;
b = 0.0;
c = 0.2 - 1.0;
d = -0.7 * 0.2;
q = (3*a*c - Math.pow(b, 2)) / (9 * Math.pow(a, 2));
r = (9*a*b*c - 27*Math.pow(a, 2)*d - 2*Math.pow(b, 3)) / (54*Math.pow(a, 3));
console.log("q = "+q);
console.log("r = "+r);
delta = Math.pow(q, 3) + Math.pow(r, 2);
console.log("delta = "+delta);
// here delta is less than zero so we use the second set of equations from the article:
rho = Math.pow((-Math.pow(q, 3)), 0.5);
theta = Math.acos(r/rho);
// For x1 the imaginary part is unimportant since it cancels out
s_real = Math.pow(rho, (1./3.)) * Math.cos( theta/3);
t_real = Math.pow(rho, (1./3.)) * Math.cos(-theta/3);
console.log("s [real] = "+s_real);
console.log("t [real] = "+t_real);
x1 = s_real + t_real - b / (3. * a);
console.log("x1 = "+x1);
console.log("should be zero: "+(a*Math.pow(x1, 3)+b*Math.pow(x1, 2)+c*x1+d));
Ici, je mets un solveur équation cubique (avec des coefficients complexes).
#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>
using namespace std;
#define PI 3.141592
long double complex_multiply_r(long double xr, long double xi, long double yr, long double yi) {
return (xr * yr - xi * yi);
}
long double complex_multiply_i(long double xr, long double xi, long double yr, long double yi) {
return (xr * yi + xi * yr);
}
long double complex_triple_multiply_r(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi) {
return (xr * yr * zr - xi * yi * zr - xr * yi * zi - xi * yr * zi);
}
long double complex_triple_multiply_i(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi) {
return (xr * yr * zi - xi * yi * zi + xr * yi * zr + xi * yr * zr);
}
long double complex_quadraple_multiply_r(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi, long double wr, long double wi) {
long double z1r, z1i, z2r, z2i;
z1r = complex_multiply_r(xr, xi, yr, yi);
z1i = complex_multiply_i(xr, xi, yr, yi);
z2r = complex_multiply_r(zr, zi, wr, wi);
z2i = complex_multiply_i(zr, zi, wr, wi);
return (complex_multiply_r(z1r, z1i, z2r, z2i));
}
long double complex_quadraple_multiply_i(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi, long double wr, long double wi) {
long double z1r, z1i, z2r, z2i;
z1r = complex_multiply_r(xr, xi, yr, yi);
z1i = complex_multiply_i(xr, xi, yr, yi);
z2r = complex_multiply_r(zr, zi, wr, wi);
z2i = complex_multiply_i(zr, zi, wr, wi);
return (complex_multiply_i(z1r, z1i, z2r, z2i));
}
long double complex_divide_r(long double xr, long double xi, long double yr, long double yi) {
return ((xr * yr + xi * yi) / (yr * yr + yi * yi));
}
long double complex_divide_i(long double xr, long double xi, long double yr, long double yi) {
return ((-xr * yi + xi * yr) / (yr * yr + yi * yi));
}
long double complex_root_r(long double xr, long double xi) {
long double r, theta;
r = sqrt(xr*xr + xi*xi);
if (r != 0.0) {
if (xr >= 0 && xi >= 0) {
theta = atan(xi / xr);
}
else if (xr < 0 && xi >= 0) {
theta = PI - abs(atan(xi / xr));
}
else if (xr < 0 && xi < 0) {
theta = PI + abs(atan(xi / xr));
}
else {
theta = 2.0 * PI + atan(xi / xr);
}
return (sqrt(r) * cos(theta / 2.0));
}
else {
return 0.0;
}
}
long double complex_root_i(long double xr, long double xi) {
long double r, theta;
r = sqrt(xr*xr + xi*xi);
if (r != 0.0) {
if (xr >= 0 && xi >= 0) {
theta = atan(xi / xr);
}
else if (xr < 0 && xi >= 0) {
theta = PI - abs(atan(xi / xr));
}
else if (xr < 0 && xi < 0) {
theta = PI + abs(atan(xi / xr));
}
else {
theta = 2.0 * PI + atan(xi / xr);
}
return (sqrt(r) * sin(theta / 2.0));
}
else {
return 0.0;
}
}
long double complex_cuberoot_r(long double xr, long double xi) {
long double r, theta;
r = sqrt(xr*xr + xi*xi);
if (r != 0.0) {
if (xr >= 0 && xi >= 0) {
theta = atan(xi / xr);
}
else if (xr < 0 && xi >= 0) {
theta = PI - abs(atan(xi / xr));
}
else if (xr < 0 && xi < 0) {
theta = PI + abs(atan(xi / xr));
}
else {
theta = 2.0 * PI + atan(xi / xr);
}
return (pow(r, 1.0 / 3.0) * cos(theta / 3.0));
}
else {
return 0.0;
}
}
long double complex_cuberoot_i(long double xr, long double xi) {
long double r, theta;
r = sqrt(xr*xr + xi*xi);
if (r != 0.0) {
if (xr >= 0 && xi >= 0) {
theta = atan(xi / xr);
}
else if (xr < 0 && xi >= 0) {
theta = PI - abs(atan(xi / xr));
}
else if (xr < 0 && xi < 0) {
theta = PI + abs(atan(xi / xr));
}
else {
theta = 2.0 * PI + atan(xi / xr);
}
return (pow(r, 1.0 / 3.0) * sin(theta / 3.0));
}
else {
return 0.0;
}
}
void main() {
long double a[2], b[2], c[2], d[2], minusd[2];
long double r, theta;
cout << "ar?";
cin >> a[0];
cout << "ai?";
cin >> a[1];
cout << "br?";
cin >> b[0];
cout << "bi?";
cin >> b[1];
cout << "cr?";
cin >> c[0];
cout << "ci?";
cin >> c[1];
cout << "dr?";
cin >> d[0];
cout << "di?";
cin >> d[1];
if (b[0] == 0.0 && b[1] == 0.0 && c[0] == 0.0 && c[1] == 0.0) {
if (d[0] == 0.0 && d[1] == 0.0) {
cout << "x1r: 0.0 \n";
cout << "x1i: 0.0 \n";
cout << "x2r: 0.0 \n";
cout << "x2i: 0.0 \n";
cout << "x3r: 0.0 \n";
cout << "x3i: 0.0 \n";
}
else {
minusd[0] = -d[0];
minusd[1] = -d[1];
r = sqrt(minusd[0]*minusd[0] + minusd[1]*minusd[1]);
if (minusd[0] >= 0 && minusd[1] >= 0) {
theta = atan(minusd[1] / minusd[0]);
}
else if (minusd[0] < 0 && minusd[1] >= 0) {
theta = PI - abs(atan(minusd[1] / minusd[0]));
}
else if (minusd[0] < 0 && minusd[1] < 0) {
theta = PI + abs(atan(minusd[1] / minusd[0]));
}
else {
theta = 2.0 * PI + atan(minusd[1] / minusd[0]);
}
cout << "x1r: " << pow(r, 1.0 / 3.0) * cos(theta / 3.0) << "\n";
cout << "x1i: " << pow(r, 1.0 / 3.0) * sin(theta / 3.0) << "\n";
cout << "x2r: " << pow(r, 1.0 / 3.0) * cos((theta + 2.0 * PI) / 3.0) << "\n";
cout << "x2i: " << pow(r, 1.0 / 3.0) * sin((theta + 2.0 * PI) / 3.0) << "\n";
cout << "x3r: " << pow(r, 1.0 / 3.0) * cos((theta + 4.0 * PI) / 3.0) << "\n";
cout << "x3i: " << pow(r, 1.0 / 3.0) * sin((theta + 4.0 * PI) / 3.0) << "\n";
}
}
else {
// find eigenvalues
long double term0[2], term1[2], term2[2], term3[2], term3buf[2];
long double first[2], second[2], second2[2], third[2];
term0[0] = -4.0 * complex_quadraple_multiply_r(a[0], a[1], c[0], c[1], c[0], c[1], c[0], c[1]);
term0[1] = -4.0 * complex_quadraple_multiply_i(a[0], a[1], c[0], c[1], c[0], c[1], c[0], c[1]);
term0[0] += complex_quadraple_multiply_r(b[0], b[1], b[0], b[1], c[0], c[1], c[0], c[1]);
term0[1] += complex_quadraple_multiply_i(b[0], b[1], b[0], b[1], c[0], c[1], c[0], c[1]);
term0[0] += -4.0 * complex_quadraple_multiply_r(b[0], b[1], b[0], b[1], b[0], b[1], d[0], d[1]);
term0[1] += -4.0 * complex_quadraple_multiply_i(b[0], b[1], b[0], b[1], b[0], b[1], d[0], d[1]);
term0[0] += 18.0 * complex_quadraple_multiply_r(a[0], a[1], b[0], b[1], c[0], c[1], d[0], d[1]);
term0[1] += 18.0 * complex_quadraple_multiply_i(a[0], a[1], b[0], b[1], c[0], c[1], d[0], d[1]);
term0[0] += -27.0 * complex_quadraple_multiply_r(a[0], a[1], a[0], a[1], d[0], d[1], d[0], d[1]);
term0[1] += -27.0 * complex_quadraple_multiply_i(a[0], a[1], a[0], a[1], d[0], d[1], d[0], d[1]);
term1[0] = -27.0 * complex_triple_multiply_r(a[0], a[1], a[0], a[1], d[0], d[1]);
term1[1] = -27.0 * complex_triple_multiply_i(a[0], a[1], a[0], a[1], d[0], d[1]);
term1[0] += 9.0 * complex_triple_multiply_r(a[0], a[1], b[0], b[1], c[0], c[1]);
term1[1] += 9.0 * complex_triple_multiply_i(a[0], a[1], b[0], b[1], c[0], c[1]);
term1[0] -= 2.0 * complex_triple_multiply_r(b[0], b[1], b[0], b[1], b[0], b[1]);
term1[1] -= 2.0 * complex_triple_multiply_i(b[0], b[1], b[0], b[1], b[0], b[1]);
term2[0] = 3.0 * complex_multiply_r(a[0], a[1], c[0], c[1]);
term2[1] = 3.0 * complex_multiply_i(a[0], a[1], c[0], c[1]);
term2[0] -= complex_multiply_r(b[0], b[1], b[0], b[1]);
term2[1] -= complex_multiply_i(b[0], b[1], b[0], b[1]);
term3[0] = complex_multiply_r(term1[0], term1[1], term1[0], term1[1]);
term3[1] = complex_multiply_i(term1[0], term1[1], term1[0], term1[1]);
term3[0] += 4.0 * complex_triple_multiply_r(term2[0], term2[1], term2[0], term2[1], term2[0], term2[1]);
term3[1] += 4.0 * complex_triple_multiply_i(term2[0], term2[1], term2[0], term2[1], term2[0], term2[1]);
term3buf[0] = term3[0];
term3buf[1] = term3[1];
term3[0] = complex_root_r(term3buf[0], term3buf[1]);
term3[1] = complex_root_i(term3buf[0], term3buf[1]);
if (term0[0] == 0.0 && term0[1] == 0.0 && term1[0] == 0.0 && term1[1] == 0.0) {
cout << "x1r: " << -pow(d[0], 1.0 / 3.0) << "\n";
cout << "x1i: " << 0.0 << "\n";
cout << "x2r: " << -pow(d[0], 1.0 / 3.0) << "\n";
cout << "x2i: " << 0.0 << "\n";
cout << "x3r: " << -pow(d[0], 1.0 / 3.0) << "\n";
cout << "x3i: " << 0.0 << "\n";
}
else {
// eigenvalue1
first[0] = complex_divide_r(complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1]), 3.0 * pow(2.0, 1.0 / 3.0) * a[0], 3.0 * pow(2.0, 1.0 / 3.0) * a[1]);
first[1] = complex_divide_i(complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1]), 3.0 * pow(2.0, 1.0 / 3.0) * a[0], 3.0 * pow(2.0, 1.0 / 3.0) * a[1]);
second[0] = complex_divide_r(pow(2.0, 1.0 / 3.0) * term2[0], pow(2.0, 1.0 / 3.0) * term2[1], 3.0 * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
second[1] = complex_divide_i(pow(2.0, 1.0 / 3.0) * term2[0], pow(2.0, 1.0 / 3.0) * term2[1], 3.0 * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
cout << "x1r: " << first[0] - second[0] - third[0] << "\n";
cout << "x1i: " << first[1] - second[1] - third[1] << "\n";
// eigenvalue2
first[0] = complex_divide_r(complex_multiply_r(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
first[1] = complex_divide_i(complex_multiply_r(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
second[0] = complex_divide_r(complex_multiply_r(1.0, sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
second[1] = complex_divide_i(complex_multiply_r(1.0, sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
cout << "x2r: " << -first[0] + second[0] - third[0] << "\n";
cout << "x2i: " << -first[1] + second[1] - third[1] << "\n";
// eigenvalue3
first[0] = complex_divide_r(complex_multiply_r(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
first[1] = complex_divide_i(complex_multiply_r(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
second[0] = complex_divide_r(complex_multiply_r(1.0, -sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, -sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
second[1] = complex_divide_i(complex_multiply_r(1.0, -sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, -sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
cout << "x3r: " << -first[0] + second[0] - third[0] << "\n";
cout << "x3i: " << -first[1] + second[1] - third[1] << "\n";
}
}
int end;
cin >> end;
}
Dans le cas quelqu'un a besoin de code C, vous pouvez utiliser ce morceau de OpenCV:
https://github.com/opencv /opencv/blob/master/modules/calib3d/src/polynom_solver.cpp
Le programme ci-dessous permet de résoudre parfaitement l'équation cubique de la forme Ax ^ 3 + Bx ^ 2 + cx + d.
Le code impeccable est écrit en C # .Net.
Les racines réelles et imaginaires sont séparés par des couleurs noires et rouges respectivement.
mettre à jour simplement la valeur des variables A, B, C et D pour voir les réponses à-dire la solution de l'équation.
Pour trouver les racines 2, solution d'une équation Quadratique s'il vous plaît visitez ici.
Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
