数值分析复习:Richardson外推和Romberg算法

2024-04-23 19:44

本文主要是介绍数值分析复习:Richardson外推和Romberg算法,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

文章目录

    • Richardson外推
    • Romberg(龙贝格)算法

本篇文章适合个人复习翻阅,不建议新手入门使用
本专栏:数值分析复习 的前置知识主要有:数学分析、高等代数、泛函分析

本节继续考虑数值积分问题

Richardson外推

命题:复合梯形公式的另一形式
f ∈ C ∞ [ a , b ] f\in C^{\infty}[a,b] fC[a,b],记 I = ∫ a b f ( x ) d x I=\int_a^bf(x)\mathrm{d}x I=abf(x)dx ,将复合梯形公式记为
T ( h ) = h 2 ∑ i = 0 n − 1 [ f ( x i ) + f ( x i + 1 ) ] T(h)=\frac{h}{2}\sum\limits_{i=0}^{n-1}[f(x_i)+f(x_{i+1})] T(h)=2hi=0n1[f(xi)+f(xi+1)]

T ( h ) = I + α 1 h 2 + α 2 h 4 + ⋯ + α l h 2 l + ⋯ T(h)=I+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_lh^{2l}+\cdots T(h)=I+α1h2+α2h4++αlh2l+

其中 α l ( l = 1 , 2 , … ) \alpha_l(l=1,2,\dots) αl(l=1,2,) h h h 无关

证明
x i + 1 2 = x i + x i + 1 2 , i = 0 , 1 , … , n − 1 x_{i+\frac{1}{2}}=\frac{x_i+x_{i+1}}{2},i=0,1,\dots,n-1 xi+21=2xi+xi+1,i=0,1,,n1

考虑 f ( x ) f(x) f(x) x = x i + 1 2 x=x_{i+\frac{1}{2}} x=xi+21 处的Taylor展开公式
f ( x ) = f ( x i + 1 2 ) + f ′ ( x i + 1 2 ) ( x − x i + 1 2 ) + f ′ ′ ( x i + 1 2 ) 2 ! ( x − x i + 1 2 ) 2 + ⋯ f(x)=f(x_{i+\frac{1}{2}})+f'(x_{i+\frac{1}{2}})(x-x_{i+\frac{1}{2}})+\frac{f''(x_{i+\frac{1}{2}})}{2!}(x-x_{i+\frac{1}{2}})^2+\cdots f(x)=f(xi+21)+f(xi+21)(xxi+21)+2!f′′(xi+21)(xxi+21)2+

若对上述 Taylor 公式代入 x = x i , x = x i + 1 x=x_{i},x=x_{i+1} x=xi,x=xi+1,则得
f ( x i + 1 ) = f ( x i + 1 2 ) + f ′ ( x i + 1 2 ) h 2 + f ′ ′ ( x i + 1 2 ) 2 ! ( h 2 ) 2 + ⋯ f(x_{i+1})=f(x_{i+\frac{1}{2}})+f'(x_{i+\frac{1}{2}})\frac{h}{2}+\frac{f''(x_{i+\frac{1}{2}})}{2!}(\frac{h}{2})^2+\cdots f(xi+1)=f(xi+21)+f(xi+21)2h+2!f′′(xi+21)(2h)2+ f ( x i ) = f ( x i + 1 2 ) + f ′ ( x i + 1 2 ) ( − h 2 ) + f ′ ′ ( x i + 1 2 ) 2 ! ( − h 2 ) 2 + ⋯ f(x_i)=f(x_{i+\frac{1}{2}})+f'(x_{i+\frac{1}{2}})(-\frac{h}{2})+\frac{f''(x_{i+\frac{1}{2}})}{2!}(-\frac{h}{2})^2+\cdots f(xi)=f(xi+21)+f(xi+21)(2h)+2!f′′(xi+21)(2h)2+

