高斯误差函数erf的数值计算方法(C++实现)
HskErf函数
前言
由于毕设的数学推导中涉及了 erf\mathrm{erf}erf 函数,关于其他函数的渐近计算推导见链接类指数级数(指数积分函数的变体)数值计算算法的C++实现。
反正闲得无聊,虽然知道这种函数肯定有现成的轮子了,然而我是情报弱者。再加上最后我的算法是要在 C++
平台上进行实现的,不如自己造一手轮子。
注意1:因为我的场景只涉及 x⩾0x\geqslant0x⩾0 的情形,所以只针对这种情况进行了考虑。事实上,根据对称性 x<0x<0x<0 ,直接用 erf(x)=1−erf(−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
注:为了与已有函数区分,函数名前加了Hsk,Hsk是我的ID的缩写(羽咲->Hasaki->Hsk)
C++代码
话不多说,先给代码再解释。这里的计算精度为 ε<10−10\varepsilon<10^{-10}ε<10−10 ,如果想要修改计算精度,应当调节 HSK_EPSILON
和 HSK_INF_SERIES
这两个宏。
#include<cmath>#define HSK_EPSILON 1e-10
#define HSK_INF_SERIES 3.5 // HSK_INF_SERIES = log10(1/HSK_EPSILON) / 4 + 1
#define HSK_PI 3.14159265358979323846double HskErf (double x) {double ans = 0;static double k = 2 / sqrt(HSK_PI);// 使用无穷展开式if (x >= HSK_INF_SERIES) {double x2 = - 1 / x / x;double an = k * exp(-x * x) / x / 2; // a1double n;for (n = 1; n < 5; ++ n) {ans += an;an *= (n - 0.5) * x2; // a{n}(x) <- a{n-1}(x)}return 1 - ans;}// 使用原点展开式else {double x2 = - x * x;double an = k * x;double n;for (n = 1; 2 * abs(an) > HSK_EPSILON; ++ n) {ans += an;an *= (1 - 1 / (n + 0.5)) * x2 / n; // a{n+1}(x) <- a{n}(x)}return ans;}
}
定义式
erf(x)=2π∫0xexp(−t2)dt(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)
原点展开
令
an(x)=(−1)nx2n+1n!(2n+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)
有递推式
an(x)=(−1)nx2n+1n!(2n+1)=−(2n−1)x2(2n+1)n(−1)n−1x2n−1(n−1)!(2n−1)=(1−22n+1)−x2nan(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)
原点展开式为
erf(x)=2π∫0xexp(−t2)dt=2π∑n=0∞∫0x(−1)nt2nn!dt=2π∑n=0∞(−1)nx2n+1n!(2n+1)=2π∑n=0∞an(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>2x2N+2> 2x^2N+2>2x2 )(详见类指数级数(指数积分函数的变体)数值计算算法的C++实现的第三章)
∣aN+1+p(x)∣=∣x∣2(N+1+p)+1(N+1+p)!=∣x∣2N+1(N+1)!(x2)p(N+2)⋯(N+1+p)⩽∣x∣2N+1(N+1)!(x2N+2)p<∣aN+1(x)∣×(12)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∣aN+1(x)∣⩾2∣aN+1(N+22)∣=2(N+22)2N+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)
无穷远点展开
渐近展开
erf(x)=2π∫0xexp(−t2)dt=2π(π2−exp(−x2)∑n=0∞(−1)n(2n−1)!!2n+11x2n+1)=1−2πexp(−x2)∑n=0∞(−1)n(2n−1)!!2n+11x2n+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=4n=4n=4 (即 x−9x^{-9}x−9 )项时效果最好。因此实际渐近展开式为
erf(x)≈1−2πexp(−x2)[12x−14x3+38x5−1516x7+10532x9](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)
其误差近似小于 10−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⩾log10(1/ε)4+1x\geqslant \frac{\log_{10}\left(1/\varepsilon\right)}{4}+1x⩾4log10(1/ε)+1 时使用无穷渐近展开,反之使用原点渐近展开。
结语
本节采用泰勒展开的方式推导了 erf
函数的展开式,这一展开式在 xxx 较小和较大时计算速度都较快,但 2<x<3.52<x<3.52<x<3.5 时需要迭代次数较多,最坏情况要迭代接近 50
次。更好的计算方法应当为全局拟合,采取预计算常数硬编码的方式,但这不在本文讨论范围内,感兴趣的读者可以参考这篇文章https://www.zhihu.com/answer/1686069171。
高斯误差函数erf的数值计算方法(C++实现)相关推荐
- 数值计算方法 线性方程组的数值解法(4)---向量和矩阵范数(norm) 高斯-赛德尔(Gauss-Seidel)迭代、共轭梯度(Conjugate Gradient)迭代
(范数部分matlab有现成函数,若有需要直接参照matlab_norm) 向量范数 设x∈Rn\boldsymbol x\in \boldsymbol R^nx∈Rn则范数||x||满足:∣∣x∣∣ ...
- 北京科技大学 数值计算方法实验代码
前言: 数值计算方法实验可以使用Matlab.C/C++.Python.Java等语言进行编程,考虑到同学期数学实验课程使用Matlab进行,建议提前熟悉Matlab编程(也效率更高). 本文中各实验 ...
- 高等数值计算方法学习笔记第4章第二部分【数值积分(数值微分)】
高等数值计算方法学习笔记第4章第二部分[数值积分(数值微分)] 四.龙贝格求积公式(第三次课) 1.梯形法的递推化 (变步长求积法) 2.龙贝格算法 五.高斯求积公式 1.一般理论(1定义1例题) 2 ...
- 数值计算方法第三章—线性方程组的数值解法知识点总结
线性方程组的数值解法 本文参考书为马东升著<数值计算方法> 高斯消去法 顺序高斯消去法 通过初等变换消去方程组系数矩阵主对角线以下的元素,而使方程组化为等价的上三角形方程组 列主元高斯消去 ...
- 数值计算方法-算法设计及其MATLAB实现
数值计算方法是一种用于解决数学问题的方法,涉及到数值近似.数值逼近.数值积分.数值微分等等.算法设计是数值计算方法的核心,它包括了数值方法的数学原理和计算机实现的算法,能够有效地解决各种数学问题. M ...
- 插值与多项式逼近的数值计算方法——《数值计算方法》
<数值计算方法>系列总目录 第一章 误差序列实验 第二章 非线性方程f(x)=0求根的数值方法 第三章 CAD模型旋转和AX=B的数值方法 第四章 插值与多项式逼近的数值计算方法 第五章 ...
- 数值计算方法 | C语言实现几个数值计算方法(实验报告版)
目录 写在前面 实验一 牛顿插值方法的实现 实验二 龙贝格求积算法的实现 实验三 高斯列主元消去法的实现 实验四 最小二乘方法的实现 写在前面 使用教材:<数值计算方法>黄云清等编著 科学 ...
- 数值计算方法上机c语言编程,数值计算方法上机实验报告.doc-资源下载在线文库www.lddoc.cn...
<数值计算方法>上机实验报告.doc 华 北 电 力 大 学实 验 报 告实验名称 数值计算方法上机实验 课程名称 数值计算方法 专业班级电力实 08 学生姓名李超然学 号20080100 ...
- 数值计算方法--线性方程组的数值解法(3) 追赶法(Thomas),选择主元(Pivoting),求逆(Inversion)
追赶法 追赶法适用于三对角矩阵,形如(XX000XXX000XXX000XXX000XX)\begin{pmatrix}X&X&0&0&0\\X &X& ...
- 高等数值计算方法学习笔记第6章【解线性代数方程组的迭代方法(高维稀疏矩阵)】
高等数值计算方法学习笔记第6章[解线性代数方程组的迭代方法(高维稀疏矩阵)] 一.引言 1.例题(说明迭代法的收敛性研究的重要性) 2.定义(迭代法,迭代法收敛)&解误差 二.基本迭代法 1. ...
最新文章
- GridView 类型公开的所有成员(公共属性、公共方法、私有属性.......)
- Swift 5.0 值得关注的特性:增加 ResultT, E: Error 枚举类型
- 讨论IM软件企业知识—会谈session的概念,附连到IM软件层次图
- Time包详解二-timer和ticket.html
- jpg 神经网络 手势识别_在STM32上跑神经网络做手势识别
- 如何做好性能压测(一)丨压测环境设计和搭建
- [XSY3381] 踢罐子(几何)
- Java内存体系结构(模型),垃圾回收和内存泄漏
- 将xml转为txt_HZ文章转短视频工具v1.0 快速将文章转为短视频 自动配音 配字幕 配图...
- ios 怎么禁止点击子视图的时候不响应父视图的点击事件
- 使用C#中的反射从字符串获取属性值
- SIM868获取LBS位置
- OA产品的技术发展过程及未来趋势
- 3-9 G: LZY的时间转化
- 基于Spring Boot的农家乐点餐系统
- beego/logs模块的使用
- 全网最全的网络安全技术栈内容梳理(持续更新中)
- DICOM医学图像处理:DICOM存储操作之 “多幅JPG图像数据存入DCM文件”
- 华为鸿蒙 HarmonyOS 2.0 手机开发者 Beta 来了,对开发者意味着什么?
- Axure RP9 进度条设置
热门文章
- 数学建模层次分析法例题及答案_【数模】层次分析法 - 全国大学生数学建模竞赛(CUMCM) - 数学建模社区-数学中国...
- 使用maven命令下载依赖jar
- 内网DNS重要使用作用
- Java:使用Java调用打印机进行打印(JPG、PDF和Word三种文件格式)实现
- C++ 二叉树求叶子结点数及输出叶子结点的路径
- camera(二) DVP接口
- 火星坐标系(高德)和84坐标系互换
- java文件复制后是乱码_复制Java源文件到MyEclipse后乱码问题怎么解决?
- 安装kafka+golang操作kafka
- java毕业设计药品管理系统Mybatis+系统+数据库+调试部署