类指数级数(指数积分函数的变体)数值计算算法的C++实现

类指数级数(指数积分函数的变体)数值计算算法的C++实现,第1张

文章目录
    • 前言
    • HskEta函数
    • HskKsi函数
    • 广义HskEta函数
      • 定义式
      • 主项分析
      • 余项分析

前言

由于毕设的数学推导中涉及了 ∑ n = 1 ∞ x n n ! × n \sum_{n=1}^{\infty}{\frac{x^n}{n!\times n}} n=1n!×nxn ∑ n = 1 ∞ x n n ! × n 2 \sum_{n=1}^{\infty}{\frac{x^n}{n!\times n^2}} n=1n!×n2xn 这两个神奇的函数,其中 ∑ n = 1 ∞ x n n ! × n \sum_{n=1}^{\infty}{\frac{x^n}{n!\times n}} n=1n!×nxn 倒是和指数积分函数 E i ( x ) \mathrm{Ei}\left(x\right) Ei(x) 有关,后面的以及推广形式就没有什么特殊的地方了。于是我给这两个函数命名为了 η ( x ) = ∑ n = 1 ∞ x n n ! × n \eta \left( x \right) =\sum_{n=1}^{\infty}{\frac{x^n}{n!\times n}} η(x)=n=1n!×nxn ξ ( x ) = ∑ n = 1 ∞ x n n ! × n 2 \xi \left( x \right) =\sum_{n=1}^{\infty}{\frac{x^n}{n!\times n^2}} ξ(x)=n=1n!×n2xn

反正闲得无聊,最后我的算法是要在 C++ 平台上进行实现的,不如自己造一手轮子。所以顺便研究了一般形式 ∑ n = 1 ∞ x n n ! × n k \sum_{n=1}^{\infty}{\frac{x^n}{n!\times n^k}} n=1n!×nkxn 的递推计算表达式,并用数学证明给出了余项估计。

注意1:因为我的场景只涉及 x ⩾ 0 x\geqslant0 x0 的情形,所以没有取绝对值,如果要应用到一般场景,需要在判断条件处加绝对值

注意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

注:为了与已有函数区分,函数名前均加了Hsk,Hsk是我的ID的缩写(羽咲->Hasaki->Hsk)

HskEta函数


a n ( x ) = x n n ! × n (1.1) a_n\left( x \right) =\frac{x^n}{n!\times n}\tag{1.1} an(x)=n!×nxn(1.1)

η ( x ) = ∑ n = 1 ∞ x n n ! × n = ∑ n = 1 N a n ( x ) + ∑ n = N + 1 ∞ a n ( x ) = η N ( x ) + R N + 1 ( x ) (1.2) \begin{aligned} \eta \left( x \right) &=\sum_{n=1}^{\infty}{\frac{x^n}{n!\times n}}\ &=\sum_{n=1}^N{a_n\left( x \right)}+\sum_{n=N+1}^{\infty}{a_n\left( x \right)}\ &=\eta _N\left( x \right) +R_{N+1}\left( x \right)\ \end{aligned}\tag{1.2} η(x)=n=1n!×nxn=n=1Nan(x)+n=N+1an(x)=ηN(x)+RN+1(x)(1.2)
其中前面的 η N ( x ) \eta _N\left( x \right) ηN(x) 称为主项, R N + 1 ( x ) R_{N+1}\left( x \right) RN+1(x) 称为余项,我们将在第3章证明在给定误差 ε \varepsilon ε 较小时有 ∣ R N ( x ) ∣ < 2 ∣ a N + 1 ( x ) ∣ \left| R_N\left( x \right) \right|<2\left| a_{N+1}\left( x \right) \right| RN(x)<2aN+1(x) 成立,并证明了递推表达式 a n − 1 ( x ) × [ ( 1 − 1 n ) 1 n ] x a_{n-1}\left( x \right)\times \left[ \left( 1-\frac{1}{n} \right) \frac{1}{n} \right] x an1(x)×[(1n1)n1]x

