数值方法求积分 详解+模板代码
什么是数值积分
数值积分可以用来求定积分的近似值。对于很多函数来说,我们是可以使用初等函数来表示出其积分的,对于这种函数,只需要求出不定积分然后代入值就能得到定积分了。
可是除此之外还有许多难求的函数和没法使用初等函数表示的函数。当我们想要求出它们的定积分的时候,需要使用数值积分来求解。
在ACM中一些题目需要使用数值积分来求解,以下列出一些求数值积分的方法,由简单到难,而对ACMer来说最重要的是复合Simpson,其精度较高,且可调精度,是乱搞积分几何的利器。
我从这学的:网易公开课 MIT 数值积分
学习的契机是想要A掉这题:ZOJ 3898我的题解
法一·黎曼和
黎曼和是用将区间等长分为n段,然后用矩形去逼近函数,每段的长为Δx。
可以选择每段左侧的函数值作为矩形的高,也可以选择每段右侧的函数值作为矩形的高。
若设n+1个函数值从左至右为x0,x1⋯xn,可得如下公式:
Left Hand Riemann=Right Hand Riemann=Δxf(x0)+Δxf(x1)+⋯+Δxf(xn−1)=Δx∑i=0n−1f(xi)Δxf(x1)+Δxf(x2)+⋯+Δxf(xn)=Δx∑i=1nf(xi)这种逼近比较粗糙,不过能比较好地传达数值积分的概念。
法二·梯形法
黎曼和虽然简单,但是精度堪忧,没法很好地模拟逼近函数,接下来介绍第二种方法。
可以看出用矩形逼近的时候有很多空缺,使用梯形去逼近,就能大大提高精度了,看上去很像了。
沿用之前的xi,由梯形公式我们可以得到如下公式:Trapezoid=Δx(f(x0)+f(x1))2+Δx(f(x1)+f(x2))2+⋯+Δx(f(xn−1)+f(xn))2=Δx(f(x0)2+f(x1)+⋯+f(xn−1)+f(xn)2)=Left Hand Riemann+Right Hand Riemann2可以看出梯形法求出的就是左右黎曼和的平均值。
公式中第一个和最后一个需要除以2其他都能合出来1个。
法三·Simpson公式
前两种都比较简单,但是精度比较差,第三种方法是用抛物线去逼近。
其实我并不知道是怎么计算的,是从公开课上学来的。
使用Simpson公式首先需要n为偶数。
将整个区间分为n2段,每段的底为2Δx,高比较难搞,但是是有公式的:Simpson height=f(xl)+4f(xm)+f(xr)6
于是有公式:
Simpson=Δx3(f(x0)+4f(x1)+2f(x2)+⋯+2f(xn−2)+4f(xn−1)+f(xn))
法四·复合Simpson
Simpson公式的精度其实已经相当不错了,可是相对于ACM中所需的精度仍然有差距。
为了提高精度,我们需要多次重复利用Simpson公式。
我们先定义被积函数为f(x),定义Simpson(l,r)=(r−l)f(l)+4f(l+r2)+f(r)6
这也就是常规的Simpson公式,再定义函数RSimpson为复合Simpson。
当我们要求RSimpson(l,r),先令m=l+r2
当Simpson(l,r)≈Simpson(l,m)+Simpson(m,r)时我们就认为精度够了返回其中一个。
当不满足的时候我们就再分段去求RSimpson(l,m)+RSimpson(m,r)
这样得到的精度就比较高了,而且通过定义≈的范围可以调整精度。
总结公式如下:
RSimpson(l,r)={Simpson(l,r) approximateRSimpson(l,m)+RSimpson(m,r) else
复合Simpson的实现
inline double getAppr(double fl, double fm, double fr, double l, double r) {return (fl+4*fm+fr)*(r-l)/6.0;
}double Simpson(double l, double r, double fl, double fr) {double m = (l+r)/2, lm = (l+m)/2, rm = (r+m)/2;double fm = f(m), flm = f(lm), frm = f(rm);double vlr = getAppr(fl, fm, fr, l, r);double vlm = getAppr(fl, flm, fm, l, m);double vrm = getAppr(fm, frm, fr, m, r);return fabs(vlr-vlm-vrm) < EPS ? vlr : Simpson(l, m, fl, fm)+Simpson(m, r, fm, fr);
}
复合Simpson的实现2(Natureal的代码)
inline double getAppr(double l,double r){return (f(l) + 4.0*f((l+r)/2) + f(r)) * (r - l) / 6.0;
}double Simpson(double l,double r){double sum = getAppr(l,r);double mid = (l+r)/2;double suml = getAppr(l,mid);double sumr = getAppr(mid,r);return (fabs(sum - suml - sumr) < EPS) ? sum : Simpson(l, mid) + Simpson(mid, r);
}
数值方法求积分 详解+模板代码相关推荐
- C++ - 类模板(class template) 详解 及 代码
类模板(class template) 详解 及 代码 本文地址: http://blog.csdn.net/caroline_wendy/article/details/16906827 类模板(c ...
- 详解模板注入漏洞(下)
作者 | 原作者gosecure,翻译整理shan66 来源 | http://gosecure.github.io/ 在上一篇文章中,我们为读者详细介绍了模版注入漏洞的概念,模版引擎的识别方法,以及 ...
- 数学建模_随机森林分类模型详解Python代码
数学建模_随机森林分类模型详解Python代码 随机森林需要调整的参数有: (1) 决策树的个数 (2) 特征属性的个数 (3) 递归次数(即决策树的深度)''' from numpy import ...
- TOPSIS(逼近理想解)算法原理详解与代码实现
写在前面: 个人理解:针对存在多项指标,多个方案的方案评价分析方法,也就是根据已存在的一份数据,判断数据中各个方案的优劣.中心思想是首先确定各项指标的最优理想值(正理想值)和最劣理想值(负理想解),所 ...
- MGN网络详解以及代码分析
MGN网络详解以及代码分析 最近阅读了云从科技最新的关于REID的论文以及相关的博客和代码,算法是基于MGN,关于网络的部分,这里记录一些自己的学习笔记. 以下是我参考的博客和代码的网址 博客: ht ...
- 三维空间刚体运动1:旋转矩阵与变换矩阵(详解加代码示例)
三维空间刚体运动1:旋转矩阵与变换矩阵(详解加代码示例) 1. 点.向量和坐标系 2.坐标系间的欧式变换 2.1 旋转 2.2 平移 3.齐次坐标和变换矩阵 4. 相似.仿射和射影变换 4.1 相似变 ...
- Go-AES算法详解与代码
目录 AES 发展史 概述 轮函数F 字节代换 行移位 列混淆 轮密钥加 密钥编排 AES和DES的不同之处 分组模式CTR AES的Go实现 aes包 cipher包 加密/解密 参考 本篇介绍分组 ...
- 模板方法模式详解附有代码案例分析(包含模板方法模式重构JDBC操作业务代码示例)
模板方法模式 一.模板方法模式的概念和角色 (一).模板方法模式的概念 (二).模板方法模式的角色 二.模板方法模式的应用场景 三. 模板方法模式的代码示例 四.模板方法模式重构JDBC操作业务 五. ...
- 算法 经典的八大排序算法详解和代码实现
算法 经典的八大排序算法详解和代码实现 排序算法的介绍 排序的分类 算法的时间复杂度 时间频度 示例 图表理解时间复杂度的特点 时间复杂度 常见的时间复杂度 空间复杂度 排序算法的时间复杂度 冒泡排序 ...
- asp.net OutputCache 详解的代码
下面的代码内容是关于asp.net OutputCache 详解的代码. <%@ OutputCache Duration="120" VaryByParam="n ...
最新文章
- 前端学习笔记之this——懂不懂由你,反正我是懂了
- 中国大学科技园市场投资规划及需求前景预测报告2022-2028年版
- linux:Too Many Open Files(打开的文件过多)
- ArrayList遍历
- 在WORD中插入带圈的数字的序号
- Linux socket编程(一) 对套接字操作的封装
- jxl.read.biff.BiffException: Unable to recognize OLE stream解决方法
- sharepoint开发流水账--sharepoint弹出窗体
- 如何理解操作系统的不确定性_如何创造可信任的机器学习模型?先要理解不确定性...
- php i方法和get的区别,浅析PHP中的i++与++i的区别及效率
- [2018.10.18 T3] 玩串
- 计算机组成原理袁春风知识点,计算机组成原理袁春风chap.pdf
- 【2016年第4期】国务院批复建立促进大数据 发展部际联席会议制度
- word文档密码破解
- python输入学生姓名_python学生信息管理系统实现代码
- 老领导调岗,你想跟他干,怎么说?
- java作业Scanner收银
- Win10关闭登录面板毛玻璃效果
- Xilinx Zynq-7000 PL端Kintex-7架构可编程逻辑资源,PS端主频可高达1GHz晶振、电源接口和拨码开关
- 北大计算机科学与技术教材,北京大学计算机科学与技术参考书目