实现指数函数

2020-06-29 04:05:19

我探索了实现指数函数$f(X)=e^x$的几种复杂的逼近技术,包括泰勒级数逼近、拉格朗日插值、切比雪夫插值、Carathéodory-Fejer逼近和MinimMax逼近。这也是对使用这些方法来近似其他函数的更一般的介绍。在此过程中,我详细介绍了每种方法的相关理论,并应用数值分析方法对各种形式的误差进行了导航。我还以双精度11包括了我自己的每个算法的Python和C++实现:请注意,这里包含的实现将以双精度输出和计算。我们仍然可以在此设置中应用复杂的方法,但作为一般规则,您需要更高的中间精度才能获得与双精度下的真实值难以区分的近似值。这是通常建议使用您的平台提供的优化的exp(X)函数而不是从头开始编写的一个(但不是唯一的)原因。带着沉重的评论。最后,我使用超优化的、平台提供的exp(X)函数作为基准,分析了每个实现的性能和准确性特征。2 2:在大多数情况下,无论您使用哪个数值库,都会在幕后映射到平台提供的exp(X)函数。例如,这就是numpy.exp在幕后做的事情。背景部分首先介绍$e^x$的激励属性和浮点系统的基础知识。每种技术的代码都包含在相应的实现标题下。

为了促进$e$的定义,我们将回顾一些微积分。设$f$是在某个实值$a$处无限可微的函数。这就是说,我们可以重复微分$f$:对于每个$k^{\text{th}}$导数$f^{(K)}$,我们可以再次微分$f^{(K)}$,得到$(k+1)^{\text{th}}$导数$f^{(k+1)}$。任何具有此性质的函数都可以唯一地表示为“以$a$为中心”的泰勒级数展开。以$a$为中心,取值为$x$的$f$的泰勒级数定义为。

$$f(X)=\frac{f(A)}{0!}(x-a)^0+\frac{f‘(A)}{1!}(x-a)^1+\ldots+\frac{f^{(K)}(A)}{k!}(x-a)^k+\ldots$$。

这是幂级数,或一元无限多项式。展开中心确定泰勒级数返回的值的邻域,每个泰勒项的系数是通过重复微分函数$f$并将其求值为$a$来确定的。一个常见的展开中心是$a$=0,在这种情况下,泰勒级数也称为Maclaurin级数,级数围绕原点为中心。这可以被认为是“默认”设置。如果在某些项$k$之后截断泰勒展开式的所有项,就会得到一个次数为$k$的多项式。级数的$k^{\text{th}}$项(如果截断,则为多项式)的系数由下式给出。

作为一个具体的例子,考虑正弦函数的泰勒级数展开。正弦函数不仅是无穷可微的,而且是循环的。

$$\BEGIN{aligned}\sin^{(1)}(X)&;=\cos(X),\sin^{(2)}(X)&;=\cos^{(1)}(X)=-\sin^{(3)}(X)&;=-\sin^{(1)}(X)=-\cos(X),\sin^{(4)}(X)&;=-\cos^{(1)}(X)=\sin(X)\end{aligned}$$。

我们确定泰勒级数的每个$k^{\text{th}}$系数的方法是在$a$处计算$f^{(K)}$,并将其除以$k$的阶乘。如果我们想围绕原点($a$=0)展开正弦函数,就可以得到循环系数。

$$\BEGIN{ALIGNED}\sin(0)&;=0,\sin^{(1)}(0)&;=\cos(0)=1,\sin^{(2)}(0)&;=\cos^{(1)}(0)=-\sin(0)=0,\sin^{(3)}(0)&;=-\sin^{(1)}(0)=-\cos(0)=-1,\sin^{(4)}(0)&;=-\cos^{(1)}(0)=\sin(0)=0\end{aligned}$$

在任意项$k$截断函数$f$的泰勒展开式,即可使用$k$次泰勒多项式有限逼近$f$。当$x$相对接近$a$时,以$a$为中心的$f$的泰勒多项式产生非常精确的$f(X)$的近似值。当绝对值$x$远离$a$时,泰勒多项式的精度迅速降低,这意味着它需要更多的泰勒级数项(即高次多项式)才能精确逼近。考虑下面的图,它显示了在区间$[-20,20]$上$\sin(X)$的值与其以$a$=0为中心的1,3,5,7和9次泰勒多项式的比较。