两式加和,得到
f ( x i ) + f ( x i + 1 ) 2 = f ( x i + 1 2 ) + h 2 8 f ′ ′ ( x i + 1 2 ) + ⋯ \frac{f(x_i)+f(x_{i+1})}{2}=f(x_{i+\frac{1}{2}})+\frac{h^2}{8}f''(x_{i+\frac{1}{2}})+\cdots 2f(xi)+f(xi+1)=f(xi+21)+8h2f′′(xi+21)+

等式两端求和,乘以 h h h 得到
T ( h ) = h ∑ i = 0 n − 1 f ( x i + 1 2 ) + h 3 8 ∑ i = 0 n − 1 f ′ ′ ( x i + 1 2 ) + ⋯ (1) T(h)=h\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{2}})+\frac{h^3}{8}\sum\limits_{i=0}^{n-1}f''(x_{i+\frac{1}{2}})+\cdots\tag 1 T(h)=hi=0n1f(xi+21)+8h3i=0n1f′′(xi+21)+(1)

另一方面,对Taylor公式从 x i x_i xi x i + 1 x_{i+1} xi+1 进行积分,得到
∫ x i x i + 1 f ( x ) d x = h ⋅ f ( x i + 1 2 ) + f ′ ( x i + 1 2 ) 2 [ ( h 2 ) 2 − ( − h 2 ) 2 ] + f ′ ′ ( x i + 1 2 ) 6 [ ( h 2 ) 3 − ( − h 2 ) 3 ] + ⋯ \int_{x_i}^{x_{i+1}}f(x)\mathrm{d}x=h\cdot f(x_{i+\frac{1}{2}})+\frac{f'(x_{i+\frac{1}{2}})}{2}[(\frac{h}{2})^2-(-\frac{h}{2})^2]+\frac{f''(x_{i+\frac{1}{2}})}{6}[(\frac{h}{2})^3-(-\frac{h}{2})^3]+\cdots xixi+1f(x)dx=hf(xi+21)+2f(xi+21)[(2h)2(2h)2]+6f′′(xi+21)[(2h)3(2h)3]+

等式两端求和得

I = ∑ i = 0 n − 1 ∫ x i x i + 1 f ( x ) d x = h ∑ i = 0 n − 1 f ( x i + 1 2 ) + h 3 24 ∑ i = 0 n − 1 f ′ ′ ( x i + 1 2 ) + ⋯ (2) I=\sum\limits_{i=0}^{n-1}\int_{x_i}^{x_{i+1}}f(x)\mathrm{d}x=h\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{2}}) +\frac{h^3}{24}\sum\limits_{i=0}^{n-1}f''(x_{i+\frac{1}{2}}) +\cdots\tag 2 I=i=0n1xixi+1f(x)dx=hi=0n1f(xi+21)+24h3i=0n1f′′(xi+21)+(2)

结合(1)(2)式,可得
T ( h ) = I + h 3 12 ∑ i = 0 n − 1 f ′ ′ ( x i + 1 2 ) + ⋯ (3) T(h)=I+\frac{h^3}{12}\sum\limits_{i=0}^{n-1}f''(x_{i+\frac{1}{2}})+\cdots\tag 3 T(h)=I+12h3i=0n1f′′(xi+21)+(3)

类似(2)式的推导,可得
∫ a b f ′ ′ ( x ) d x = h ∑ i = 0 n − 1 f ′ ′ ( x i + 1 2 ) + h 3 24 ∑ i = 0 n − 1 f ( 4 ) ( x i + 1 2 ) + ⋯ \int_a^bf''(x)\mathrm{d}x=h\sum\limits_{i=0}^{n-1}f''(x_{i+\frac{1}{2}})+\frac{h^3}{24}\sum\limits_{i=0}^{n-1}f^{(4)}(x_{i+\frac{1}{2}})+\cdots abf′′(x)dx=hi=0n1f′′(xi+21)+24h3i=0n1f(4)(xi+21)+