因此写出代码

#define HSK_EPSILON 0.000001

double HskEta (double x) {
    double ans = 0;
    double an = x;      // a1(x)
    
    double n;
    for (n = 2; 2 * an > HSK_EPSILON; ++ n) {
        ans += an;
        double coff = 1 / n;
        an *= coff * (1 - coff) * x;  // a{n}(x) <- a{n-1}(x)
    }
    return ans;
}
HskKsi函数


a n ( x ) = x n n ! × n 2 (2.1) a_n\left( x \right) =\frac{x^n}{n!\times n^2}\tag{2.1} an(x)=n!×n2xn(2.1)

ξ ( x ) = ∑ n = 1 ∞ x n n ! × n = ∑ n = 1 N a n ( x ) + ∑ n = N + 1 ∞ a n ( x ) = ξ N ( x ) + R N + 1 ( x ) (2.2) \begin{aligned} \xi \left( x \right) &=\sum_{n=1}^{\infty}{\frac{x^n}{n!\times n}}\ &=\sum_{n=1}^N{a_n\left( x \right)}+\sum_{n=N+1}^{\infty}{a_n\left( x \right)}\ &=\xi _N\left( x \right) +R_{N+1}\left( x \right)\ \end{aligned}\tag{2.2} ξ(x)=n=1n!×nxn=n=1Nan(x)+n=N+1an(x)=ξN(x)+RN+1(x)(2.2)
其中前面的 ξ N ( x ) \xi _N\left( x \right) ξN(x) 称为主项, R N + 1 ( x ) R_{N+1}\left( x \right) RN+1(x) 称为余项,我们将在第3章证明在给定误差 ε \varepsilon ε 较小时有 ∣ R N ( x ) ∣ < 2 ∣ a N + 1 ( x ) ∣ \left| R_N\left( x \right) \right|<2\left| a_{N+1}\left( x \right) \right| RN(x)<2aN+1(x) 成立,并证明了递推表达式 a n − 1 ( x ) × [ ( 1 − 1 n ) 2 1 n ] x a_{n-1}\left( x \right)\times \left[ \left( 1-\frac{1}{n} \right) ^2\frac{1}{n} \right] x an1(x)×[(1n1)2n1]x

因此写出代码

#define HSK_EPSILON 0.000001

double HskKsi (double x) {
    double ans = 0;
    double an = x;  // a1(x)
    
    double n;
    for (n = 2; 2 * an > HSK_EPSILON; ++ n) {
        ans += an;
        double coff = 1 / n;
        an *= coff * (1 - coff) * (1 - coff) * x;  // a{n}(x) <- a{n-1}(x)
    }
    return ans;
}
广义HskEta函数

这部分给出上面内容的数学推导与严格证明。

定义式


a n ( x ) = x n n ! × n k (3.1) a_n\left( x \right) =\frac{x^n}{n!\times n^k}\tag{3.1} an(x)=n!×nkxn(3.1)

η ( x ) = ∑ n = 1 ∞ x n n ! × n k = ∑ n = 1 N a n ( x ) + ∑ n = N + 1 ∞ a n ( x ) = η N ( x ) + R N + 1 ( x ) (3.2) \begin{aligned} \eta \left( x \right) &=\sum_{n=1}^{\infty}{\frac{x^n}{n!\times n^k}}\ &=\sum_{n=1}^N{a_n\left( x \right)}+\sum_{n=N+1}^{\infty}{a_n\left( x \right)}\ &=\eta _N\left( x \right) +R_{N+1}\left( x \right)\ \end{aligned}\tag{3.2} η(x)=n=1n!×nkxn=n=1Nan(x)+n=N+1an(x)=ηN(x)+RN+1(x)(3.2)

主项分析

