Comments 22
Мать мои шаблоны!
А что у вас за проект такой, что приходится вычислять производные на этапе компиляции?
А что у вас за проект такой, что приходится вычислять производные на этапе компиляции?
Проекта никакого нет. Просто после прочтения статьи показалось интересной идея выведения точного выражения для производной с помощью шаблонов.
UFO just landed and posted this here
UFO just landed and posted this here
Операция численного дифференцирования не устойчивая, в отличии от численного интегрирования. Если простым языком, то при интегрировании гладкость функции возрастает, а при дифференцировании падает.
Поэтому предлагаю следующий трюк, который конечно не снимает проблемы, но позволяет лучше с ней справится в некоторых случаях. Вычислять производную через интеграл
int_{x_{i}}^{x_{i+1}} u_x dx = u(x_{i+1}) — u(x_{i})
Беря, например, для численного интегрирования метод трапеций получим разностное соотношение
(u_x(i+1) + u_x(i)) *(x_{i+1} — x_{i})/2 = u(x_{i+1}) — u(x_{i})
Которое дает n-2 линейных уравнений для i=2..n-1. Добавив к нему по условию на краях получим замкнутую систему, которую можно решить методом прогонки. Если будут проблемы при решении линейной системы, значит в данном случае нужно использовать какие-то другие методы и/или предположения.
PS. Тут u_x и u связаны, только через интегральные соотношения, т.е. u_x не является формально производной u по x. Она является дискретной функцией связанной с другой дискретной функции u определенной зависимостью. Обозначения выбраны только для удобства.
Поэтому предлагаю следующий трюк, который конечно не снимает проблемы, но позволяет лучше с ней справится в некоторых случаях. Вычислять производную через интеграл
int_{x_{i}}^{x_{i+1}} u_x dx = u(x_{i+1}) — u(x_{i})
Беря, например, для численного интегрирования метод трапеций получим разностное соотношение
(u_x(i+1) + u_x(i)) *(x_{i+1} — x_{i})/2 = u(x_{i+1}) — u(x_{i})
Которое дает n-2 линейных уравнений для i=2..n-1. Добавив к нему по условию на краях получим замкнутую систему, которую можно решить методом прогонки. Если будут проблемы при решении линейной системы, значит в данном случае нужно использовать какие-то другие методы и/или предположения.
PS. Тут u_x и u связаны, только через интегральные соотношения, т.е. u_x не является формально производной u по x. Она является дискретной функцией связанной с другой дискретной функции u определенной зависимостью. Обозначения выбраны только для удобства.
Вообще-то функция derivative выводит аналитическое выражение для производной. Приближенное вычисление используется только в случае дифференцирования произвольного функтора или лямбда-выражения.
1. Посчитайте производную от abs(x) при x=0. Поскольку у Вас для нее нет шаблона, значит будет считаться численно.
2. Вычисление значение многочлена в точке для операций с плавающей точкой тоже неустойчиво. Хотя производная легко считается в виде аналитического выражения, подсчет значения производной для конкретного x может уже при степенях > 20 может отличаться на порядки.
3. Я не увидел у Вас системы упрощения выражений. Т.е. могут возникать очень плохие выражения для операций с плавающей точки типа 0 — 0. Нарушение ассоциативности операций сложения (или вычитания) за счет округления. Часто
(a + b) + c != a + (b + c)
поэтому может помочь «правильное» переписывание формул по которым считаете. С умножением все тип-топ. Я понимаю что у всех в башке матанализ, а должна быть алгебра. Но это исторический дефект нашего математического образования.
Читайте Кнута «Искусство программирования» т.2 «Получисленные алгоритмы».
4. У Вас все сделано не плохо, но если делать правильно математически, то получится небольшая система компьютерной алгебры. Я помню с каким-то прологом поставился пример аналитического дифференцирования функций. Буквально несколько строчек, но потом огромный код просто по упрощению полученного выражения. Он при этом часто работал некорректно. Поскольку есть целый набор теорем о том, что для некоторых типов выражений нет канонического представления и даже не существует алгоритм нормальной формы.
5. Поэтому не смотря на все недостатки операций с плавающей точкой, все таки и проще и лучше для C++ выбрать численный вариант.
2. Вычисление значение многочлена в точке для операций с плавающей точкой тоже неустойчиво. Хотя производная легко считается в виде аналитического выражения, подсчет значения производной для конкретного x может уже при степенях > 20 может отличаться на порядки.
3. Я не увидел у Вас системы упрощения выражений. Т.е. могут возникать очень плохие выражения для операций с плавающей точки типа 0 — 0. Нарушение ассоциативности операций сложения (или вычитания) за счет округления. Часто
(a + b) + c != a + (b + c)
поэтому может помочь «правильное» переписывание формул по которым считаете. С умножением все тип-топ. Я понимаю что у всех в башке матанализ, а должна быть алгебра. Но это исторический дефект нашего математического образования.
Читайте Кнута «Искусство программирования» т.2 «Получисленные алгоритмы».
4. У Вас все сделано не плохо, но если делать правильно математически, то получится небольшая система компьютерной алгебры. Я помню с каким-то прологом поставился пример аналитического дифференцирования функций. Буквально несколько строчек, но потом огромный код просто по упрощению полученного выражения. Он при этом часто работал некорректно. Поскольку есть целый набор теорем о том, что для некоторых типов выражений нет канонического представления и даже не существует алгоритм нормальной формы.
5. Поэтому не смотря на все недостатки операций с плавающей точкой, все таки и проще и лучше для C++ выбрать численный вариант.
Цель поста — показать, что с помощью шаблонов на С++ в принципе можно обрабатывать несложные математические формулы. Например, дифференцировать. Я ни в коем случае не собирался делать собственную систему компьютерной алгебры. Вы можете обратить внимание на название библиотеки — CrazyMath :) Разве это подходящее имя для системы компьютерной алгебры?
Оптимизации простенькие я делал, но они были направлены против практически экспоненциального роста сложности вычисления производных высших порядков. Уже при вычислении производной 3-го порядка с помощью неоптимизированного кода мой компилятор (VC++ 2010) крепко задумывался. Упрощение выражений — безусловно очень сложная задача. Но я считаю, что ее следует рассматривать отдельно. Можно попытаться написать что-то вроде template class Simplify;
Оптимизации простенькие я делал, но они были направлены против практически экспоненциального роста сложности вычисления производных высших порядков. Уже при вычислении производной 3-го порядка с помощью неоптимизированного кода мой компилятор (VC++ 2010) крепко задумывался. Упрощение выражений — безусловно очень сложная задача. Но я считаю, что ее следует рассматривать отдельно. Можно попытаться написать что-то вроде template class Simplify;
Неужели ни автор статьи, ни читатели не видят, что вычисления, приведенные выше, выполняются в run-time, а не в compile-time? Что все эти функторы, описанные автором, вызываются в run-time? Автор не имеет базового понятия о шаблонном мета-программировании.
Непосредственно функтор, соответствующий аналитическому выражению производной, вычисляется в компайл-тайме. А вот в этом коде в компайл-тайме будут вычислены и конкретные цифры:
Конечно, для этого придется добавить модификатор const к оператору «скобочки», а также к методу expression(). Да, это я не учел. Но модифицированный код работает и вычисляет производную на этапе компиляции.
const auto func = Sqr(X) * Sqrt(X);
const auto d_func = derivative(func);
const double d_func_value = d_func(1);
Конечно, для этого придется добавить модификатор const к оператору «скобочки», а также к методу expression(). Да, это я не учел. Но модифицированный код работает и вычисляет производную на этапе компиляции.
Всем, кто пытается найти производную алгебраическими методами, а уж тем более методами ООП, следует попробовать найти производную от функции (3x+1)/(2x+1) в точке 10^20.
Верное значение — 0,25 * 10^-40. Между прочим, это значительно больше, чем double.MIN_VALUE.
А если получился ноль, как у автора статьи — есть повод задуматься.
PS а если эту производную вы успешно вычислили, то замените один из иксов в формуле на 1/sin(1/x)) и попробуйте снова.
Верное значение — 0,25 * 10^-40. Между прочим, это значительно больше, чем double.MIN_VALUE.
А если получился ноль, как у автора статьи — есть повод задуматься.
PS а если эту производную вы успешно вычислили, то замените один из иксов в формуле на 1/sin(1/x)) и попробуйте снова.
>> Верное значение — 0,25 * 10^-40. Между прочим, это значительно больше, чем double.MIN_VALUE.
Да, но 10^20 в квадрате уже больше, чем может поместиться в double. Вы прекрасно понимаете, что как раз в этом примере точности double недостаточно для вычисления «в лоб» без предварительного преобразования выражений. Функтор, полученный с помощью моего кода, соответствует корректному математическому выражению для вычисления производной. То, что его надо упростить, а затем преобразовать для снижения погрешностей вычисления — это уже отдельная задача. Задача далеко не тривиальная и может послужить темой целой книги, а не поста на Хабре. Мой же пост был посвящен одному из вариантов применения шаблонов С++.
Да, но 10^20 в квадрате уже больше, чем может поместиться в double. Вы прекрасно понимаете, что как раз в этом примере точности double недостаточно для вычисления «в лоб» без предварительного преобразования выражений. Функтор, полученный с помощью моего кода, соответствует корректному математическому выражению для вычисления производной. То, что его надо упростить, а затем преобразовать для снижения погрешностей вычисления — это уже отдельная задача. Задача далеко не тривиальная и может послужить темой целой книги, а не поста на Хабре. Мой же пост был посвящен одному из вариантов применения шаблонов С++.
Вы ошибаетесь!
10^20 в квадрате есть 10^40, что помещается в double без проблем (максимальное значение — более 10^300).
В том-то и дело, что алгоритм упрощения выражения — нетривиален. И никто и никогда не реализует его на шаблонах (разве что кто-нибудь придумает транслятор с нормального языка в шаблоны). А потому, сферой применимости шаблонов C++ дифференцирование не является!
10^20 в квадрате есть 10^40, что помещается в double без проблем (максимальное значение — более 10^300).
В том-то и дело, что алгоритм упрощения выражения — нетривиален. И никто и никогда не реализует его на шаблонах (разве что кто-нибудь придумает транслятор с нормального языка в шаблоны). А потому, сферой применимости шаблонов C++ дифференцирование не является!
>> Вы ошибаетесь!
Ну да, я ошибся. Не стоит это драматизировать. По-моему вы употребляете восклицательные знаки чаще, чем нужно.
>> 10^20 в квадрате есть 10^40, что помещается в double без проблем (максимальное значение — более 10^300).
Дело тут в сложении чисел разных порядков. Точности мантиссы для этого просто не хватает. Выражение для вычисления производной, а еще лучше саму функцию, нужно предварительно преобразовать с целью минимизации количества операций сложения чисел разных порядков.
Не надо меня убеждать в том, что шаблоны С++ неприменимы для дифференцирования. И уж тем более тот факт, что «никто и никогда не реализует его на шаблонах», не является для меня аргументом. Из того, что алгоритм нетривиален, не следует то, что алгоритм нереализуем.
Ну да, я ошибся. Не стоит это драматизировать. По-моему вы употребляете восклицательные знаки чаще, чем нужно.
>> 10^20 в квадрате есть 10^40, что помещается в double без проблем (максимальное значение — более 10^300).
Дело тут в сложении чисел разных порядков. Точности мантиссы для этого просто не хватает. Выражение для вычисления производной, а еще лучше саму функцию, нужно предварительно преобразовать с целью минимизации количества операций сложения чисел разных порядков.
Не надо меня убеждать в том, что шаблоны С++ неприменимы для дифференцирования. И уж тем более тот факт, что «никто и никогда не реализует его на шаблонах», не является для меня аргументом. Из того, что алгоритм нетривиален, не следует то, что алгоритм нереализуем.
Sign up to leave a comment.
Вычисление производных с помощью шаблонов на С++