При вычислении корней квадратного уравнения ax^2 + bx + c = 0 в случае, когда |b| >> |a| и |b| >> |c|, прямая формула может давать потерю значащих цифр; проанализируйте причину численной нестабильности и предложите устойчивую на практике схему вычисления корней, объяснив её преимущества.
Коротко — причина и устойчивый приём. Причина численной нестабильности. - Стандартная формула x1,2=−b±b2−4ac2a \;x_{1,2}=\dfrac{-b\pm\sqrt{b^2-4ac}}{2a}\;x1,2=2a−b±b2−4ac содержит разность двух больших близких по модулю чисел −b-b−b и b2−4ac\sqrt{b^2-4ac}b2−4ac при ∣b∣≫∣a∣,∣c∣|b|\gg|a|,|c|∣b∣≫∣a∣,∣c∣. Это приводит к катастрофическому округлению (потере значащих цифр) при вычислении меньшей по модулю части числителя. - Приближённо b2−4ac≈∣b∣(1−2acb2)\sqrt{b^2-4ac}\approx|b|\bigl(1-\tfrac{2ac}{b^2}\bigr)b2−4ac≈∣b∣(1−b22ac), поэтому например для b>0b>0b>0 числитель −b+D≈−2ac/b-b+\sqrt{D}\approx -2ac/b−b+D≈−2ac/b — маленькое число, получаемое как разность двух величин порядка ∣b∣|b|∣b∣, что вычёркивает значащие цифры. Устойчивая на практике схема (стандартный надёжный приём). 1. Вычислить дискриминант D=b2−4ac \;D=b^2-4ac\;D=b2−4ac и его корень s=D \;s=\sqrt{D}\;s=D (в комплексной арифметике — тоже работает). 2. Задать q=−12(b+sign(b) s),
q=-\tfrac{1}{2}\bigl(b+\operatorname{sign}(b)\,s\bigr), q=−21(b+sign(b)s),
где sign(b)=1\operatorname{sign}(b)=1sign(b)=1 при b≥0b\ge0b≥0, иначе −1-1−1. 3. Тогда устойчиво вычислить корни как x1=qa,x2=cq.
x_1=\dfrac{q}{a},\qquad x_2=\dfrac{c}{q}. x1=aq,x2=qc.
(Если нужно — можно поменять индексы; при q=0q=0q=0 вернуть двойной корень −b/(2a)-b/(2a)−b/(2a).) Почему это устойчиво. - Выражение b+sign(b) sb+\operatorname{sign}(b)\,sb+sign(b)s складывает числа одного знака (вместо их вычитания), так что величина в qqq не образуется как малое остаточное значение двух больших близких чисел — поэтому не теряются значащие цифры. - Второй корень получается через отношение c/qc/qc/q, что использует соотношение произведения корней x1x2=c/ax_1x_2=c/ax1x2=c/a и не вводит новых вычитаний больших близких величин. - Работает и для комплексных корней (с корректной реализацией D\sqrt{D}D и правилом выбора знака). Краткое замечание об особых случаях. - Если D=0D=0D=0: корень −b/(2a)-b/(2a)−b/(2a). - Если a=0a=0a=0: это не квадратичное уравнение — решать как линейное. - Если qqq очень близко к нулю по машинной точности, используйте симметричную схему или повышенную точность.
Причина численной нестабильности.
- Стандартная формула x1,2=−b±b2−4ac2a \;x_{1,2}=\dfrac{-b\pm\sqrt{b^2-4ac}}{2a}\;x1,2 =2a−b±b2−4ac содержит разность двух больших близких по модулю чисел −b-b−b и b2−4ac\sqrt{b^2-4ac}b2−4ac при ∣b∣≫∣a∣,∣c∣|b|\gg|a|,|c|∣b∣≫∣a∣,∣c∣. Это приводит к катастрофическому округлению (потере значащих цифр) при вычислении меньшей по модулю части числителя.
- Приближённо b2−4ac≈∣b∣(1−2acb2)\sqrt{b^2-4ac}\approx|b|\bigl(1-\tfrac{2ac}{b^2}\bigr)b2−4ac ≈∣b∣(1−b22ac ), поэтому например для b>0b>0b>0 числитель −b+D≈−2ac/b-b+\sqrt{D}\approx -2ac/b−b+D ≈−2ac/b — маленькое число, получаемое как разность двух величин порядка ∣b∣|b|∣b∣, что вычёркивает значащие цифры.
Устойчивая на практике схема (стандартный надёжный приём).
1. Вычислить дискриминант D=b2−4ac \;D=b^2-4ac\;D=b2−4ac и его корень s=D \;s=\sqrt{D}\;s=D (в комплексной арифметике — тоже работает).
2. Задать
q=−12(b+sign(b) s), q=-\tfrac{1}{2}\bigl(b+\operatorname{sign}(b)\,s\bigr),
q=−21 (b+sign(b)s), где sign(b)=1\operatorname{sign}(b)=1sign(b)=1 при b≥0b\ge0b≥0, иначе −1-1−1.
3. Тогда устойчиво вычислить корни как
x1=qa,x2=cq. x_1=\dfrac{q}{a},\qquad x_2=\dfrac{c}{q}.
x1 =aq ,x2 =qc . (Если нужно — можно поменять индексы; при q=0q=0q=0 вернуть двойной корень −b/(2a)-b/(2a)−b/(2a).)
Почему это устойчиво.
- Выражение b+sign(b) sb+\operatorname{sign}(b)\,sb+sign(b)s складывает числа одного знака (вместо их вычитания), так что величина в qqq не образуется как малое остаточное значение двух больших близких чисел — поэтому не теряются значащие цифры.
- Второй корень получается через отношение c/qc/qc/q, что использует соотношение произведения корней x1x2=c/ax_1x_2=c/ax1 x2 =c/a и не вводит новых вычитаний больших близких величин.
- Работает и для комплексных корней (с корректной реализацией D\sqrt{D}D и правилом выбора знака).
Краткое замечание об особых случаях.
- Если D=0D=0D=0: корень −b/(2a)-b/(2a)−b/(2a).
- Если a=0a=0a=0: это не квадратичное уравнение — решать как линейное.
- Если qqq очень близко к нулю по машинной точности, используйте симметричную схему или повышенную точность.