Комплексная функция ошибок в Mathematica
Функция комплексной ошибки w(z) определяется как e^(-x^2) erfc(-ix), Проблема с использованием w(z), как определено выше, заключается в том, что erfc имеет тенденцию взрываться при большем x (в дополнение к экспоненте, идущей в 0, так что все остается маленьким), так что Mathematica возвращается к вычислениям произвольной точности, которые делают жизнь ОЧЕНЬ медленной. Эта функция используется для реализации профиля voigt - формы линии, обычно используемой в спектроскопии и других смежных областях. Сейчас я возвращаюсь к расчету формы линии один раз и использую интерполяцию для ускорения, но это не позволяет мне легко изменять параметры формы линии (или подгонять их).
Scipy имеет хорошую и быструю реализацию w(z) как scipy.special.wofzи мне было интересно, есть ли эквивалент в Mathematica.
5 ответов
Сложная функция ошибок может быть записана в терминах полинома Эрмита H_{-1}(x):
In[1]:= FullSimplify[2 HermiteH[-1,I x] == Sqrt[Pi] Exp[-x^2] Erfc[I x]]
Out[1]= True
И оценка не страдает так много недостатков и переполнений
In[68]:= 2 HermiteH[-1, I x] /. x -> 100000.
Out[68]= 6.12323*10^-22 - 0.00001 I
In[69]:= Sqrt[Pi] E^-x^2 Erfc[I x] /. x -> 100000.
During evaluation of In[69]:= General::unfl: Underflow occurred in computation. >>
During evaluation of In[69]:= General::ovfl: Overflow occurred in computation. >>
Out[69]= Indeterminate
Тем не менее, некоторые быстрые тесты показывают, что скорость оценки функции Эрмита будет медленнее, чем скорость произведения экспоненциальной функции и ошибки...
Действительная и мнимая части комплексной функции ошибки на действительной линии могут быть явно и эффективно вычислены в Mathematica с использованием интеграла Доусона:
In[9]:= Sqrt[Pi] Exp[-x^2] Erfc[I x] ==
E^-x^2 Sqrt[\[Pi]] - 2 I DawsonF[x] // FullSimplify
Out[9]= True
Это примерно в 4 раза быстрее, чем при использовании HermiteH[-1,z],
In[10]:= w1[x_] := E^-x^2 Sqrt[\[Pi]] - 2 I DawsonF[x]
w2[x_] := 2 HermiteH[-1, I x]
In[15]:= AbsoluteTiming[w1 /@ Range[-5.0, 5.0, 0.001];]
Out[15]= {2.3272327, Null}
In[16]:= AbsoluteTiming[w2 /@ Range[-5.0, 5.0, 0.001];]
Out[16]= {10.2400239, Null}
Расширение серии на бесконечности показывает, что действительная и мнимая части имеют очень разные масштабы. Я бы посоветовал вычислять их отдельно, а не добавлять. Ниже я использую первые несколько терминов расширения серии, чтобы получить мнимую часть.
In[186]:=
w[x_?NumericQ] := {N[Exp[-SetPrecision[x, 25]^2], 20],
N[(3 /(4 Sqrt[\[Pi]] x^5) + 1/(2 Sqrt[\[Pi]] x^3) + 1/(
Sqrt[\[Pi]] x))]}
In[187]:= w[11]
Out[187]= {2.8207700884601354011*10^-53, 0.05150453151309212}
In[188]:= w[1000]
Out[188]= {3.296831478088558579*10^-434295, 0.0005641898656429712}
Не уверен, как сильно вы хотите эту очень маленькую реальную часть. Если вы можете уронить его, это будет держать цифры в разумных пределах. В некоторых диапазонах (или если требуется точность, превышающая точность станка), вы можете использовать больше терминов из расширения этой воображаемой части.
Даниэль Лихтблау Вольфрам Исследования
Программа на языке C для сложной функции ошибок (она же функция Фаддеева), которую можно запустить из Mathematica, также доступна в RooFit. Прочитайте статью Karbach et al. arXiv: 1407.0748 для получения дополнительной информации.