主项迭代式
a n ( x ) = x n − 1 ( n − 1 ) ! × ( n − 1 ) k × ( n − 1 ) k n k + 1 x = x n − 1 ( n − 1 ) ! × ( n − 1 ) k × [ ( 1 − 1 n ) k 1 n ] x = a n − 1 ( x ) × [ ( 1 − 1 n ) k 1 n ] x (3.3) \begin{aligned} a_n\left( x \right) &=\frac{x^{n-1}}{\left( n-1 \right) !\times \left( n-1 \right) ^k}\times \frac{\left( n-1 \right) ^k}{n^{k+1}}x\ &=\frac{x^{n-1}}{\left( n-1 \right) !\times \left( n-1 \right) ^k}\times \left[ \left( 1-\frac{1}{n} \right) ^k\frac{1}{n} \right] x\ &=a_{n-1}\left( x \right)\times \left[ \left( 1-\frac{1}{n} \right) ^k\frac{1}{n} \right] x\ \end{aligned}\tag{3.3} an(x)=(n1)!×(n1)kxn1×nk+1(n1)kx=(n1)!×(n1)kxn1×[(1n1)kn1]x=an1(x)×[(1n1)kn1]x(3.3)
其中主项
η N ( x ) = ∑ n = 1 N a n ( x ) (3.4) \eta _N\left( x \right) =\sum_{n=1}^N{a_n\left( x \right)}\tag{3.4} ηN(x)=n=1Nan(x)(3.4)

余项分析

我们要求余项 $\left| R_{N+1}\left( x \right) \right|<\varepsilon $ ,其中 ε \varepsilon ε 代表误差是一个很小的量,不妨设 ε < 0.1 \varepsilon<0.1 ε<0.1

我们将证明:若 2 ∣ a N + 1 ( x ) ∣ < 0.1 2\left| a_{N+1}\left( x \right) \right|<0.1 2aN+1(x)<0.1 k = 1 , 2 k=1,2 k=1,2 ∣ R N ( x ) ∣ < 2 ∣ a N + 1 ( x ) ∣ \left| R_N\left( x \right) \right|<2\left| a_{N+1}\left( x \right) \right| RN(x)<2aN+1(x) 成立。

证明:

①若 N + 2 > 2 ∣ x ∣ N+2> 2\left| x \right| N+2>2x

∣ x ∣ N + 2 < 1 2 \frac{\left| x \right|}{N+2}< \frac{1}{2} N+2x<21 ,则有下放缩式成立
∣ a N + 1 + p ( x ) ∣ = ∣ x ∣ N + 1 + p ( N + 1 + p ) ! ( N + 1 + p ) k = ∣ x ∣ N + 1 ( N + 1 ) ! ( N + 1 + p ) k ∣ x ∣ p ( N + 2 ) ⋯ ( N + 1 + p ) ⩽ ∣ x ∣ N + 1 ( N + 1 ) ! ( N + 1 ) k ( ∣ x ∣ N + 2 ) p < ∣ a N + 1 ( x ) ∣ × ( 1 2 ) p (3.5) \begin{aligned} \left| a_{N+1+p}\left( x \right) \right| &=\frac{\left| x \right|^{N+1+p}}{\left( N+1+p \right) !\left( N+1+p \right) ^k}\ &=\frac{\left| x \right|^{N+1}}{\left( N+1 \right) !\left( N+1+p \right) ^k}\frac{\left| x \right|^p}{\left( N+2 \right) \cdots \left( N+1+p \right)}\ &\leqslant \frac{\left| x \right|^{N+1}}{\left( N+1 \right) !\left( N+1 \right) ^k}\left( \frac{\left| x \right|}{N+2} \right) ^p\ &<\left| a_{N+1}\left( x \right) \right| \times \left( \frac{1}{2} \right) ^p\ \end{aligned} \tag{3.5} aN+1+p(x)=(N+1+p)!(N+1+p)kxN+1+p=(N+1)!(N+1+p)kxN+1(N+2)(N+1+p)xp(N+1)!(N+1)kxN+1(N+2x)p<aN+1(x)×(21)p(3.5)
故可对余项进行放缩
∣ R N + 1 ( x ) ∣ = ∑ n = N + 1 ∞ a n ( x ) = ∑ p = 0 ∞ a N + 1 + p ( x ) < ∑ p = 0 ∞ ∣ a N + 1 ( x ) ∣ × ( 1 2 ) p = 2 ∣ a N + 1 ( x ) ∣ (3.6) \begin{aligned} \left| R_{N+1}\left( x \right) \right| &=\sum_{n=N+1}^{\infty}{a_n\left( x \right)}\ &=\sum_{p=0}^{\infty}{a_{N+1+p}\left( x \right)}\ &< \sum_{p=0}^{\infty}{\left| a_{N+1}\left( x \right) \right| \times \left( \frac{1}{2} \right) ^p}\ &=2\left| a_{N+1}\left( x \right) \right|\ \end{aligned} \tag{3.6} RN+1(x)=n=N+1an(x)=p=0aN+1+p(x)<p=0aN+1(x)×(21)p=2aN+1(x)(3.6)
②若 N + 2 ⩽ 2 ∣ x ∣ N+2\leqslant 2\left| x \right| N+22x

