Здесь представлен алгоритм для решения кубического уравнения методом Виета-Кардано. Программа написана для случая действительных коэффициентов (корни могут быть комплексными).
Кубическое уравнение записывается в виде:
x3+a*x2+b*x+c=0.
Для нахождения его корней, в случае действительных коэффициентов, вначале вычисляются:
Q=(a2-3b)/9, R=(2a3-9ab+27c)/54.
Далее, если R2<Q3, то уравнение имеет три действительных корня, вычисляющихся по формулам (Виета):
t=acos(R/sqrt(Q3))/3,
x1=-2*sqrt(Q)cos(t)-a/3,
x2=-2*sqrt(Q)cos(t+(2*pi/3))-a/3,
x3=-2*sqrt(Q)cos(t-(2*pi/3))-a/3.
В том случае, когда R2>=Q3, то действительных корней один (общий случай) или два (вырожденные случаи). Кроме действительного корня, имеется два комплексно-сопряженных. Для их нахождения вычисляются (формула Кардано):
A=-sign(R)[|R|+sqrt(R2-Q3)]1/3, B=Q/A при A!=0 или B=0 при A=0.
Действительный корень будет:
x1=(A+B)-a/3.
Комплексно-сопряженные корни:
x2,3=-(A+B)/2-a/3 + i*sqrt(3)*(A-B)/2
В том случае, когда A=B, то комплексно-сопряженные корни вырождаются в действительный: x2=-A-a/3.
Формулы Кардано и Виета требуют применения специальных функций, и в том случае, когда требуется провести большую серию вычислений корней кубического уравнения с не слишком сильно меняющимися коэффициентами, более быстрым алгоритмом является использование метода Ньютона или других итерационных методов (с нахождением начального приближения по формулам Кардано-Виета).
Ниже расположена программа для нахождения корней кубического уравнения с действительными коэффициентами.
/* Cubic equation solution. Real coefficients case.
int Cubic(double *x,double a,double b,double c); Parameters: x - solution array (size 3). On output: 3 real roots -> then x is filled with them; 1 real + 2 complex -> x[0] is real, x[1] is real part of complex roots, x[2] - non-negative imaginary part. a, b, c - coefficients, as described Returns: 3 - 3 real roots; 1 - 1 real root + 2 complex; 2 - 1 real root + complex roots imaginary part is zero (i.e. 2 real roots). */
#include <math.h> /* for sqrt(), fabs(), pow(), cos(), acos(). */ #define M_PI (3.141592653589793) #define M_2PI (2.*M_PI)
int Cubic(double *x,double a,double b,double c) { double q,r,r2,q3; q=(a*a-3.*b)/9.; r=(a*(2.*a*a-9.*b)+27.*c)/54.; r2=r*r; q3=q*q*q; if(r2<q3) { double t=acos(r/sqrt(q3)); a/=3.; q=-2.*sqrt(q); x[0]=q*cos(t/3.)-a; x[1]=q*cos((t+M_2PI)/3.)-a; x[2]=q*cos((t-M_2PI)/3.)-a; return(3); } else { double aa,bb; if(r<=0.) r=-r; aa=-pow(r+sqrt(r2-q3),1./3.); if(aa!=0.) bb=q/aa; else bb=0.; a/=3.; q=aa+bb; r=aa-bb; x[0]=q-a; x[1]=(-0.5)*q-a; x[2]=(sqrt(3.)*0.5)*fabs(r); if(x[2]==0.) return(2); return(1); } }
8 8 8
|