Кратко — устойчивые стратегии, алгоритм и оценка погрешности.
1) Устойчивый численный алгоритм (однопроходный, избегает катастрофического вычитания)
- Вычислить дискриминант и его корень:
$$D=b^2-4ac,\sqrt{D}=\mathrm{sqrt}(D)(\text{в комплексной области при}D<0).$$ - Выбрать знак так, чтобы избежать вычитания близких по модулю чисел:
$$q=-\frac{1}{2}(b+\mathrm{sign}(b)\sqrt{D}),\text{(при}b=0\mathrm{sign}(b):=1).$$ Тогда корни вычисляются как
$$r_1=\frac{q}{a},r_2=\frac{c}{q}.$$ Если $q=0$ (двойной корень), взять $r=-\frac{b}{2a}$.
Пояснение: формула через $q$ эквивалентна классическому $(-b\pm \sqrt{D})/(2a)$, но устраняет потерю значащих цифр при вычитании большого близких по модулю слагаемых.
2) Нюансы масштабирования и точности
- Нормализация на $a$ (приведение к моническому многочлену $x^2+(b/a)x+(c/a)=0$) снижает число операций, но при очень малых $a$ (в интервале $(0,1]$ $a$ может быть очень мала) деление усиливает погрешности. Если $a$ не чрезмерно мало — нормализация допустима; при малых $a$ лучше вычислять без предварительного умножения/деления, использовать удвоенную точность или компенсированные операции (FMA).
- Для комплексных корней применяют ту же схему с комплексными sqrt и sign.
3) Оценка погрешностей (линейная аппроксимация)
- Для малого приращения коэффициентов $\Delta a,\Delta b,\Delta c$ линейная аппроксимация изменения корня $r$ даёт
$$\Delta r\approx -\frac{r^2\Delta a+r\Delta b+\Delta c}{2ar+b},$$ поэтому абсолютная оценка
$$|\Delta r|\le \frac{|r{|}^2|\Delta a|+|r||\Delta b|+|\Delta c|}{|2ar+b|}.$$ Для относительных возмущений коэффициентов можно перейти к соответствующим формам (подставить $\Delta a=a\delta a$ и т.д.).
- Оценка через обусловленность (приближённо): если вычисления имеют относительную ошибку порядка машинного эпсилон $\epsilon$, то
$$|\Delta r|\lesssim \kappa (r)\epsilon ,$$ где приближённая обусловленность
$$\kappa (r)\approx \frac{|a||r{|}^2+|b||r|+|c|}{|2ar+b|}.$$ Это даёт оценку прямой погрешности через известное относительное погрешности коэффициентов/вычислений.
4) Практические рекомендации для данного интервала $a\in (0,1],b\in [-2,2],c\in [-1,1]$ - Всегда использовать формулу с $q$ для минимизации потерь значащих цифр.
- Вычислять $b^2$ и $4ac$ с FMA, если доступно, чтобы уменьшить погрешность дискриминанта.
- При очень малом $|2ar+b|$ (производная многочлена близка к нулю в корне) корень плохо обусловлен — требуется увеличение точности (quad, long double) или метод интервальной арифметики для гарантированных границ.
- Для гарантированных интервалов корней применять интервальную арифметику (outward rounding) или методы Кравчик/интервальный Ньютона: сначала получать численное приближение устойчивым алгоритмом, затем уточнять и контролировать радиусы интервалов.
5) Итоговая схема устойчивого вычисления и оценки
- Вычислить $D=b^2-4ac$ (FMA), $\sqrt{D}$.
- Найти $q=-\frac{1}{2}(b+\mathrm{sign}(b)\sqrt{D})$, затем $r_1=q/a,r_2=c/q$ (обработать случаи $q=0$ и комплексные $\sqrt{D}$).
- Оценить чувствительность по формуле
$$|\Delta r|\le \frac{|r{|}^2|\Delta a|+|r||\Delta b|+|\Delta c|}{|2ar+b|}$$ и/или использовать $\kappa (r)\epsilon$ с $\kappa$ выше, где $\epsilon$ — машинный эпсилон или оценка относительной ошибки вычислений.
- При плохой обусловленности увеличить точность или применить интервальную арифметику.
Этого достаточно, чтобы надёжно получать корни для заданных интервалов коэффициентов и оценивать их погрешности.