Представитель Шуры Люберецкого в ЖЖ (brat_luber) wrote,
Представитель Шуры Люберецкого в ЖЖ
brat_luber

Вы делаете это неправильно

По-моему, почти в каждый курс программирования входит задачка вроде “напишите программу, решающую квадратные уравнения”. Обычно это второе или третье задание после “Hello, world!” – считается, что это хороший способ продемонстрировать нетривиальные инструкции ветвления. “Хорошее” решение сводится к вычислению дискриминанта и в случае, если он неотрицательный – вычислению корней квадратного трехчлена ax2+bx+c по формуле, известной из курса алгебры за седьмой, что ли, класс:

square1

Вот такой вариант решения обычно считается более-менее приемлемым (хотя его можно/нужно обвешать еще несколькими проверками – например, не равен ли коэффициент a нулю?):

int solve(double a, double b, double c, double *x1, double *x2){
	double d = b*b - 4.*a*c;
	if( d >= 0 ){
		d = sqrt(d);
		*x1 = (-b + d)/(2.*a);
		*x2 = (-b - d)/(2.*a);
		return 0;
	}
	return -1;
}

В чем проблема? На первый взгляд все более-менее хорошо, но… Давайте для тестирования будем подсовывать уравнения с известными корнями – используя для этого теорему Виета. А именно, зафиксируем коэффициент a=1, тогда уравнение с корнями x1 и x2 будет иметь коээффициенты b=-(x1+x2) и c=x1x2. Сравнивая корни, полученные при решении уравнения, с известными нам, оценим “качество” решения.

Если корни “нормальные” – те, с которыми справится шестиклассник – то все хорошо. Но что будет, если взять два “нехороших” корня – к примеру, x1=1, x2=10-14 (это для double; если вы пользуетесь float – то возьмите второй корень, равный 10-5)? Проверьте – не забыв включить вывод максимально возможного количества значащих цифр (в printf лучше всего использовать форматный спецификатор %g, при использовании вывода в стиле C++, через iostream, такой вывод включен по умолчанию). Ошибка при вычислении второго корня возникнет уже в четвертой значащей цифре, это, на самом деле, уже довольно неприлично.

В чем дело? Если обратить внимание на “обычную” формулу для вычисления корней квадратного уравнения, то мы увидим, что при вычислении меньшего корня корень из дискриминанта вычитается из b – а так как они в этом случае очень близки, то возникающая при этом вычислительная ошибка оказывается слишком большой.

Метод, разумеется, можно улучшить. Для начала – можно вспомнить о существовании еще одной формулы для корней квадратного уравнения:

square2

Выводится она абсолютно аналогичным образом, от “классической” отличается тем, что “не работает” при c=0.

Если переписать программу, чтобы она использовала эту формулу – то меньший корень “нехорошего” уравнения будет вычисляться точно, а проблемы возникнут с большим корнем. Причина та же самая – вычитание двух близких по величине чисел. Но ведь если вычислять больший корень по первой формуле, а меньший – по второй, то эта проблема исчезнет! Поэтому более правильный метод решения квадратного уравнения должен выглядеть так:

- вычисляем дискриминант D=b2-4ac
- если дискриминант неотрицателен, то вычисляем

q

- корни уравнения равны q/a и c/q.

Как реализовать это в программе – довольно очевидно, это особо ее не усложнит.

Какая здесь мораль? Численные методы и программирование – это две совершенно разных области человеческого знания. Математические задачки – вроде решения квадратного уравнения – могут показаться интересными с точки зрения обучения программированию, но это “чужая территория”, и можно столкнуться с совершенно непредсказуемыми вещами. Признайтесь, многие ли слышали о сложностях, возникающих при решении на компьютере квадратных уравнений – хотя казалось бы, что может быть проще?

Запись опубликована в блоге Шуры Люберецкого. Вы можете оставлять свои комментарии там, используя свое имя пользователя из ЖЖ (вход по OpenID).

Subscribe

Recent Posts from This Journal

Comments for this post were disabled by the author