Nombres complexes D: Equation du troisième degré

Description

Bonjour,

Autre nom: Équation cubique (en: Cubic equation, de: Kubische Gleichung).

Selon WikipédiA, ces équations étaient connues des anciens Babyloniens, Grecs, Chinois, Indiens et Égyptiens.

La résolution est basée ici sur des méthodes algébriques classiques.

La forme générale de l'équation du troisième degré (ou cubique) est:
  a₃·x³ + a₂·x² + a₁·x + a₀ = 0;  avec a₃ ≠ 0.
Nous la désignerons par CE (cubic equation).
 
Dans la forme générale, remplaçons les aᵢ par bᵢ = aᵢ/a₃:
  ► x³ + b₂·x² + b₁·x + b₀ = 0
C'est l'équation cubique "simplifiée" (simplified cubic equation: SCE).

Posons x = y-c:
  ► y³ - 3·y²·c + 3·y·c² - c³ + b₂·(y² - 2·y·c + c²) + b₁·(y-c) + b₀ = 0
  ► y³ + (b₂ - 3·c)·y² + (b₁ - 2·b₂·c + 3·c²)·y + b₀ - b₁·c + b₂·c² - c³ = 0

En posant c = b₂/3, on annule le monôme de degré 2:
  ► y³ + p·y + q = 0  (depressed cubic equation: DCE)
avec:
  p = b₁ - 2·b₂·c + 3·c² = b₁ - 2·b₂²/3 + b₂²/3 = b₁ - b₂²/3
  q = b₀ - b₁·c + b₂·c² - c³ = b₀ - b₁·b₂/3 + b₂·b₂²/9 - b₂³/27 = b₀ - b₁·b₂/3 + 2·b₂³/27
void SCEtoDCE(const Cmplx b[3], Cmplx& p, Cmplx& q) {
  p = b[1] - b[2]*b[2]/3;
  q = b[0] - b[1]*b[2]/3 + b[2]*b[2]*b[2]/13.5;
} // simplified cubic equation to depressed cubic equation

L'équation DCE est bien de la forme:  y³ + p·y + q = 0, c'est-à-dire sans terme en y².

Si p = 0:
  b₁ = b₂²/3 et q = b₀ - b₁·b₂/3 + 2·b₂³/27 = b₀ - b₂³/9 + 2·b₂³/27 = b₀ - b₂³/27
  yᵢ = ³√-q  et  xᵢ = yᵢ - b₂/3

Si p ≠ 0:
posons  y = z - p/(3·z):
  ► z³ - 3·z²·p/(3·z) + 3·z·p²/(9·z²) - p³/(27·z³) + p·z - p²/(3·z) + q = 0
  ► z³ - p³/(27·z³) + q = 0
Amplifions l'équation par z³ = w:
  ► z⁶ + q·z³ - p³/27 = w² + q·w - p³/27 = 0
c'est une équation complexe du deuxième degré que savons résoudre:
  w = (√(q*q + p*p*p/6.75) - q)/2
  zᵢ = ³√w, yᵢ = z - p/(3·z)
Il suffit maintenant de calculer les 3 racines cubiques zᵢ de w, puis les valeurs des yᵢ et xᵢ.

Le code dans Equation.cpp utilise la classe Cmplx contenue dans les fichiers Cmplx.h et Cmplx.cpp, juste complétée des deux opérateurs de comparaison:
  bool operator==(const Cmplx y,const Reel d) {return y.R==d && y.I==0;}
  bool operator!=(const Cmplx y,const Reel d) {return y.R!=d || y.I!=0;}
et de la correction:
  Cmplx sqrt(const Cmplx& z)
    {if (z==0) return z; Reel y=sqrt(2*(abs(z)-z.R)); return Cmplx(z.I/y,y/2);}

