Кратко — набор рабочих приёмов и замечаний о численной устойчивости при $a{x}^2+bx+c=0$ и малом $|b|$ по сравнению с $|a|$ и $|c|$.
1) Стандартная формула (опасность вычитания)
- Формула: $$x_{1,2}=\frac{-b\pm \sqrt{b^2-4ac}}{2a}.$$ - Проблема: если $|b|\ll \sqrt{|4ac|}$ и корни имеют противоположные знаки по модулю, то при вычислении того корня, где знак в числителе приводит к вычитанию схожих по величине чисел, возникает катастрофическое вычитание и потеря значащих цифр.
2) Устойчивое вычисление через вспомогательную величину $q$ - Вычислить дискриминант $D=b^2-4ac$ и $\sqrt{D}$ (для комплексных корней — комплексный корень).
- Возьмём $\mathrm{sign}(b)$ как $1$ если $b\ge 0$, иначе $-1$. Затем: $$q=-\frac{1}{2}(b+\mathrm{sign}(b)\sqrt{D}).$$ - Тогда корни: $$x_1=\frac{q}{a},x_2=\frac{c}{q}.$$ - Почему лучше: выбор знака в $q$ делает сумму в скобках большой по модулю, избегая вычитания близких чисел; второй корень получается через деление по произведению (устойчиво).
3) Использование Виетовых соотношений (если один корень вычислен точно)
- Если точно получили корень с большим модулем $x_{\text{big}}$, второй корень вычисляется как $$x_{\text{other}}=\frac{c}{a{x}_{\text{big}}}.$$ - Это часто даёт высокую относительную точность для второго корня, когда прямая формула даёт погрешность.
4) Балансировка (масштабирование переменной)
- Сделать уравнение моническим: $$x^2+\frac{b}{a}x+\frac{c}{a}=0.$$ - Масштабировать $x=y\cdot s$ с s=∣ca∣s=\sqrt{\left|\frac{c}{a}\right|}s= ac (чтобы константа стала $\pm 1$). Тогда коэффициенты становятся более сбалансированными и уменьшается риск вычитания при вычислении корней. После решения для $y$ перевести обратно $x=sy$.
5) Асимптотические/разложенческие приближения при очень малом $b$ - Для $b\to 0$ корни стремятся к $\pm \sqrt{-\frac{c}{a}}$ (если $ac<0$ — действительные) с малой поправкой. Можно использовать разложение в ряд по $b$ для начального приближения и затем уточнить (см. пункт 6).
6) Итеративное уточнение (Ньютон) для повышения точности
- Если получили приближённый корень $x^{(0)}$, одну итерацию Ньютона: $$x^{(1)}=x^{(0)}-\frac{a{x}^{(0)2}+b{x}^{(0)}+c}{2a{x}^{(0)}+b}$$ - Одна-две итерации часто восстанавливают полную машинную точность (если начальное приближение достаточно близко).
7) Комплексные корни
- Если $D<0$, вычисляйте $\sqrt{D}$ как комплексное число и применяйте тот же устойчивый приём с $q$. Для вещественной и мнимой частей лучше отделять: $$\mathrm{Re}(x)=-\frac{b}{2a},\mathrm{Im}(x)=\pm \frac{\sqrt{4ac-b^2}}{2|a|}.$$ - Следует избегать явного вычитания при вычислении малой реальной части; лучше получить её как $-b/(2a)$.
8) Применение повышенной точности или специальных алгоритмов
- Если проблема сильно плохо обусловлена (например, корни одинаковы по модулю и противоположны по знаку), единственный надёжный путь — арифметика с большей точностью (long double, quad) или использование специализированных библиотек/методов (например, алгоритмы для полиномиальных корней — Jenkins–Traub), которые дают устойчивые результаты.
Краткие указания по выборам в практике
- Используйте формулу с $q$ и затем $c/q$ — это самое простое и надёжное решение в большинстве случаев.
- Если один корень получается с плохой точностью, пересчитайте второй как $c/(a{x}_1)$.
- Если нужно гарантировать абсолютную/относительную точность при сильно неустойчивых коэффициентах — повышайте точность вычислений и/или применяйте итеративное уточнение.
Замечание о обусловленности
- При маленьком $|b|$ сумма корней $\frac{-b}{a}$ мала относительно их модулей (корни почти одинаковы по модулю, противоположны по знаку), поэтому задача найти оба корня с высокой относительной точностью одновременно плохо обусловлена: обычно можно получить оба с хорошей абсолютной точностью, но только один с хорошей относительной точностью; второй лучше получать через произведение.