Как избежать катастрофической отмены для небольших чисел в f(x) = (1-cos(x))/x**2 в Python 3.7?

Как избежать катастрофической отмены для небольших чисел в f(x) = (1-cos(x))/x**2 в Python 3.7?

Это то, что я пытался до сих пор (ключ, я знаю, это некоторая тригонометрическая идентичность, которая позволяет вам избежать сокращения, и я также знаю, используя правило Лопиталя, что предел → 0 для f (x) равен 0,5 , поэтому правильный вывод программы очень близок к 0,5, что вы и получите, если, например, используете x = 1,2e-4, но вы получите отмену с меньшими числами, такими как 1,2e-8, и мне нужно сделать это так не бывает).

from math import *
def f(x):     #these are all the same function using different identities  
   a = (1-(sin(x)/tan(x)))/(x**2)
   b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
   c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
   d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
   e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
   return a, b, c, d, e

print(k(1.2e-8))
#Output: (0.0, 0.7709882115452477, 0.0, 0.0, 0.0) - whereas need 0.5000...

python-3.x python function trigonometry floating-accuracy
person Corwin of Amber    schedule 21.02.2021    source источник
comment
вы могли бы рассмотреть возможность использования numpy с 64-битным float64. Это может помочь с отменой.   -  person Tim    schedule 21.02.2021
comment
Есть еще одно тождество: sin(x)/tan(x) = cos(x)   -  person Stefan    schedule 21.02.2021

Ответы (4)


arrow_upward
4
arrow_downward

Как это:

sin(x)/x * tan(x/2)/x

Он выполняет свою работу до конца, x = 1e-308 все еще в порядке.

К сожалению, я не могу подробно рассказать, почему это работает хорошо.

person harold    schedule 21.02.2021
comment
Вероятно, это связано с реализацией sin против cosine. Обратите внимание, что sin(1e-100)=1e-100, а cos(1e-100)=1.0. Точно так же tan(1e-100)=1e-100. - person Tim; 21.02.2021
comment
Я предполагаю, что они использовали приближение sin(x)=x, когда x мало. Вы можете показать с помощью L'Hopital, что sin(x)/x приближается к 1, когда x приближается к 0. - person Tim; 21.02.2021
comment
Вы можете заменить sin(x)/x на numpy.sinc(x/numpy.pi) - sinc имеет предел 1, когда x приближается к 0. Tan(x)/x также можно выразить с помощью sinc - person Stefan; 21.02.2021
comment
Идеальный! Это потрясающе, большое спасибо @harold, Тиму и Стефану. Ваш вклад очень полезен и ценится. - person Corwin of Amber; 21.02.2021
comment
Используйте Херби! herbie.uwplse.aorg/demo - person soegaard; 22.02.2021

arrow_upward
2
arrow_downward

Вместо этого используйте numpy.float128. numpy — это стандарт для анализа данных и более сложной математики. Его можно установить с помощью следующей команды в терминале.

pip install numpy

from numpy import *
def f(x):     #these are all the same function using different identities  
   a = (1-(sin(x)/tan(x)))/(x**2)
   b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
   c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
   d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
   e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
   return a, b, c, d, e
print(f(float128(1.2e-8)))

Это печатает

(0.5003141275115400736, 0.49956120933620291774, 0.49993766842387149567, 0.49993766842387149567, 0.49956120933620291774)
person Tim    schedule 21.02.2021
comment
float128 эквивалентен float64 на некоторых платформах. - person harold; 21.02.2021
comment
изменил numpy.float128 на float128, так как мы все еще хотим использовать numpy версии триггерных функций. - person Tim; 21.02.2021
comment
float64, к сожалению, не работает, поэтому, возможно, придется отложить ответ @harold, если ваша платформа по умолчанию использует float128 для float64. - person Tim; 21.02.2021
comment
Большое спасибо за ваш ответ! - person Corwin of Amber; 21.02.2021

arrow_upward
1
arrow_downward

