惯性聚合 高效追踪和阅读你感兴趣的博客、新闻、科技资讯
阅读原文 在惯性聚合中打开

推荐订阅源

美团技术团队
Blog — PlanetScale
Blog — PlanetScale
阮一峰的网络日志
阮一峰的网络日志
M
MIT News - Artificial intelligence
月光博客
月光博客
钛媒体:引领未来商业与生活新知
钛媒体:引领未来商业与生活新知
U
Unit 42
博客园_首页
WordPress大学
WordPress大学
H
Hackread – Cybersecurity News, Data Breaches, AI and More
J
Java Code Geeks
F
Fortinet All Blogs
腾讯CDC
罗磊的独立博客
IT之家
IT之家
I
InfoQ
V
V2EX
博客园 - 叶小钗
A
About on SuperTechFans
Y
Y Combinator Blog
C
Check Point Blog
量子位
Martin Fowler
Martin Fowler
Vercel News
Vercel News

你是下雨天

老派约会之必要 AI的逗号 相机配件 烘焙初体验——蛋挞和泡芙 摩托驾驶证科目一易错点 摩托驾驶证科目一易错点 三维旋转、欧拉角、四元数 维生素每日摄入量 高斯分布常用公式 爱得太迟 《娱乐至死》读后感 科目三 无线通信标准杂谈 信号调制与解调 IEEE论文爬虫及数据统计 SLAM综述 我的电脑(不时更新) 【高等数值分析】Krylov子空间方法 【随机过程2】连续参数马尔可夫链 【随机过程2】马尔可夫过程2 | 状态空间 【随机过程2】马尔可夫过程1 | 基本概念 【随机过程2】布朗运动2 | 推广 Hello World 【高等数值分析】常微分方程数值解
【高等数值分析】数值积分和数值微分
Glooow · 2022-01-09 · via 你是下雨天

1. 预备理论

根据Newton-Leibniz公式有 \(\int_a^x f(t)dt = F(x)-F(a)\),但是绝大部分情况很难解析求解,需要数值积分。例如中点公式 \[ \int_a^b f(x)dx \approx f(\frac{a+b}{2})(b-a) \]\(f(x)\in C^2[a,b]\),则中点公式截断误差\[ \int_a^b f(x)dx - f(\frac{a+b}{2})(b-a) = \frac{(b-a)^2}{24} f''(\xi), \quad \xi\in(a,b) \]

2. 插值型求积公式

顾名思义,就是先插值再求积分。方法为给定求积节点 \(x_k,k=0,1,...,n\) 和求积系数 \(A_k,k=0,1,...,n\),插值型求积公式表示为 \[ \int_a^b f(x)dx \approx \sum_{k=0}^n A_k f(x_k) \] 根据多项式插值中的理论,余项可表示为 \[ E_n(f) = \frac{1}{(n+1)!}\int_a^b f^{(n+1)}(\xi(x)) w_{n+1}(x) dx \]\(f(x)\in{\mathcal P}_n\) 时都有 \(E_n(f)=0\)。由此可以定义代数精度,如果对所有 \(p\in{\mathcal P}_m\)\(E_n(p)=0\),而对某个 \(q\in{\mathcal P}_{m+1}\)\(E_n(q)\ne 0\),称求积公式具有 \(m\)代数精度

插值型求积公式主要分为两类:Newton-Cotes求积公式和Gauss型求积公式。前者等距选取插值节点,后者则未必。

2.1 Newton-Cotes求积公式

方法是将区间 \([a,b]\) \(n\) 等分,得到 \(h = \frac{b-a}{n},x_k=a+kh,k=0,...,n\),再利用Lagrange插值公式 \(l_k(x)\),得到 \[ \int_a^b f(x)dx = (b-a)\sum_{k} \frac{f(x_k)}{b-a}\int_a^b l_k(x)dx = (b-a)\sum_k C_k^{(n)} f(x_k) \] 其中 \(C_k^{(n)}=\frac{1}{b-a}\int_a^b \prod_{j\ne k}\frac{x-x_j}{x_k-x_j}dx = \frac{1}{n}\int_a^b \prod_{j\ne k}\frac{t-j}{k-j}dx\) 称为 Cotes 求积系数,不仅与被积函数无关与求积区间也无关

