Python SciPy special.sici用法及代碼示例

本文簡要介紹 python 語言中 scipy.special.sici 的用法。

用法:

scipy.special.sici(x, out=None) = <ufunc 'sici'>#

正弦和餘弦積分。

正弦積分是

\[\int_0^x \frac{\sin{t}}{t}dt\]

餘弦積分是

\[\gamma + \log(x) + \int_0^x \frac{\cos{t} - 1}{t}dt\]

其中\(\gamma\)是歐拉常數，\(\log\)是對數的主分支[1]。

參數 ：：

x： array_like

計算正弦和餘弦積分的實數或複數點。

out： ndarray 的元組，可選

函數結果的可選輸出數組

返回 ：：

si： 標量或 ndarray

x 處的正弦積分

ci： 標量或 ndarray

x 處的餘弦積分

注意：

對於 x < 0 的實數參數，ci 是餘弦積分的實部。對於這樣的點 ci(x) 和 ci(x + 0j) 相差 1j*pi 倍。

對於實參，該函數是通過調用 Cephes’ 來計算的[2] 思慈常規。對於複雜的參數，該算法基於 Mpmath[3] 斯和詞例程。

參考：

[1] (1,2)

米爾頓·阿布拉莫維茨和艾琳·A·斯特根編輯。包含公式、圖表和數學表格的數學函數手冊。紐約：多佛，1972 年。(參見第 5.2 節。)

[2]

Cephes 數學函數庫，http://www.netlib.org/cephes/

[3]

弗雷 Delhi 克約翰遜等人。 “mpmath：arbitrary-precision 浮點運算的 Python 庫”(0.19 版)http://mpmath.org/

例子：

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.special import sici, exp1

sici 接受實數或複數輸入：

>>> sici(2.5)
(1.7785201734438267, 0.2858711963653835)
>>> sici(2.5 + 3j)
((4.505735874563953+0.06863305018999577j),
(0.0793644206906966-2.935510262937543j))

對於右半平麵中的 z，正弦和餘弦積分與指數積分 E1 相關(在 SciPy 中實現為 scipy.special.exp1 )：

• Si(z) = (E1(i*z) - E1(-i*z))/2i + pi/2

• Ci(z) = -(E1(i*z) + E1(-i*z))/2

請參見 [1](方程 5.2.21 和 5.2.23)。

我們可以驗證這些關係：

>>> z = 2 - 3j
>>> sici(z)
((4.54751388956229-1.3991965806460565j),
(1.408292501520851+2.9836177420296055j))

>>> (exp1(1j*z) - exp1(-1j*z))/2j + np.pi/2  # Same as sine integral
(4.54751388956229-1.3991965806460565j)

>>> -(exp1(1j*z) + exp1(-1j*z))/2            # Same as cosine integral
(1.408292501520851+2.9836177420296055j)

繪製在實軸上評估的函數；水平虛線位於 pi/2 和 -pi/2 處：

>>> x = np.linspace(-16, 16, 150)
>>> si, ci = sici(x)

>>> fig, ax = plt.subplots()
>>> ax.plot(x, si, label='Si(x)')
>>> ax.plot(x, ci, '--', label='Ci(x)')
>>> ax.legend(shadow=True, framealpha=1, loc='upper left')
>>> ax.set_xlabel('x')
>>> ax.set_title('Sine and Cosine Integrals')
>>> ax.axhline(np.pi/2, linestyle=':', alpha=0.5, color='k')
>>> ax.axhline(-np.pi/2, linestyle=':', alpha=0.5, color='k')
>>> ax.grid(True)
>>> plt.show()