结合 ∫ a b f ′ ′ ( x ) d x = f ′ ( b ) − f ′ ( a ) \int_a^bf''(x)\mathrm{d}x=f'(b)-f'(a) abf′′(x)dx=f(b)f(a),可将(3)式化为
T ( h ) = I + α 1 h 2 + h 5 c 4 ∑ i = 0 n − 1 f ( 4 ) ( x i + 1 2 ) + ⋯ T(h)=I+\alpha_1h^2+h^5c_4\sum\limits_{i=0}^{n-1}f^{(4)}(x_{i+\frac{1}{2}})+\cdots T(h)=I+α1h2+h5c4i=0n1f(4)(xi+21)+

重复上述操作,考虑 ∫ a b f ( 4 ) ( x ) d x \int_a^bf^{(4)}(x)\mathrm{d}x abf(4)(x)dx,消去 h 5 h^5 h5 的项,得到 h 4 h^4 h4 的项,继续重复操作,可得
T ( h ) = I + α 1 h 2 + α 2 h 4 + ⋯ + α l h 2 l + ⋯ T(h)=I+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_lh^{2l}+\cdots T(h)=I+α1h2+α2h4++αlh2l+

定义:Richardson外推
从低阶精度格式的截断误差的渐近展开式出发,做简单线性计算从而得到高阶精度格式的方法称为Richardson外推

例:
考虑复合梯形公式 T ( h ) T(h) T(h) 满足的式子
T ( h ) = I + α 1 h 2 + α 2 h 4 + ⋯ + α l h 2 l + ⋯ T(h)=I+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_lh^{2l}+\cdots T(h)=I+α1h2+α2h4++αlh2l+

此时截断误差量级为 O ( h 2 ) O(h^{2}) O(h2)

取步长为 h 2 \frac{h}{2} 2h,则有
T ( h 2 ) = I + α 1 h 2 4 + α 2 h 4 16 + ⋯ + α l h 2 l 2 2 l + ⋯ T(\frac{h}{2})=I+\alpha_1\frac{h^2}{4}+\alpha_2\frac{h^4}{16}+\cdots+\alpha_l\frac{h^{2l}}{2^{2l}}+\cdots T(2h)=I+α14h2+α216h4++αl22lh2l+

结合这两个式子,消去 h 2 h^{2} h2项,得
4 T ( h 2 ) − T ( h ) 3 = I − 1 4 α 2 h 4 + ⋯ + α l 3 ( 1 2 2 l − 1 ) h 2 l + ⋯ \frac{4T(\frac{h}{2})-T(h)}{3}=I-\frac{1}{4}\alpha_2h^4+\cdots+\frac{\alpha_l}{3}(\frac{1}{2^{2l}}-1)h^{2l}+\cdots 34T(2h)T(h)=I41α2h4++3αl(22l11)h2l+
T 1 ( h ) = 4 T ( h 2 ) − T ( h ) 3 T_1(h)=\frac{4T(\frac{h}{2})-T(h)}{3} T1(h)=34T(2h)T(h),且
T 1 ( h ) = I + β 2 h 4 + β 3 h 6 + ⋯ + β l h 2 l + ⋯ T_1(h)=I+\beta_2h^4+\beta_3h^6+\cdots+\beta^lh^{2l}+\cdots T1(h)=I+β2h4+β3h6++βlh2l+
若用 T 1 ( h ) T_1(h) T1(h) 估计 I I I ,则截断误差量级提高到 O ( h 4 ) O(h^{4}) O(h4)
类似地,可继续做……

注:只要截断误差可表示为 h h h 的幂级数,均可使用 Richardson外推提高精度

Romberg(龙贝格)算法

在上述对复合梯形公式的截断误差进行Richardson外推的过程中,记复合梯形公式 T 0 ( h ) = T ( h ) = h 2 ∑ i = 0 n − 1 [ f ( x i ) + f ( x i + 1 ) ] T_0(h)=T(h)=\frac{h}{2}\sum\limits_{i=0}^{n-1}[f(x_i)+f(x_{i+1})] T0(h)=T(h)=2hi=0n1[f(xi)+f(xi+1)]