此时(对于 k = 1 , 2 k=1,2 k=1,2
2 ∣ a N + 1 ( x ) ∣ ⩾ 2 ∣ a N + 1 ( N + 2 2 ) ∣ = 2 ( N + 2 2 ) N + 1 ( N + 1 ) ! ( N + 1 ) k > 0.1 (3.7) \begin{aligned} 2\left| a_{N+1}\left( x \right) \right|&\geqslant 2\left| a_{N+1}\left( \frac{N+2}{2} \right) \right|\\ &=\frac{2\left( \frac{N+2}{2} \right) ^{N+1}}{\left( N+1 \right) !\left( N+1 \right) ^k}\\ &>0.1\\ \end{aligned} \tag{3.7} 2aN+1(x)2aN+1(2N+2)=(N+1)!(N+1)k2(2N+2)N+1>0.1(3.7)
关于这里第3行的不等式的证明会十分繁琐,我们可以用 mathematica 画图来展示这一点的正确性,并且展示其他 k k k 值对 ε \varepsilon ε 的要求。

Table[NMinimize[{(2 ((x + 2)/2)^(x + 1))/((x + 1)! (x + 1)^k), x > 0, 
   x < 50}, x], {k, 0, 10}]
Plot[(2 ((x + 2)/2)^(x + 1))/((x + 1)! (x + 1)^1), {x, 0, 10}]

这与前提 2 ∣ a N + 1 ( x ) ∣ < 0.1 2\left| a_{N+1}\left( x \right) \right|<0.1 2aN+1(x)<0.1 相矛盾,舍去。

证毕。

这意味着在对于 ε < 0.1 \varepsilon<0.1 ε<0.1 的情形总有 ∣ R N ( x ) ∣ < 2 ∣ a N + 1 ( x ) ∣ \left| R_N\left( x \right) \right|<2\left| a_{N+1}\left( x \right) \right| RN(x)<2aN+1(x) 成立,所以我们可以用 2 ∣ a N + 1 ( x ) ∣ 2\left| a_{N+1}\left( x \right) \right| 2aN+1(x) 作为余项的上界。事实上,可以用斯塔林公式证明最紧的渐近上界为 e e − 1 ∣ a N + 1 ( x ) ∣ \frac{\mathrm{e}}{\mathrm{e}-1}\left| a_{N+1}\left( x \right) \right| e1eaN+1(x) ,与这里的结论相差无几。

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/langs/799453.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2022-05-06
下一篇 2022-05-06

发表评论

登录后才能评论

评论列表(0条)

保存