Вы можете условно вернуть предел для маленького x, например:

import matplotlib.pyplot as plt
from numpy import *
epsilon=1e-8

def f(x):
   if x<epsilon
      return 0.5

   return (1-sin(x)/tan(x))/x**2
   #note: same as (1-cos(x))/x**2

x=arange(0,6,0.01)
y=vectorize(f)
plt.plot(x,y(x))
plt.show()

Нанесенная кривая выглядит гладкой

Примечание. Я предпочитаю numpy математике. Vectorize позволяет вызывать функцию с массивом (не очень эффективно, но удобно).

person Stefan    schedule 21.02.2021
comment
Большое спасибо за ваш ответ! - person Corwin of Amber; 21.02.2021

arrow_upward
0
arrow_downward

Проблема заключается в ограниченной точности float и double. Вы должны использовать арифметику с большей точностью, такую ​​как mpfr. Его можно использовать в Python через привязку, например

https://pypi.org/project/gmpy2/

Вот пример, когда я использую его в среде более высокого уровня под названием Sagemath: я использую 100-битную точность:

sage: R = RealField(100)  // 100 bits of precision
sage: def f(x):     #these are all the same function using different identities  
....:    a = (1-(sin(x)/tan(x)))/(x**2)
....:    b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
....:    c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
....:    d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
....:    e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
....:    return a, b, c, d, e
....: 
sage: f(R(1.2e-8))
(0.50000000000000264647624775223,
 0.49999999999999716827551705076,
 0.49999999999999716827551705076,
 0.49999999999999716827551705076,
 0.49999999999999169007478634929)
person hivert    schedule 21.02.2021
comment
Вежливо сказать, в чем проблема, когда минусуют! - person hivert; 21.02.2021
comment
Большое спасибо за ваш ответ! (я не минусовал, а плюсовал) - person Corwin of Amber; 21.02.2021
comment
«Использовать больше точности» — не самый лучший ответ для чисел с плавающей запятой. Надо разбираться в вопросах. Уильям Кахан приписывает следующий пример Жану-Мишелю Мюллеру примерно в 1980 году. Пусть f(y, z) = 108−(815−1500/z)/y. Учитывая x[0] = 4 и x[1] = 4¼, определите x[n+1] = f(x[n], x[n-1]). Вычислите x[80]. Правильный ответ около 5. Все числовые реализации ошибаются; большинство сообщает 100. См. связанный документ, чтобы понять, почему. Хороший ответ на вопрос OP — понять источники ошибок и разработать расчет, который позволяет их избежать или уменьшить. - person Eric Postpischil; 28.02.2021
comment
Чтобы было ясно, смысл «Все числовые реализации ошибаются» заключается в том, что использование большего количества битов не поможет. Использование арифметики с большей точностью не даст правильного ответа. Этот пример, конечно же, был создан для этого, и использование большей точности часто помогает при рутинных вычислениях. Но иногда это просто пинает банку дальше по дороге. - person Eric Postpischil; 28.02.2021
comment
@EricPostpischil: спасибо за указатель. Я попробовал ваш встречный пример с точностью 1000 бит и получил 4,99999999999999999785572107704779534792177... Так что это сложнее. Хотя я согласен с вами в этом принципе, я предположил, что этот вопрос является вопросом реализации, а не вопросом математического/численного анализа, который выходит за рамки SO. - person hivert; 03.03.2021
comment
Математика и численный анализ не выходят за рамки Stack Overflow, по крайней мере как таковые, хотя вопросы почти в любой области могут оказаться слишком общими, чтобы ответить на них. Конечно, если возможны краткие ответы, математический и числовой анализ в том виде, в каком они применимы к программированию. (Есть ли у людей какое-то представление о том, что программирование — это просто кодирование? Информатика связана с компьютерами так же, как астрономия связана с телескопами.) И в вопросе предлагается избежать катастрофической отмены и, в частности, упоминается использование тригонометрического тождества, чтобы избежать этого. - person Eric Postpischil; 04.03.2021