观察到,当$x$接近$a$=0时,$\sin(X)$的泰勒近似更精确,但很快就偏离真实值$\sin(X)$,远离0。1次泰勒多项式只是$\sin(X)$在原点附近很小的区间内的精确逼近,而9次泰勒多项式在$[-5,5]$内是非常精确的。近似在变得极不准确之前能维持多久取决于泰勒多项式的项数。高次多项式将在更长的时间内保持$\sin(X)$的较好逼近,但任何有限多项式最终都将变得非常不准确。

数学常数$e$(几乎)完全是由它在幂运算下表现出的非常好的性质所驱动的。特别地,$e$的定义是出于寻找一个连续函数的愿望,该函数是它自己的导数,并且将加法恒等式0映射到乘法恒等式1。这是因为用这样的函数来解决困难的积分和微分问题要方便得多。推而广之,应用数学和物理中的很大一部分问题归结为求解微分方程,对微分方程来说,这样的函数是基本的。

碰巧的是,$f(X)=e^x$唯一地满足这一性质。我们可以展示这一点,并在这个过程中直接定义$e$,从任意函数$f$在$a$=0处无限可微的泰勒级数表示开始。设$a_0,a_1,\ldots$是以$a$为中心的$f$的泰勒级数的系数。然后我们有泰勒级数。

从微分的线性推导出一阶导数$f‘$的泰勒级数展开式为。

为了确定一个函数是它自己的导数,我们求解满足$f=f‘$的系数$a_0,a_1,\ldots$:

$$a_0+a_1 x+a_2 x^2+\ldots=a_1+2a_2 x+3a_3 x^2+\ldots$$。

给定$a_0$=1,我们发现作为其自身导数的函数的泰勒级数是

我们用$e^x$表示此函数,其中$e$被定义为此函数在$x$=1时的值。

下图更直观地说明了$e^x$为什么是特殊的,其中具有不同基数的指数函数与它们的导数一起绘制。基数小于$e$的指数函数(如$b$=2)比其导数增长得更快。但是当基数大于$e$时,比如$b$=4,它的增长速度就没有它的导数那么快了。

有一种内在的张力,因为我们希望在不做太多工作的情况下确定$e^x$的精确值。在我们考虑算法的效率之前,我们需要考虑它的准确性。这导致我们定义了各种类型的错误,其中最重要的错误来自我们近似实数的方式。通常不可能计算任意函数$f$的$f(X)$的精确值,因为计算机不能处理任意实数。33:几乎所有的实数都不是可计算的。可计算的实数通常不能精确地表示出所需的精度水平,因为它们要么是无理的(因此具有无限的小数展开式),要么是具有非常长的小数展开式的有理数。我们最多只能把这个值近似到某种可以接受的精度。

IEEE754浮点标准通过将给定社区中所有附近的实值映射到单个舍入值,将实数间隔离散化成可计算的形式。在内部,IEEE 754二进制浮点数$N$使用规范化形式表示。

其中第一位分配给符号(符号位),$p$位$b_1.b_2b_3\ldots b_p$包括尾数或有效数,$E_{k}$是由$k$位组成的整数指数。请注意,由于此形式是标准化的,$b_1=1$,而$b_2,\l点b_p$中的每一个可以等于0或1。IEEE 754单精度二进制浮点数的总大小为$32$位:$8$分配给[-126,127]$中的指数$E\,$23$分配给尾数($p$=24表示归一化位)。因此,您可以用单精度浮点数表示$2^{32}$不同的值,下溢和溢出限制分别为$2^{127}\约3.4倍10^{38}$和$2^{-126}\约1.2倍10^{-38}$。同样,IEEE 754双精度浮点值的总大小为$64$位:$11$位分配给[-1022,1024]$中的指数$E\,$52$位分配给尾数($p=53$)。您可以双精度表示$2^{64}$DISTINCT值,下溢和溢出限制分别为$2^{1024}\约1.8\乘以10^{308}$和$2^{-1022}\约2.2\乘以10^{-308}$。

浮点值分布不均匀。下面的密度图说明了这一点,它绘制了所有16位浮点值与它们的基2指数的关系图。

