Хабр Курсы для всех
РЕКЛАМА
Практикум, Хекслет, SkyPro, авторские курсы — собрали всех и попросили скидки. Осталось выбрать!
Хабр доступен 24/7 благодаря поддержке друзей
Комментарии 36
Сайт, посвящённый проблеме точности вычислений в разных языках: http://0.30000000000000004.com/
Делай раз — смотрим на формулу: у вас 3*10^4 возводится в степень 8.
Делай два — результат этого возведения — 3*10^32 — упс, вышли за приделы точности double
Делай три — открываем тетрадочку с лекциями по численным методам и читаем «Для обеспечения устойчивого счета, требуется заложить минимум 4 порядка точности арифметики». Суммируем 32+4=36. Ой, что-то теория работает.
Это я к чему. Вычислительная математика работает строго по принципу GIGO — мусор на входе — мусор на выходе. Вы заранее засунули компьютеру ерунду, которую он обсчитать не может, и удивляетесь, что он выдал вам ерунду. Так математический софт не пишут.
Делай два — результат этого возведения — 3*10^32 — упс, вышли за приделы точности double
Делай три — открываем тетрадочку с лекциями по численным методам и читаем «Для обеспечения устойчивого счета, требуется заложить минимум 4 порядка точности арифметики». Суммируем 32+4=36. Ой, что-то теория работает.
Это я к чему. Вычислительная математика работает строго по принципу GIGO — мусор на входе — мусор на выходе. Вы заранее засунули компьютеру ерунду, которую он обсчитать не может, и удивляетесь, что он выдал вам ерунду. Так математический софт не пишут.
Плюсанусь, однако замечу, что не столь важен максимальный порядок, сколь разница между максимальным и минимальным. Собственно, вангую, что не будь «а/2бе», который таки порядка 10^0, то и устойчивость достигалась бы на меньшей точности, чем 36.
«Устойчивость» вообще термин здесь малоуместный.
«Устойчивость» означает, что малые изменения параметров системы (обычно — коэффициентов диффура, или граничных условий), ведут к малым изменениям решения.
И этот термин относится к алгоритму счета, там и критерий специальный есть, и теория, которая этот критерий строго доказывает, и так далее.
Здесь же у нас чисто арифметическая задача «сколько брать знаков на расчет» и ответ к ней — давным давно эмпирически предложен — запас должен быть 4 (это как самый минимум!) знака.
Вы совершенно правы, указывая на то, что нужно брать разницу между самой большой и самой меньшей степенями десятки.
Однако, при численном решении реальных задач (в матфизике, например), задачу сначала нормируют (чтобы самое большое число по модулю было равно 1), а затем уже с нормированной и безразмерной задачей работают. Тогда критерий становится простым — «хотим N знаков точности — берем 4+N знака арифметики (а лучше — больше)».
«Устойчивость» означает, что малые изменения параметров системы (обычно — коэффициентов диффура, или граничных условий), ведут к малым изменениям решения.
И этот термин относится к алгоритму счета, там и критерий специальный есть, и теория, которая этот критерий строго доказывает, и так далее.
Здесь же у нас чисто арифметическая задача «сколько брать знаков на расчет» и ответ к ней — давным давно эмпирически предложен — запас должен быть 4 (это как самый минимум!) знака.
Вы совершенно правы, указывая на то, что нужно брать разницу между самой большой и самой меньшей степенями десятки.
Однако, при численном решении реальных задач (в матфизике, например), задачу сначала нормируют (чтобы самое большое число по модулю было равно 1), а затем уже с нормированной и безразмерной задачей работают. Тогда критерий становится простым — «хотим N знаков точности — берем 4+N знака арифметики (а лучше — больше)».
Прежде всего, f(а, b) — отображение, заданное в виде арифметического выражения, относительно простого, да. Но оно очень чувствительно к изменениям входящих параметров; это — фактически — явление неустойчивости.
Ваша формула — полином. Это гладкая функция, непрерывная вместе со всеми своими производными. Изменение параметров на бесконечно малую величину в данном случае приведет к изменению результата также на бесконечно малую величину.
В очередной раз повторяю, что термин «устойчивость» относится к конкретному алгоритму вычислений, а не к формуле. Вы же не приводите алгоритм (как конкретно интерпретатор пайтона это все считал? Прямо с точностью до действий?).
Например, для задачи об охлаждении неоднородного стержня, алгоритм прямого счета является условно устойчивым (нужно придерживаться определенного соотношения шагов разностной схемы, иначе на выходе — мусор), а схема Кранка-Николсон (не НиколсонА, Филлис Николсон — женщина!) — абсолютно устойчива и работает для любого соотношения шагов.
В очередной раз повторяю, что термин «устойчивость» относится к конкретному алгоритму вычислений, а не к формуле. Вы же не приводите алгоритм (как конкретно интерпретатор пайтона это все считал? Прямо с точностью до действий?).
Например, для задачи об охлаждении неоднородного стержня, алгоритм прямого счета является условно устойчивым (нужно придерживаться определенного соотношения шагов разностной схемы, иначе на выходе — мусор), а схема Кранка-Николсон (не НиколсонА, Филлис Николсон — женщина!) — абсолютно устойчива и работает для любого соотношения шагов.
Вычислительная математика работает и по другим принципам, иногде получая вполне корректные данные на входе, но выдавая мусор на выходе, и дело часто отнюдь не только в том, что где-то мы вылазим за пределы допустимых значений или теряем разряды. Проблема в том, что некоторые вычисления являются вычислительно нестабильными, например элементарное вычитание двух чисел либо суммирование ряда.
Я не буду сейчас приводить и объяснять конкретные примеры, но в одном из них выполнение серии безобидных операций без выходов за пределы чисел на выходе дает число, в котором все разряды кроме старшего неправильные. Есть и контр-примеры, в которых изменение порядка вычислений сохраняет вычислительную стабильность даже при переполнениях.
Я не буду сейчас приводить и объяснять конкретные примеры, но в одном из них выполнение серии безобидных операций без выходов за пределы чисел на выходе дает число, в котором все разряды кроме старшего неправильные. Есть и контр-примеры, в которых изменение порядка вычислений сохраняет вычислительную стабильность даже при переполнениях.
В данном случае, под входным мусором понимается как алгоритм, так и данные, которые ему скормили. Так что GIGO и здесь продолжает работать.
вылазим за пределы допустимых значений или теряем разряды
Есть ещё ошибки округления. Но это вроде бы полный набор. Какие ещё нестабильности вы знаете?
В вычислительной математике выделяют четыре категории ошибок:
-ошибки входных данных (неточности измерений физических величин)
-ошибки, вносимые математической моделью (несоответствие между моделью и реальностью)
-погрешность численного метода (неустойчивость из-за «плохой» геометрии, погрешность аппроксимации численной схемы и так далее)
-погрешность арифметики
-ошибки, вносимые математической моделью (несоответствие между моделью и реальностью)
-погрешность численного метода (неустойчивость из-за «плохой» геометрии, погрешность аппроксимации численной схемы и так далее)
-погрешность арифметики
Не совсем правильно считать данный пример ерундой; Для иллюстрации влияния погрешностей округления, действительно, можно привести пример гораздо проще, но смысл данного примера в другом: допустим мы решили посчитать данное выражение с одинарной точностью, посчитали; мы понимаем, что здесь степени и возможны влияния погрешностей округления. Для проверки результата мы можем пересчитать все это в удвоенной точности, и получим такой же результат! Вычисляя выражение с учетверенной точностью (все в рамках IEEE 754), снова получаем тот же результат. Так легко и поверить, что он правильный, но он-то оказывается не верным!
Вывод из примера, что получение мало отличающихся результатов при увеличении точности не гарантирует, что результат правильный.
Вывод из примера, что получение мало отличающихся результатов при увеличении точности не гарантирует, что результат правильный.
Прежде чем что-либо считать, всегда можно априорно оценить результат.
Более того, всегда можно заставить компьютер работать с символьными выражениями, то есть сохранять все знаки (симулируя старательного школьника, который решает пример).
В очередной раз повторяю — никакие реальные вычисления не делаются без априорных теоретических оценок. Доверять можно только сочетанию из формального доказательства и результата работы программы. Иначе — получите GIGO. В вычислительной практике никто не опирается только на то, что насчитал компьютер.
Более того, всегда можно заставить компьютер работать с символьными выражениями, то есть сохранять все знаки (симулируя старательного школьника, который решает пример).
В очередной раз повторяю — никакие реальные вычисления не делаются без априорных теоретических оценок. Доверять можно только сочетанию из формального доказательства и результата работы программы. Иначе — получите GIGO. В вычислительной практике никто не опирается только на то, что насчитал компьютер.
Wolfram Alpha считает что ответ примерно равен -4.27229×10^20
А мне Wolfram Alpha сказал, что -54767/66192. Это как раз примерно -0.82739.
Я вставлял ссылку WA, но ее почемуто удалили ;( не знаю почему у нас разные ответы
Я вставлял в форму WA вот такую строку: «a = 77617, b = 33096, c = (333.75 — a * a) * b ** 6 + a * a * (11 * a * a * b * b — 121 * b ** 4 — 2) + 5.5 * b ** 8 + a/(2.0*b)». Возможно это не совсем правильный способ, но по другому я не знаю как. Ответ получил как писал выше.
Еще я пытался вводить строку с уже подставленными значениями a и b: "(333.75 — 77617 * 77617) * 33096 ** 6 + 77617 * 77617 * (11 * 77617 * 77617 * 33096 * 33096 — 121 * 33096 ** 4 — 2) + 5.5 * 33096 ** 8 + 77617/(2.0*33096)". В этом случае WA ответил без дробей, сразу числом из поста.
Еще я пытался вводить строку с уже подставленными значениями a и b: "(333.75 — 77617 * 77617) * 33096 ** 6 + 77617 * 77617 * (11 * 77617 * 77617 * 33096 * 33096 — 121 * 33096 ** 4 — 2) + 5.5 * 33096 ** 8 + 77617/(2.0*33096)". В этом случае WA ответил без дробей, сразу числом из поста.
Мне такое сказал Maple — но не сразу. Пришлось представить 333.75 как 333+3/4, и т.п…
Вообще пример хороший — например в прикладном софте такое запросто можно проморгать. Формула есть? Есть. Выглядит вроде адекватно — никаких мелких знаменателей нет. Восьмая степень — это все ж не двадцать восьмая. Уравнение состояния Бенедикта — Вебба — Рубина содержит шестую степень — и ничего. Так что вопросы откуда взялась такая странная формула возникнут не у всех, а значит «тестирования» формулы может и не быть.
Вообще пример хороший — например в прикладном софте такое запросто можно проморгать. Формула есть? Есть. Выглядит вроде адекватно — никаких мелких знаменателей нет. Восьмая степень — это все ж не двадцать восьмая. Уравнение состояния Бенедикта — Вебба — Рубина содержит шестую степень — и ничего. Так что вопросы откуда взялась такая странная формула возникнут не у всех, а значит «тестирования» формулы может и не быть.
Я не работал с WA, но он может и «схитрить» и все символьно посчитать, даже, пожалуй, подменив 333.75=333+75/100.
Ну прям секрет открыли. А еще unsigned int нельзя использовать для вычислений с результатом, который выходит за пределы 2^32.
Если посчитать первый множитель в обычном калькуляторе, то получается результат -7,9171109033773850490791882372801e+36, откуда видно, что для точных вычислений надо как минимум 37 разрядов.
Если посчитать первый множитель в обычном калькуляторе, то получается результат -7,9171109033773850490791882372801e+36, откуда видно, что для точных вычислений надо как минимум 37 разрядов.
Плюс-минус 1 — это офигительная точность для того диапазона значений, которые принимает данная функция))
Мы ж не деньги считаем (их в 8ю степень не возводят).
Мы ж не деньги считаем (их в 8ю степень не возводят).
Благодарю всех за комментарии, но, к сожалению, практически все здесь увидели пример, лишь как обычную иллюстрацию влияния погрешностей округления. Смысл примера — в другом. Мы получаем практически одинаковые неправильные результаты вычислений последовательно повышая точность представления чисел, и это интересно и показательно, т.к. получение мало отличающихся результатов на различных разрядных сетках является «доказательством» его правильности.
Я бы сказал, это смысл комментов в другом. То, что вы повысили точность, не означает, что она стала достаточной. Не знаю, почему вы считаете, что разные разрядные сетки что-то доказывают. Нельзя получить 1 в выражении 1000000000877 — 876 — 1000000000000, если разрядная сетка 1, 2, 4, или 8 знаков, и младшие разряды в первом числе превращаются в 0. Смысл не в количестве разрядов, а в том, помещается ли туда результат.
Поэтому фраза
Поэтому фраза
Но откуда знать, что при дальнейшем увеличении точности не произойдет еще какого-либо скачкатехнически абсолютно неграмотна. Оттуда и знать, что точность теперь достаточная. Почему достаточная? Потому что порядок множителей меньше или равен 10^36.
Не знаю, почему вы считаете, что разные разрядные сетки что-то доказывают...
Конечно не доказывают, но могут придавать уверенности, особенно в сложных вычислениях. (реально сталкивался с тем, что люди выполняют «такую проверку», осуществляя расчеты в академических учреждениях...)
Если бы наша задача состояла в вычислении собственных чисел какой-либо несимметричной матрицы, процедура вычислений не была бы столь прозрачна, и было бы очень замачиво «проверить» результат, вычислив его, например, на удвоенной и расширенной удвоенной точностях. И вот, мы проделали вычисления, и результаты «почти» совпали; вдобавок еще провели вычисления на учетвернной точности, получили «такие же» значения… так недолго и поверить, что результат «правильный» (понимаю, что вычисления на различных разрядных сетках — не основание для проверки результатов вычислений собственных значений несим. матриц, для этого есть другие подходы — например, построение их спектральных портретов — псевдоспектров).
Пример показывает, что «проверка» на различных разрядных сетках ничего не значит, т.к. почти одинаковые неправильные значения могут получаться и при увеличении точности.
Про наличие еще «одного скачка» — применительно к примеру, согласен, что фраза не учитывает, что в этом простом выражении можно легко оценить необходимое число разрядов; а если бы процедура вычислений была бы достаточно сложна и получение оценки требуемого числа разрядов было бы проблематично…
Если бы наша задача состояла в вычислении собственных чисел какой-либо несимметричной матрицы, процедура вычислений не была бы столь прозрачна, и было бы очень замачиво «проверить» результат, вычислив его, например, на удвоенной и расширенной удвоенной точностях.
Совершенно бессмысленная затея. В ситуации с матрицами и собственными числами нужно оценивать число обусловленности матрицы, а не гонять компьютер почем зря. Более того, для сколько-нибудь «интересной» задачи расчет без аппаратной поддержки арифметики займет чрезвычайно длительное время, а на выходе даст пшик.
Число обусловленности здесь не особо подходит, например, для вырожденой матрицы число обусловленности (спектральное) есть бесконечность, и что делать дальше, если собственные значения нужно найти (одно из этих значений будет 0, а остальные)… Но здесь не об этом…
В MatLab есть тестовая матрица gallery(5), следующего вида:
А = [[-9, 11, -21, 63, -252], [70, -69, 141, -421, 1684], [-575, 575, -1149, 3451, -13801], [3891, -3891, 7782, -23345, 93365], [1024, -1024, 2048, -6144, 24572]]. Нужно найти ее собственные значения (несимвольно);
Есть даже ее совсем неправильный спетральный портрет, который подталкивает поверить в неправильно вычисленные собственные значения.
Вот и неправильный результат: eigs (A) = ( -0.03697710+0.02749668j, -0.03697710-0.02749668j, 0.01470225+0.04265918j, 0.01470225-0.04265918j, 0.04454971+0.j ) — (вычислено при помощи numpy.linalg.eig); MatLab раньше (R2008a, например) давал также неправильный результат, сейчас — не знаю, т.к. не пользуюсь MatLab давно.
Какой же результат правильный?! — 0.0 (т.е. 5-ти кратное нулевое соб. значение)
Правда этот пример не такой красивый (для игр с точностями), т.к. на различных точностях дает различные наборы собственных значений.
В MatLab есть тестовая матрица gallery(5), следующего вида:
А = [[-9, 11, -21, 63, -252], [70, -69, 141, -421, 1684], [-575, 575, -1149, 3451, -13801], [3891, -3891, 7782, -23345, 93365], [1024, -1024, 2048, -6144, 24572]]. Нужно найти ее собственные значения (несимвольно);
Есть даже ее совсем неправильный спетральный портрет, который подталкивает поверить в неправильно вычисленные собственные значения.
Вот и неправильный результат: eigs (A) = ( -0.03697710+0.02749668j, -0.03697710-0.02749668j, 0.01470225+0.04265918j, 0.01470225-0.04265918j, 0.04454971+0.j ) — (вычислено при помощи numpy.linalg.eig); MatLab раньше (R2008a, например) давал также неправильный результат, сейчас — не знаю, т.к. не пользуюсь MatLab давно.
Какой же результат правильный?! — 0.0 (т.е. 5-ти кратное нулевое соб. значение)
Правда этот пример не такой красивый (для игр с точностями), т.к. на различных точностях дает различные наборы собственных значений.
Будем последовательны. Число обусловленности прекрасно «работает» в случае вырожденной матрицы — оно сообщает нам, что матрица поганая, а следовательно, стандартный метод для нее — не подходит.
В вычислительной математик1е нет одного универсального метода, который бы оптимально подходил для всех задач, к каждому методу жестко прибит гвоздями-сотками критерий его применимости (как и вообще в науке, кстати). Вы в данной публикации критерии применимости систематически игнорируете, в результате закономерно получаете лажу.
В вычислительной математик1е нет одного универсального метода, который бы оптимально подходил для всех задач, к каждому методу жестко прибит гвоздями-сотками критерий его применимости (как и вообще в науке, кстати). Вы в данной публикации критерии применимости систематически игнорируете, в результате закономерно получаете лажу.
Всё правильно вы сделали. В современном мире люди без спец математического образования могут, что-то считать по формулам и не знать ничего про коварные точности.
А все ваши критики в этой статье просто матемаческие снобы, считающие людей с не таким образованием просто более низшей кастой.
А все ваши критики в этой статье просто матемаческие снобы, считающие людей с не таким образованием просто более низшей кастой.
Мне кажется, вы переусложняете без необходимости.
У каждого программиста есть более простой пример перед глазами:
И что это доказывает? Что 2 ** 36 равно нулю, поскольку все ваши вычисления дали нуль?
Так и здесь: поскольку вы отнимаете друг от друга числа порядка 10**36, то для точного результата мантисса должна вмещать в себя по крайней мере 36 десятичных разрядов. Сколько разрядов вмещает мантисса у double, помните? Всего лишь 15-17.
Это даже в википедии написано: English Wikipedia: Double-precision_floating-point_format
У каждого программиста есть более простой пример перед глазами:
(int8) (2 ** 36) совпадает с(int16) (2 ** 36) и с(int32) (2 ** 36).И что это доказывает? Что 2 ** 36 равно нулю, поскольку все ваши вычисления дали нуль?
Так и здесь: поскольку вы отнимаете друг от друга числа порядка 10**36, то для точного результата мантисса должна вмещать в себя по крайней мере 36 десятичных разрядов. Сколько разрядов вмещает мантисса у double, помните? Всего лишь 15-17.
Это даже в википедии написано: English Wikipedia: Double-precision_floating-point_format
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Вычисления с плавающей точкой: можно ли доверять результатам?