Задачи с плавающей точкой в асимптотических функциях, приближающихся к нулю - Python
Новичок в питоне от MATLAB.
Я использую гиперболическое касательное усечение функции масштаба. Я сталкиваюсь с моей проблемой при применении 0.5 * math.tanh(r/rE-r0) + 0.5 функция на массив значений диапазона r = np.arange(0.1,100.01,0.01), Я получаю несколько 0.0 Значения для функции на стороне приближаются к нулю, что вызывает проблемы с доменом при выполнении логарифма:
P1 = [ (0.5*m.tanh(x / rE + r0 ) + 0.5) for x in r] # truncation function
Я использую этот обходной путь:
P1 = [ -m.log10(x) if x!=0.0 else np.inf for x in P1 ]
этого достаточно для того, что я делаю, но это всего лишь бинтовое решение.
Как просили для математической ясности:
В астрономии шкала величин работает примерно так:
mu = -2.5log(flux) + mzp # apparent magnitude
где mzp - величина, при которой можно увидеть 1 фотон в секунду. Следовательно, большие потоки равняются меньшей (или более отрицательной) видимой величине. Я делаю модели для источников, которые используют многокомпонентные функции. Ex. две sersic функции с различными sersic индексами с P1 внешнее усечение внутреннего компонента и 1-P1 внутреннее усечение внешнего компонента. Таким образом, при добавлении функции усечения к каждому компоненту величина, определяемая радиусом, станет очень большой из-за того, что mu1-2.5 * log (P1) получает как P1 асимптотически приближается к нулю.
TLDR: я хотел бы знать, есть ли способ сохранить числа с плавающей запятой, точность которых недостаточна, чтобы их можно было отличить от нуля (в частности, в результатах функций, которые асимптотически приближаются к нулю). Это важно, потому что при логарифме таких чисел возникает ошибка домена.
Последнее число перед выходом в логарифмическом P1 начинает читать ноль 5.551115123125783e-17, который является обычным результатом арифметической ошибки округления с плавающей запятой, где желаемый результат должен быть равен нулю.
Любой вклад будет принята с благодарностью.
@user: Дэн, не помещая весь мой сценарий:
xc1,yc1 = 103.5150,102.5461;
Ee1 = 23.6781;
re1 = 10.0728*0.187;
n1 = 4.0234;
# radial brightness profile (magnitudes -- really surface brightness but fine in ex.)
mu1 = [ Ee1 + 2.5/m.log(10)*bn(n1)*((x/re1)**(1.0/n1) - 1) for x in r];
# outer truncation
rb1 = 8.0121
drs1 = 11.4792
P1 = [ (0.5*m.tanh( (2.0 - B(rb1,drs1) ) * x / rb1 + B(rb1,drs1) ) + 0.5) for x in r]
P1 = [ -2.5*m.log10(x) if x!=0.0 else np.inf for x in P1 ] # band-aid for problem
mu1t = [x+y for x,y in zip(P1,mu1)] # m1 truncated by P1
где bn(n1)=7,72 и B(rb1,drs1) = 2,65 - 4,98 * (r_b1 / (-drs1));
mu1 - профиль величины усеченного компонента. P1 - функция усечения. Многие из последних записей для P1 равны нулю, что связано с тем, что числа с плавающей запятой неотличимы от нуля из-за точности с плавающей запятой.
Простой способ увидеть проблему:
>>> r = np.arange(0,101,1)
>>> P1 = [0.5*m.tanh(-x)+0.5 for x in r]
>>> P1
[0.5, 0.11920292202211757, 0.01798620996209155, 0.002472623156634768, 0.000335350130466483, 4.539786870244589e-05, 6.144174602207286e-06, 8.315280276560699e-07, 1.1253516207787584e-07, 1.5229979499764568e-08, 2.0611536366565986e-09, 2.789468100949932e-10, 3.775135759553905e-11, 5.109079825871277e-12, 6.914468997365475e-13, 9.35918009759007e-14, 1.2656542480726785e-14, 1.7208456881689926e-15, 2.220446049250313e-16, 5.551115123125783e-17, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Обратите внимание также на числа с плавающей точкой перед нулями.
2 ответа
Напомним, что гиперболический тангенс может быть выражен как (1-e^{-2x})/(1+e^{-2x}). Немного алгебры, мы можем получить, что 0.5*tanh(x)-0.5 (отрицание вашей функции) равно e^{-2x}/(1+e^{-2x}). Логарифм этого будет -2*x-log(1+exp(-2*x)), который будет работать и быть стабильным везде.
То есть я рекомендую вам заменить:
P1 = [ (0.5*m.tanh( (2.0 - B(rb1,drs1) ) * x / rb1 + B(rb1,drs1) ) + 0.5) for x in r]
P1 = [ -2.5*m.log10(x) if x!=0.0 else np.inf for x in P1 ] # band-aid for problem
С помощью этого более простого и стабильного способа сделать это:
r = np.arange(0.1,100.01,0.01)
#r and xvals are numpy arrays, so numpy functions can be applied in one step
xvals=(2.0 - B(rb1,drs1) ) * r / rb1 + B(rb1,drs1)
P1=2*xvals+np.log1p(np.exp(-2*xvals))
Вы можете попробовать две вещи.
(1) метод грубой силы: найдите пакет с плавающей арифметикой переменной точности и используйте его вместо встроенной фиксированной точности. Я играю с вашей проблемой в Maxima [1] и обнаружил, что мне нужно значительно повысить точность с плавающей точкой, чтобы избежать недополнения, но это возможно. Я могу опубликовать код Maxima, если хотите. Я хотел бы представить, что есть некоторая подходящая библиотека с плавающей точкой переменной точности для Python.
(2) приблизительное значение log((1/2)(1 + tanh(-x)) с рядом Тейлора или какое-либо иное приближение, чтобы вообще избежать log (tanh (...)).