Pas à pas, avec  a[3] = [-5+i·4], a[2] = [-6-i·3], a[1] = [1+i·7], a[0] = [2-i·8],
transformé en  b[2] = a[2]/a[3], b[1] = a[1]/a[3], b[0] = a[0]/a[3],
on effectue (d'une manière non optimisée) les calculs, impressions et contrôles correspondants à la "théorie" ci-dessus.

Ensuite, on utilise la fonction:
// simplified cubic complex equation: x*x*x + b[2]*x*x + b[1]*x + b[0] = 0
void ResolSCE(const Cmplx b[3], Cmplx x[3]) {
  Cmplx b2d3 = b[2]/3, p = b[1] - b[2]*b2d3;
  if (p == 0) { // b₀ - b₂³/27
    cbrt(b[0] - b[2]*b[2]*b[2]/27, x);
    x[0] -= b2d3; x[1] -= b2d3; x[2] -= b2d3;
  } else {
    Cmplx q = b[0] - b[1]*b2d3 + b[2]*b[2]*b[2]/13.5;
    cbrt((sqrt(q*q + p*p*p/6.75) - q)/2, x);
    x[0] -= (p/x[0]+b[2])/3; x[1] -= (p/x[1]+b[2])/3; x[2] -= (p/x[2]+b[2])/3;
  }
} // 3 solutions: x[0], x[1], x[2]
pour résoudre en une seule fois une équation donnée, sans entrer dans les détails.

Ce code est simplifié, optimisé et tient compte de la possibilité p == 0.
Cette éventualité apparait par exemple lorsque les trois solutions sont identiques; et ne pas en tenir compte, peut amener à des erreurs de calcul.

Observez, à la fin, les trois tests qui donnent 3 solutions distinctes, 2 identiques et 3 identiques:
Equation complexe du troisieme degre: a[3]*x*x*x + a[2]*x*x + a[1]*x + a[0] = 0
  avec: a[3] = [-5+i·4], a[2] = [-6-i·3], a[1] = [1+i·7], a[0] = [2-i·8]

Divison de l'equation par a[3]: x*x*x + b[2]*x*x + b[1]*x + b[0] = 0
  avec: b[2] = [0.439024+i·0.95122], b[1] = [0.560976-i·0.95122], b[0] = [-1.02439+i·0.780488]

Substitution: x = y - b[2]/3: y*y*y + p*y + q = 0
  avec: p = [0.798334-i·1.22963], q = [-1.4901+i·0.718808]

Substitution y = z - p/(3*z): z*z*z - p*p*p/(27*z*z*z) + q = 0
Substitution w = z*z*z: w - p*p*p/(27*w) + q = 0

Multiplication par w: w*w + q*w - p*p*p/27 = 0
Quadratic equation solution: w = [0.0574432+i·0.0432726]
  test: w*w + q*w - p*p*p/27 = [-8.32667e-17-i·1.249e-16]

Substitution: y[i] = z[i] - p/z[i]/3
  y[0] = [-0.00840461+i·1.18833], y[1] = [0.879511+i·0.116787], y[2] = [-0.871106-i·1.30512]
  test: y[0]*y[0]*y[0] + p*y[0] + q = [-1.77636e-15-i·1.11022e-15]
  test: y[1]*y[1]*y[1] + p*y[1] + q = [-1.9984e-15-i·9.99201e-16]
  test: y[2]*y[2]*y[2] + p*y[2] + q = [-1.9984e-15-i·1.11022e-15]

Substitution: x[i] = y[i] - b[2]/3
  x[0] = [-0.154746+i·0.871258], x[1] = [0.73317-i·0.200286], x[2] = [-1.01745-i·1.62219]
  test: x[0]*x[0]*x[0] + b[2]*x[0]*x[0] + b[1]*x[0] + b[0] = [-2.22045e-15-i·7.77156e-16]
  test: x[1]*x[1]*x[1] + b[2]*x[1]*x[1] + b[1]*x[1] + b[0] = [-2.22045e-15-i·9.99201e-16]
  test: x[2]*x[2]*x[2] + b[2]*x[2]*x[2] + b[1]*x[2] + b[0] = [-1.77636e-15-i·9.99201e-16]
  Cubic equation:
  test: a[3]*x[0]*x[0]*x[0] + a[2]*x[0]*x[0] + a[1]*x[0] + a[0] = [1.33227e-14-i·2.66454e-15]
  test: a[3]*x[1]*x[1]*x[1] + a[2]*x[1]*x[1] + a[1]*x[1] + a[0] = [1.46549e-14-i·4.44089e-15]
  test: a[3]*x[2]*x[2]*x[2] + a[2]*x[2]*x[2] + a[1]*x[2] + a[0] = [1.95399e-14-i·1.77636e-15]


Function ResolSCE(b, x): x*x*x + b[2]*x*x + b[1]*x + b[0] = 0
  avec: b[2] = [0.439024+i·0.95122], b[1] = [0.560976-i·0.95122], b[0] = [-1.02439+i·0.780488]
  x[0] = [-0.154746+i·0.871258], x[1] = [0.73317-i·0.200286], x[2] = [-1.01745-i·1.62219]
  test: x[0]*x[0]*x[0] + b[2]*x[0]*x[0] + b[1]*x[0] + b[0] = [-8.88178e-16-i·2.22045e-16]
  test: x[1]*x[1]*x[1] + b[2]*x[1]*x[1] + b[1]*x[1] + b[0] = [-8.88178e-16+i·2.22045e-16]
  test: x[2]*x[2]*x[2] + b[2]*x[2]*x[2] + b[1]*x[2] + b[0] = [-8.88178e-16+i·6.66134e-16]


Function ResolSCE(c, x), solution double: x*x*x + c[2]*x*x + c[1]*x + c[0] = 0
  (x-a[0])*(x-a[0])*(x-a[1]) = 0 ==> c[2] = [-5+i·9], c[1] = [56-i·20], c[0] = [-164+i·452]
  x[0] = [2-i·8], x[1] = [1+i·7], x[2] = [2-i·8]
  test: x[0]*x[0]*x[0] + c[2]*x[0]*x[0] + c[1]*x[0] + c[0] = [2.84217e-14-i·2.27374e-13]
  test: x[1]*x[1]*x[1] + c[2]*x[1]*x[1] + c[1]*x[1] + c[0] = [0+i·0]


Function ResolSCE(c, x), solution triple: x*x*x + c[2]*x*x + c[1]*x + c[0] = 0
  (x-a[0])*(x-a[0])*(x-a[0]) = 0 ==> c[2] = [-6+i·24], c[1] = [-180-i·96], c[0] = [376-i·416]
  x[0] = [2-i·8], x[1] = [2-i·8], x[2] = [2-i·8]
  test: x[0]*x[0]*x[0] + c[2]*x[0]*x[0] + c[1]*x[0] + c[0] = [0+i·0]


Bonne lecture ...
 

Liens

WikipédiA: Équation cubique
WikipediA: Cubic function
Sergei Khashin: Solution of cubic and quartic equations C++
CodeS-SourceS: Nombres complexes A: avec la bibliothèque logicielle std complex
CodeS-SourceS: Nombres complexes B: avec struct Cmplx
CodeS-SourceS: Nombres complexes C: avec class Cmplx
 

Codes Sources

A voir également

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.