我们可以看到,值相对密集地聚集在$0$附近,并且越远离原点,值就越稀疏。作为另一个例子,所有$32$位浮点数的一半驻留在实区间$[-1,1]$中。同时,最小和最大的$32$位浮点数分别为$-3.4\x 10^{38}$和$3.4\x 10^{38}$。更一般地,在每个间隔$[2^n,2^{n+1}]$中,可用的浮点值在它们之间以$2^{n-p}$的间距分布,其中$p$是所考虑精度的尾数位数。随着您远离原点的距离越远,您的计算就越有可能在可用浮点值之间撞上实际值,取而代之的是四舍五入为最接近的可用值。

对于任何二进制浮点系统,我们可以使用尾数中可用的位数$p$推导出该系统中可用的最大精度。由于浮点值分布不均匀,因此可用精度会随着远离原点而降低(这会导致近似错误!)。为此,我们定义机器epsilon,用$\epsilon_{M}$表示,它是$1$与大于$1$的最小可表示浮点值之间的差值。在单精度浮点运算中,大于$1$的最小可分值为$2^{-23}$,因此$\epsilon_{M}=2^{-23}$。同样,在双精度下,我们得到$\epsilon_{M}=2^{-52}$。这些系统的最大小数精度可以通过将$\epsilon_{M}$转换为十进制值来获得。综上所述,我们看到$32$位浮点的最大小数精度约为$7$位,而$64$位浮点的最大小数精度约为$16$位。您只能为可用浮点值的子集实现此精度;在距离原点最远的边界,您只能获得$6$十进制位的精度

根本错误。当模型本质上以一种不平凡的方式偏离现实时,就会出现这种情况。例如,流行病学中的Kermack-McKendrick SIR模型没有默认考虑再感染的可能性,因此如果恢复的人群也容易感染,那么它将显示出重大的根本性错误。就我们的目的而言,我们不需要担心根本错误。

浮点错误。这种类型的错误可能以多种方式发生。在最直接的意义上,您无法获得比浮点系统的精度所提供的精度更高的精度。任何小数位数超过16位的实数都不能完全精确地表示为双精度。您可能会发现通过相等比较之类的操作引入了更细微的错误,因为每个浮点数对应着无限多个四舍五入到它的实值。如果在计算过程中,中间结果超出下溢/上溢限制,即使完成计算会将其带回可用值范围,也会发生浮点错误。

离散化错误。这也称为截断误差。每当你取一个连续函数,并以有限多个值进行近似时,它就会发生。离散化误差的一个示例是多项式插值,其中使用$n+1$点用次数为$n$的多项式$P_{n}$逼近连续函数$f(X)$。当使用函数的泰勒级数近似函数时,会出现一个相关的例子。泰勒级数不能在有限时间内求值,因为它有无限多项,所以你必须在某个项处截断级数$n$,这将导致一个次数为$n$的多项式$T_n$。截断后切下的泰勒级数的无穷多剩余项构成了计算中的离散化误差。

收敛错误。这也被称为逐渐失去意义。当您的计算执行太多工作并返回超出可用精度的结果时,就会发生这种情况。这与浮点错误不同,因为它与迭代重复有关。初始结果可能是可表示的(但四舍五入),但重复迭代会合成初始误差,直到它超过精度阈值,此时它不再增加精度,而是降低精度。

$$1+\frac{1}{2}\epsilon_{M}-1=0\neq 1.11022\x 10^{-16}=1-1+\frac{1}{2}\epsilon_{M}$$。

$$0.1+(0.2+0.3)=0.600\ldots 0 0 1\NEQ 0.6=(0.1+0.2)+0.3$$。

x=10/9对于范围(20)内的k:x-=1 x*=10打印(X)。

这是一个失去意义的教科书例子。每个结果总共有16位小数位,但结果的精度会随着时间的推移而降低。请注意,使用起始值$x=\frac{10}{9}$,在一次迭代之后,我们将拥有

$$f(X)=(x-1)\CDOT 10=\Left(\FRAC{10}{9}-1\f25 Right)\CDOT 10=\FRAC{1}{9}\CDOT 10=\FRAC{10}{9}$$。

因此,$x$的值应该在每次迭代后保持不变。但是取而代之的是,随着每一次迭代,我们的值都会迅速偏离真实值。