加速一次(即进行一次Richardson外推)后的估计式记为
T 1 ( h ) = 4 T ( h 2 ) − T ( h ) 3 T_1(h)=\frac{4T(\frac{h}{2})-T(h)}{3} T1(h)=34T(2h)T(h)

记加速 n n n 次的估计式为 T n ( h ) T_n(h) Tn(h),则有递推式
T n ( h ) = 4 n 4 n − 1 T n − 1 ( h 2 ) − 1 4 n − 1 T n − 1 ( h ) T_n(h)=\frac{4^n}{4^n-1}T_{n-1}(\frac{h}{2})-\frac{1}{4^n-1}T_{n-1}(h) Tn(h)=4n14nTn1(2h)4n11Tn1(h)

若记 T m ( k ) = T m ( h 2 k ) , k = 0 , 1 , 2 , … T_m^{(k)}=T_m(\frac{h}{2^k}),k=0,1,2,\dots Tm(k)=Tm(2kh),k=0,1,2,,则有递推式
T n ( k ) = 4 n 4 n − 1 T n − 1 ( k + 1 ) − 1 4 n − 1 T n − 1 ( k ) T_n^{(k)}=\frac{4^n}{4^n-1}T_{n-1}^{(k+1)}-\frac{1}{4^n-1}T_{n-1}^{(k)} Tn(k)=4n14nTn1(k+1)4n11Tn1(k)

定理:
设被积函数 f ( x ) f(x) f(x) 充分光滑

  1. lim ⁡ k → ∞ T m ( k ) = I \lim\limits_{k\to\infty}T_m^{(k)}=I klimTm(k)=I
  2. lim ⁡ m → ∞ T m ( k ) = I \lim\limits_{m\to\infty}T_m^{(k)}=I mlimTm(k)=I

注:证明略去,第一个结论说明当节点数目无穷多时, T m ( k ) T_m^{(k)} Tm(k) 收敛于准确的积分值;第二个结论说明随着Richardson外推的进行, T m ( k ) T_m^{(k)} Tm(k) 也收敛于准确的积分值

上述递推式和收敛定理给出了如下的Romberg算法

定义:Romberg算法
对预先给定的精度 ε \varepsilon ε,求 I = ∫ a b f ( x ) d x I=\int_a^bf(x)\mathrm{d}x I=abf(x)dx 的近似值,算法如下:

初始取 k = 0 , m = 0 , h = b − a k=0,m=0,h=b-a k=0,m=0,h=ba

  1. 代入梯形公式,求 T 0 ( k ) ( k = 0 , 1 , 2 , … ) T_0^{(k)}(k=0,1,2,\dots) T0(k)(k=0,1,2,)
  2. 加速一次,由递推公式求 T 1 ( k ) T_1^{(k)} T1(k)
  3. 直至 ∣ T k ( 0 ) − T k − 1 ( 0 ) ∣ < ε |T_k^{(0)}-T_{k-1}^{(0)}|<\varepsilon Tk(0)Tk1(0)<ε,则取 T k ( 0 ) ≈ I T_{k}^{(0)}\approx I Tk(0)I

注:具体求解顺序如下表

在这里插入图片描述

参考书籍:《数值分析》李庆扬 王能超 易大义 编

这篇关于数值分析复习:Richardson外推和Romberg算法的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/929709

相关文章

不懂推荐算法也能设计推荐系统

本文以商业化应用推荐为例,告诉我们不懂推荐算法的产品,也能从产品侧出发, 设计出一款不错的推荐系统。 相信很多新手产品,看到算法二字,多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法(注:传统是指科班出身的产品都会接触过)。但对于推荐算法,多数产品对着网上搜到的资源,都会无从下手。特别当某些推荐算法 和 “AI”扯上关系后,更是加大了理解的难度。 但,不了解推荐算法,就无法做推荐系

康拓展开(hash算法中会用到)

康拓展开是一个全排列到一个自然数的双射(也就是某个全排列与某个自然数一一对应) 公式: X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! 其中,a[i]为整数,并且0<=a[i]<i,1<=i<=n。(a[i]在不同应用中的含义不同); 典型应用: 计算当前排列在所有由小到大全排列中的顺序,也就是说求当前排列是第