该方法有如下性质:

  • \(\sum_{k=0}^n C_k^{(n)} = 1\)(取 \(f\equiv1\) 即可得证);
  • \(E_n(f) = \frac{1}{(n+1)!}\int_a^b f^{(n+1)}(\xi(x)) w_{n+1}(x) dx\)
  • \(n\) 为偶数时,代数精度为 \(n+1\);当 \(n\) 为奇数时,代数精度为 \(n\)。(直观理解是因为奇次多项式的奇对称性积分后恰好为 0)

下面是 \(n\) 取不同数值的特例。

\(n=1\) 时为梯形公式,代数精度为 \(1\)\[ \begin{align} &\int_a^b f(x)dx \approx \frac{b-a}{2}[f(a)+f(b)] \\ &E_1(f) = -\frac{(b-a)^3}{12} f''(\xi), \xi\in[a,b] \end{align} \] \(n = 2\) 时为 Simpson 公式,代数精度为 \(3\)\[ \begin{align} &\int_a^b f(x)dx \approx \frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)] \\ &E_1(f) = -\frac{(b-a)^5}{2880} f^{(4)}(\xi), \xi\in[a,b] \end{align} \]

2.2 Gauss型求积公式

求积公式也表示为 \[ \int_a^b f(x)dx \approx \sum_{k=0}^n A_k f(x_k) \] 但与Newton-Cotes方法不同的是求积节点并非等距选取,而是将参数待定,解出使代数精度最高的参数 \(A_k\)\(x_k\)。共有 \(2n+2\) 个参数待定,分别取 \(f(x)=1,x,x^2,...,x^{2n+1}\) 列方程组,因此代数精度最高可以达到 \(2n+1\)

从理论上也可以证明代数精度至多为 \(2n+1\)

证明:取 \(f(x)=\prod_{k=0}^n(x-x_k)^2\),那么上面等式左边一定有 \(\int_a^b f(x)dx > 0\),但是等式右边有 \(\sum_{k} A_k f(x_k)=0\),显然不相等。证毕。

除此之外,Gauss型求积公式还可以拓展到带权积分的情况,也即 \[ \int_a^b \rho(x)f(x)dx = \sum_{k=0}^n A_k f(x_k) + \frac{1}{(n+1)!}\int_a^b f^{(n+1)}(\xi(x)) w_{n+1}(x) dx \] 关键问题是如何求解 \(A_k,x_k\) 呢?

由于代数精度为 \(2n+1\),对 \(\forall f\in{\mathcal P}_{2n+1}\),截断误差应满足 \(\frac{1}{(n+1)!}\int_a^b f^{(n+1)}(\xi(x)) w_{n+1}(x) dx=0\)。又由于 \(f^{(n+1)}\in{\mathcal P}_n\),那么只需要取 \(w_{n+1}(x)\)\(n+1\) 阶正交多项式(权函数为 \(\rho\))即可满足该条件。因此Gauss型求积公式的求解方法为:

  1. 权函数为 \(\rho\),求 \(n+1\) 阶正交多项式;
  2. \(n+1\) 个根 \(x_k\)(在逼近一章中已经证明了一定存在 \(n+1\) 个不同的根);
  3. \(f(x)=1,x,...,x^n\) 列线性方程组求 \(A_k\)

Gauss型求积公式有如下性质:

  • \(A_k > 0,k=0,1,...,n\)(取 \(f(x)=l_i(x)\) 即可证明);
  • \(\sum_{k=0}^n A_k = \int_a^b \rho(x)dx\)(取 \(f\equiv\) 即可证明);
  • 余项 \(E_n(f) = \frac{1}{(2n+2)!}f^{(2n+2)}(\eta)\int_a^b \rho(x) [w_{n+1}(x)]^2 dx\)

3. 数值微分

数值微分和数值积分类似,也有插值型微分公式,即先进性多项式插值,再求导数值。

还可以利用Richardson外推公式提高精度。