最后,所有这些类型的错误都出现在同一近似值中的情况并不少见。假设你正在用它的泰勒级数逼近一个连续函数$f(X)$。通过为级数选择一个截断点,您会立即遇到离散化误差,从而为您的近似的总误差设定一个下限。然后,根据泰勒级数近似的实现,您对单个泰勒项的中间计算可能会随着时间的推移而受到复合误差的影响。您的结果可能不能精确地表示为浮点数,这将通过舍入误差进一步降低计算的准确性。

我们有两种测量误差的方法。给定函数$f$及其近似值$F$,在某个$x$处的绝对误差由下式给出。

相对误差取决于并归一化绝对误差。在评估非常广泛的可能值范围内的近似精度时,它非常有用。例如,假设我们有一个函数$f$的近似值$F$,该函数返回值介于0和4,294,967,296之间。对于该范围内的小值,绝对误差$\epsilon$=10基本上毫无价值,但对于接近上限的较大值,则相当不错。为了更好地了解近似函数在这样的范围内有多精确,我们可以考虑相对误差。这伴随着一个警告,即当$x$处的$f$的真值正好为0时,相对误差在$x$处是未定义的,因为否则我们将除以$0$。幸运的是,我们在这里不会遇到这个问题,因为没有满足$e^x$=0的$x$的值。

绝对误差和相对误差给我们提供了一种严格量化近似精度的方法。如果我们可以证明其中一个误差度量的上界,我们就可以断言近似的最坏情况下的精度。此外,这些误差度量适用于常见的汇总统计数据,因此我们可以使用它们来确定平均精度,或按分位数确定精度。相对误差也为我们提供了一种用精度位数来评估我们计算的准确性的方法。对于任何相对误差$\eta$,相应的精度是满足以下条件的最大整数$p$。

其中$\beta$是正在考虑的基数。十进制精度指定$\beta$=10,二进制精度指定$\beta$=2。

为了避免处理下溢和溢出限制,我们将在域$[-709,709]$中为$x$实现$e^x$。这些值接近可以双精度计算的极限(尽管从技术上讲,如果我们不关心是否具有对称域,我们可以将$x$计算到-740左右)。在开始实现之前,我们需要建立一个良好的基准函数来与之进行比较。我们将在Python中使用NumPy提供的exp(X)函数,在C++中使用标准库。

接下来,我们将编写一个REL_ERROR报告函数,它将返回我们的实现相对于基准的相对错误。由于我们正在评估值的区间,因此此函数将使用基准给定的“真”函数$f$、我们的近似值$F(X)$、要运算的区间$[a,b]$的下确界$a$和上确界$b$,以及要在区间内计算的值$x$的个数$n$作为其输入。它将返回我们评估的每个$x$值的相对误差向量。此函数的输出还将用于使用matplotlib以图形方式评估相对误差。

我们还需要一些方便的函数,从linspace开始。linspace函数将接受区间的界限$a、b$和正整数$n$作为输入,然后返回由区间中的$n$线性间隔浮点值组成的向量。NumPy提供了这个开箱即用的函数,但是我们将为C++测试编写我们自己的函数。我们还将编写几个C++函数来计算NumPy提供的汇总统计信息。

Python导入numpy为npdef rel_error(f,F,a,b,n,lot=false,save=false,file=NONE):";";";";评估两个函数相对误差的基准函数。首先通过理想函数,然后通过测试中的函数。:param f:要比较的基准函数。:param F:我们要测试的函数的实现。:param a:要考虑的区间的左边界。:param b:要考虑的区间的右边界。:param n:要查看的间隔中的值数。:param Plot:布尔值,是否绘制matplotlib中的错误。:param save:boolean,是将绘图保存到磁盘还是在标准输出中显示。:param file:要保存绘图的文件位置的字符串名称。如果save为true,file为:param float64:表示返回精度的布尔标志。用法:RETURNS:包含以下内容的列表:区间[a,b]内k个线性间隔的x值中每个值的相对误差数组|f(X)-F(X)|/f(X);最大相对误差;平均相对误差;最大迭代次数;平均迭代次数。";";";x=np.linspace(a,b,n)Trials=np.zeros(x.form)in Enumerate(X):Trials[i]=F(。

..