性能分析之MySQL索引实战案例

文章目录 一、前言二、准备三、MySQL索引优化四、MySQL 索引知识回顾五、总结 一、前言 在上一讲性能工具之 JProfiler 简单登录案例分析实战中已经发现SQL没有建立索引问题,本文将一起从代码层去分析为什么没有建立索引? 开源ERP项目地址:https://gitee.com/jishenghua/JSH_ERP 二、准备 打开IDEA找到登录请求资源路径位置

csu 1446 Problem J Modified LCS (扩展欧几里得算法的简单应用)

这是一道扩展欧几里得算法的简单应用题,这题是在湖南多校训练赛中队友ac的一道题,在比赛之后请教了队友,然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意:把题意转变下就变成了:求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点:求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个

综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

LntonAIServer视频质量诊断功能中的抖动检测是一个专门针对视频稳定性进行分析的功能。抖动通常是指视频帧之间的不必要运动,这种运动可能是由于摄像机的移动、传输中的错误或编解码问题导致的。抖动检测对于确保视频内容的平滑性和观看体验至关重要。 优势 1. 提高图像质量 - 清晰度提升:减少抖动,提高图像的清晰度和细节表现力,使得监控画面更加真实可信。 - 细节增强:在低光条件下,抖

【数据结构】——原来排序算法搞懂这些就行,轻松拿捏

前言:快速排序的实现最重要的是找基准值,下面让我们来了解如何实现找基准值 基准值的注释:在快排的过程中,每一次我们要取一个元素作为枢纽值,以这个数字来将序列划分为两部分。 在此我们采用三数取中法,也就是取左端、中间、右端三个数,然后进行排序,将中间数作为枢纽值。 快速排序实现主框架: //快速排序 void QuickSort(int* arr, int left, int rig

poj 3974 and hdu 3068 最长回文串的O(n)解法(Manacher算法)

求一段字符串中的最长回文串。 因为数据量比较大,用原来的O(n^2)会爆。 小白上的O(n^2)解法代码:TLE啦~ #include<stdio.h>#include<string.h>const int Maxn = 1000000;char s[Maxn];int main(){char e[] = {"END"};while(scanf("%s", s) != EO

秋招最新大模型算法面试,熬夜都要肝完它

💥大家在面试大模型LLM这个板块的时候,不知道面试完会不会复盘、总结,做笔记的习惯,这份大模型算法岗面试八股笔记也帮助不少人拿到过offer ✨对于面试大模型算法工程师会有一定的帮助,都附有完整答案,熬夜也要看完,祝大家一臂之力 这份《大模型算法工程师面试题》已经上传CSDN,还有完整版的大模型 AI 学习资料,朋友们如果需要可以微信扫描下方CSDN官方认证二维码免费领取【保证100%免费

dp算法练习题【8】

不同二叉搜索树 96. 不同的二叉搜索树 给你一个整数 n ,求恰由 n 个节点组成且节点值从 1 到 n 互不相同的 二叉搜索树 有多少种?返回满足题意的二叉搜索树的种数。 示例 1: 输入:n = 3输出:5 示例 2: 输入:n = 1输出:1 class Solution {public int numTrees(int n) {int[] dp = new int

SWAP作物生长模型安装教程、数据制备、敏感性分析、气候变化影响、R模型敏感性分析与贝叶斯优化、Fortran源代码分析、气候数据降尺度与变化影响分析

查看原文>>>全流程SWAP农业模型数据制备、敏感性分析及气候变化影响实践技术应用 SWAP模型是由荷兰瓦赫宁根大学开发的先进农作物模型,它综合考虑了土壤-水分-大气以及植被间的相互作用;是一种描述作物生长过程的一种机理性作物生长模型。它不但运用Richard方程,使其能够精确的模拟土壤中水分的运动,而且耦合了WOFOST作物模型使作物的生长描述更为科学。 本文让更多的科研人员和农业工作者