Решение квадратного уравнения

Получение корней квадратного уравнения относится к элементарным понятиям алгебры; формулы для корней известны со времен царя Гороха. Однако, при численной реализации этих формул могут встретиться неприятности, обход одной из которых будет здесь разобран. Программа написана для случая действительных коэффициентов (корни могут быть комплексными).

Квадратное уравнение записывается в виде:
a*x2+b*x+c=0,   a!=0.
Классическая формула для нахождения его корней (как в действительном, так и в комплексном случае):

x1=(-b-sqrt(b2-4ac))/(2a),    x2=(-b+sqrt(b2-4ac))/(2a).
Существует, однако, и альтернативный способ получения корней:
x1=2c/(-b-sqrt(b2-4ac)),    x2=2c/(-b+sqrt(b2-4ac)).
Трудность, о которой будет идти здесь речь (она не единственная), связана с тем, что в том случае, когда b2 значительно превышает |4ac|, то при пользовании классическими формулами один из корней будет формироваться с использованием вычитания двух близких чисел, что может привести к потере точности. Чтобы этого избежать, по классическим формулам вычисляется только один корень (при вычислении которого этого не происходит), а второй корень находится по формуле Виета.

Случай, когда корни комплексные, для действительных корней уравнения указанной трудности не представляет.

Ниже расположена программа для решения квадратного уравнения с действительными коэффициентами.

/* Quadratic equation solution. Real coefficients case.

   int Quadratic(double *x,double a,double b,double c);
   Parameters:
   x - solution array (size 2). On output:
       2 real roots -> then x is filled with them;
       2 complex-conjugate roots -> x[0] is real part, 
         x[1] is non-negative imaginary part. 
       other cases -> x[0] unique valid root if (-1) returned,
                     no valid roots otherwise.
   a, b, c - coefficients: ax^2 + bx + c = 0. 
   Returns: 2 - 2 real and distince roots;
            1 - 1 real distinct root (x[0]=x[1]);
            0 - 2 complex roots;
            -1 - one real root in case a==0;
            -2 - no roots in case a=0, b=0;
            -3 - infinite number of roots (a=b=c=0). 
*/

#include <math.h>   /* for sqrt() */

int Quadratic(double *x,double a,double b,double c) {
  double d;
  /* degenerated cases */
  if(a==0.) {
    if(b==0.) {
      if(c==0.) return(-3);
      return(-2);
    }
    x[0]=-c/b; return(-1);
  }
  /* the main case */
  d=b*b-4.*a*c;       /* the discriminant */
  /* one distinct root */
  if(d==0.) {
    x[0]=x[1]=-b/(2.*a); return(1);
  }
  /* conjugate complex roots */
  if(d<0.) {
    double t=0.5/a;
    x[0]=-b*t; x[1]=sqrt(-d)*t;
    return(0);
  }
  /* 2 real roots: avoid subtraction of 2 close numbers */
  if(b>=0.) d=(-0.5)*(b+sqrt(d));
  else d=(-0.5)*(b-sqrt(d));
  x[0]=d/a; x[1]=c/d;
  return(2);
}
Скачать текст программы quadratic.c
Полиномы
Нелинейные уравнения
Индекс
Hosted by uCoz