由于毕设的数学推导中涉及了 e r f \mathrm{erf} erf 函数，关于其他函数的渐近计算推导见链接类指数级数（指数积分函数的变体）数值计算算法的C++实现。
反正闲得无聊，虽然知道这种函数肯定有现成的轮子了，然而我是情报弱者。再加上最后我的算法是要在 C++ 平台上进行实现的，不如自己造一手轮子。
注意1：因为我的场景只涉及 x ⩾ 0 x\geqslant0 x⩾0 的情形，所以只针对这种情况进行了考虑。事实上，根据对称性 x < 0 x<0 x<0 ，直接用 e r f ( x ) = 1 − e r f ( − x ) \mathrm{erf}\left(x\right)=1-\mathrm{erf}\left(-x\right) erf(x)=1−erf(−x) 即可
注意2：这里我采用的是局部展开，所以最坏复杂度可能会很高。正经的做法使用Remez算法进行展开，详见链接https://baike.baidu.com/item/%E9%9B%B7%E7%B1%B3%E5%85%B9%E7%AE%97%E6%B3%95/23665981?fr=aladdin，也可参考FDLIBM的C源码实现，见http://gnuwin32.sourceforge.net/packages/fdlibm.htm
C++代码注：为了与已有函数区分，函数名前加了Hsk，Hsk是我的ID的缩写(羽咲->Hasaki->Hsk)
话不多说，先给代码再解释。这里的计算精度为
ε
<
1
0
−
10
\varepsilon<10^{-10}
ε<10−10 ，如果想要修改计算精度，应当调节 HSK_EPSILON 和 HSK_INF_SERIES 这两个宏。
#includee r f ( x ) = 2 π ∫ 0 x exp ( − t 2 ) d t (4.1) \mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi}}\int_0^x{\exp\left(-t^2\right)\mathrm{d}t}\tag{4.1} erf(x)=π 2∫0xexp(−t2)dt(4.1)
原点展开令
a
n
(
x
)
=
(
−
1
)
n
x
2
n
+
1
n
!
(
2
n
+
1
)
(4.2)
a_n\left( x \right) =\frac{\left( -1 \right) ^nx^{2n+1}}{n!\left(2n+1\right)}\tag{4.2}
an(x)=n!(2n+1)(−1)nx2n+1(4.2)
有递推式
a
n
(
x
)
=
(
−
1
)
n
x
2
n
+
1
n
!
(
2
n
+
1
)
=
−
(
2
n
−
1
)
x
2
(
2
n
+
1
)
n
(
−
1
)
n
−
1
x
2
n
−
1
(
n
−
1
)
!
(
2
n
−
1
)
=
(
1
−
2
2
n
+
1
)
−
x
2
n
a
n
(
x
)
(4.3)
\begin{aligned} a_n\left( x \right) &=\frac{\left( -1 \right) ^nx^{2n+1}}{n!\left( 2n+1 \right)}\ &=-\frac{\left( 2n-1 \right) x^2}{\left( 2n+1 \right) n}\frac{\left( -1 \right) ^{n-1}x^{2n-1}}{\left( n-1 \right) !\left( 2n-1 \right)}\ &=\left( 1-\frac{2}{2n+1} \right) \frac{-x^2}{n}a_n\left( x \right)\ \end{aligned}\tag{4.3}
an(x)=n!(2n+1)(−1)nx2n+1=−(2n+1)n(2n−1)x2(n−1)!(2n−1)(−1)n−1x2n−1=(1−2n+12)n−x2an(x)(4.3)
原点展开式为
e
r
f
(
x
)
=
2
π
∫
0
x
exp
(
−
t
2
)
d
t
=
2
π
∑
n
=
0
∞
∫
0
x
(
−
1
)
n
t
2
n
n
!
d
t
=
2
π
∑
n
=
0
∞
(
−
1
)
n
x
2
n
+
1
n
!
(
2
n
+
1
)
=
2
π
∑
n
=
0
∞
a
n
(
x
)
(4.4)
\begin{aligned} \mathrm{erf}\left( x \right) &=\frac{2}{\sqrt{\pi}}\int_0^x{\exp \left( -t^2 \right) \mathrm{d}t}\ &=\frac{2}{\sqrt{\pi}}\sum_{n=0}^{\infty}{\int_0^x{\frac{\left( -1 \right) ^nt^{2n}}{n!}\mathrm{d}t}}\ &=\frac{2}{\sqrt{\pi}}\sum_{n=0}^{\infty}{\frac{\left( -1 \right) ^nx^{2n+1}}{n!\left(2n+1\right)}}\ &=\frac{2}{\sqrt{\pi}}\sum_{n=0}^{\infty}{a_n\left( x \right)}\ \end{aligned} \tag{4.4}
erf(x)=π
2∫0xexp(−t2)dt=π
2n=0∑∞∫0xn!(−1)nt2ndt=π
2n=0∑∞n!(2n+1)(−1)nx2n+1=π
2n=0∑∞an(x)(4.4)
余项估计（假设
N
+
2
>
2
x
2
N+2> 2x^2
N+2>2x2 ）(详见类指数级数（指数积分函数的变体）数值计算算法的C++实现的第三章）
∣
a
N
+
1
+
p
(
x
)
∣
=
∣
x
∣
2
(
N
+
1
+
p
)
+
1
(
N
+
1
+
p
)
!
=
∣
x
∣
2
N
+
1
(
N
+
1
)
!
(
x
2
)
p
(
N
+
2
)
⋯
(
N
+
1
+
p
)
⩽
∣
x
∣
2
N
+
1
(
N
+
1
)
!
(
x
2
N
+
2
)
p
<
∣
a
N
+
1
(
x
)
∣
×
(
1
2
)
p
(4.5)
\begin{aligned} \left| a_{N+1+p}\left( x \right) \right| &=\frac{\left| x\right| ^{2\left(N+1+p\right)+1}}{\left( N+1+p \right) !}\ &=\frac{\left| x \right|^{2N+1}}{\left( N+1 \right) !}\frac{\left( x^2 \right)^p}{\left( N+2 \right) \cdots \left( N+1+p \right)}\ &\leqslant \frac{\left| x \right|^{2N+1}}{\left( N+1 \right) !}\left( \frac{x^2}{N+2} \right) ^p\ &<\left| a_{N+1}\left( x \right) \right| \times \left( \frac{1}{2} \right) ^p\ \end{aligned} \tag{4.5}
∣aN+1+p(x)∣=(N+1+p)!∣x∣2(N+1+p)+1=(N+1)!∣x∣2N+1(N+2)⋯(N+1+p)(x2)p⩽(N+1)!∣x∣2N+1(N+2x2)p<∣aN+1(x)∣×(21)p(4.5)
2 ∣ a N + 1 ( x ) ∣ ⩾ 2 ∣ a N + 1 ( N + 2 2 ) ∣ = 2 ( N + 2 2 ) 2 N + 3 ( N + 1 ) ! > 2 (4.6) \begin{aligned} 2\left| a_{N+1}\left( x \right) \right|&\geqslant 2\left| a_{N+1}\left( \sqrt{\frac{N+2}{2}} \right) \right|\\ &=\frac{2\left( \sqrt{\frac{N+2}{2}} \right) ^{2N+3}}{\left( N+1 \right) !}\\ &>2\\ \end{aligned} \tag{4.6} 2∣aN+1(x)∣⩾2∣∣∣∣∣aN+1(2N+2 )∣∣∣∣∣=(N+1)!2(2N+2 )2N+3>2(4.6)
无穷远点展开渐近展开
e
r
f
(
x
)
=
2
π
∫
0
x
exp
(
−
t
2
)
d
t
=
2
π
(
π
2
−
exp
(
−
x
2
)
∑
n
=
0
∞
(
−
1
)
n
(
2
n
−
1
)
!
!
2
n
+
1
1
x
2
n
+
1
)
=
1
−
2
π
exp
(
−
x
2
)
∑
n
=
0
∞
(
−
1
)
n
(
2
n
−
1
)
!
!
2
n
+
1
1
x
2
n
+
1
(4.8)
\begin{aligned} \mathrm{erf}\left( x \right) &=\frac{2}{\sqrt{\pi}}\int_0^x{\exp \left( -t^2 \right) \mathrm{d}t}\ &=\frac{2}{\sqrt{\pi}}\left( \frac{\sqrt{\pi}}{2}-\exp \left( -x^2 \right) \sum_{n=0}^{\infty}{\left( -1 \right) ^n\frac{\left( 2n-1 \right) !!}{2^{n+1}}\frac{1}{x^{2n+1}}} \right)\ &=1-\frac{2}{\sqrt{\pi}}\exp \left( -x^2 \right) \sum_{n=0}^{\infty}{\left( -1 \right) ^n\frac{\left( 2n-1 \right) !!}{2^{n+1}}\frac{1}{x^{2n+1}}}\ \end{aligned}\tag{4.8}
erf(x)=π
2∫0xexp(−t2)dt=π
2(2π
−exp(−x2)n=0∑∞(−1)n2n+1(2n−1)!!x2n+11)=1−π
2exp(−x2)n=0∑∞(−1)n2n+1(2n−1)!!x2n+11(4.8)
注意到无穷展开式并非展开阶数越高越好，当阶数较大时，系数会以比指数更快的速度发散。经过 mathematica 测试，发现展开到
n
=
4
n=4
n=4 （即
x
−
9
x^{-9}
x−9 ）项时效果最好。因此实际渐近展开式为
e
r
f
(
x
)
≈
1
−
2
π
exp
(
−
x
2
)
[
1
2
x
−
1
4
x
3
+
3
8
x
5
−
15
16
x
7
+
105
32
x
9
]
(4.9)
\mathrm{erf}\left( x \right) \approx 1-\frac{2}{\sqrt{\pi}}\exp \left( -x^2 \right) \left[ \frac{1}{2x}-\frac{1}{4x^3}+\frac{3}{8x^5}-\frac{15}{16x^7}+\frac{105}{32x^9} \right] \tag{4.9}
erf(x)≈1−π
2exp(−x2)[2x1−4x31+8x53−16x715+32x9105](4.9)
其误差近似小于
1
0
−
4
(
x
−
1
)
10^{-4\left(x-1\right)}
10−4(x−1) ，我们可以通过 mathematica 画图证明这一近似误差估计式，其拟合效果非常好
Plot[{Log10[
Abs[Erf[x] - (1 -
2/Sqrt[Pi]
Exp[-x^2] (1/(2 x) - 1/(4 x^3) + 3/(8 x^5) - 15/(16 x^7) +
105/(32 x^9)))]], -4 x + 4}, {x, 1, 5}]
也就是说，当 x ⩾ log 10 ( 1 / ε ) 4 + 1 x\geqslant \frac{\log_{10}\left(1/\varepsilon\right)}{4}+1 x⩾4log10(1/ε)+1 时使用无穷渐近展开，反之使用原点渐近展开。
结语本节采用泰勒展开的方式推导了 erf 函数的展开式，这一展开式在
x
x
x 较小和较大时计算速度都较快，但
2
<
x
<
3.5
250 次。更好的计算方法应当为全局拟合，采取预计算常数硬编码的方式，但这不在本文讨论范围内，感兴趣的读者可以参考这篇文章https://www.zhihu.com/answer/1686069171。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）