Romberg 数值积分算法+P3779 题解

Romberg 算法

吊打 Simpson（？） 的且不玄学（没有什么十五倍）的数值积分算法。

缺点是过程复杂一点，但是只体现在证明上，代码很短。

铺垫算法

梯形求积公式

公式

\[\int _a^b f(x)dx\approx \frac{(f(a)+f(b))(b-a)}2\\ \text{令 }(1)=\frac{(f(a)+f(b))(b-a)}2 \]

计算梯形求积公式的误差

注意到

\[\int _a^bf(x)dx=\int_a^b 1\cdot f(x)d(x-a)\\ =(b-a)f(b)-\int_a^b (x-a)f'(x)dx\\ =(b-a)f(a)-\int_a^b (x-b)f'(x)dx \]

上两式相加得到

\[\int_a^b f(x)dx=\frac{(f(a)+f(b))(b-a)}2-\frac{1}{2}\int_a^b((x-a)+(x-b))f'(x)dx\\ =(1)-\frac{1}{2}\int_a^b((x-a)+(x-b))f'(x)dx\\ =(1)+\frac 12\int _a^b(x-a)(x-b)f''(x)dx \]

根据定积分第一中值定理：

\[\exists \xi \in [a,b],\int_a^bf(x)g(x)=f(\xi)\int_a^bg(x)dx \]

\[\int_a^bf(x)dx=(1)+\frac{f(\xi)}2\int_a^b(x-a)(x-b)dx\\ =(1)-\frac{f''(\xi)(b-a)^3}{12} \]

误差为

\[-\frac{f''(\xi)(b-a)^3}{12} \]

复化梯形求积公式

把他分割成 \(n\) 段，每段应用梯形求积公式就得到了复化梯形求积公式。

公式

设 \(x_k=a+\dfrac{k(b-a)}n\)。

\[\int _a^bf(x)dx=\sum _{k=0}^{n-1}\int _{x_k}^{x_{k+1}}f(x)dx\approx \frac{b-a}{2n}(f(a)+f(b)+2\sum _{k=1}^{n-1}f(x_k))\\ \text{令 }(2)=\frac{b-a}{2n}(f(a)+f(b)+2\sum _{k=1}^{n-1}f(x_k)) \]

计算复化梯形求积公式的误差

\[\int_a^b f(x)dx-(2)=\sum _{k=0}^{n-1}\int_{x_k}^{x_{k+1}}f(x)dx-\sum_{k=0}^{n-1} \frac{a-b}{2n}(f(x_k)+f(x_{k+1}))\\ =\sum _{k=0}^{n-1}\int_{x_k}^{x_{k+1}}f(x)dx-\frac{a-b}{2n}(f(x_k)+f(x_{k+1}))\\ =-\frac{1}{12}\sum_{k=0}^{n-1}f''(\xi_k)(x_{k+1}-x_k)^3\\ =-\frac{(b-a)^3}{12n^3}\sum_{k=0}^{n-1}f''(\xi_k)\\ \]

由于 \(\xi_k\in[a,b]\)，那么

\[\exists\ \eta\in[a,b],\sum_{k=0}^{n-1}\frac{f''(\xi_k)}{n}=f''(\eta) \]

再设 \(h=\dfrac{b-a}n\)：

\[\int_a^b f(x)dx-(2)=-\frac{h^2}{12}(b-a)f''(\eta) \]

那么

\[\int_a^b f(x)dx-(2)=\sum _{k\ge 1}\alpha_kh^{2k} \]

略，可参考蒋尔雄等著《数值逼近》。

没找到证明。

理查德森外推法

考虑一步长为 \(h\) 的级数：

\[x=F(h)+ah^p+bh^q+\dots \]

可以采取方法，去掉 \(ah^p\) ，以期在取前几项的计算中得到更高精度。

考虑一个步长 \(h_2\)，

\[x=F(h_2)+ah_2^q+qh_2^q+\dots \]

设 \(r=\dfrac{h^p}{h_2^p}\)。

\[(1-r)x=F(h)-rF(h_2)+b(h^q-rh_2^q)+\dots\\ x=\frac{F(h)-rF(h_2)}{1-r}+\frac{b(1-(\frac{h_2}{h})^{q-p})}{1-r}h^q+\dots \]

新的

\[F^*(h)=\frac{F(h)-rF(h_2)}{1-r} \]

应用在减少复化梯形求积公式的误差上会发挥作用。

Romberg 积分法

每次外推取 \(h_2=2h\)，增加同一 \(n\) 的精度；而不断把 \(n\) 变为原来的 \(2\) 倍，增大 \(n\) 带来的精度。基于这样的想法，定义龙贝格函数：

\[R(n,m)=\left\{\begin{matrix} \dfrac{(b-a)(f(a)+f(b))}2 & (n=m=0)\\ \dfrac{R(n-1,0)}{2}+\dfrac{b-a}{2^n}\displaystyle{\sum_{k=1}^{2^n-1} f(a+\frac{(b-a)(2k-1)}{2^n})} & (n\neq 0,m=0)\\ \dfrac{4^mR(n,m-1)-R(n-1,m-1)}{4^m-1} \end{matrix}\right. \]

然后依照上面的内容，

\[\lim_{(n,m)\to \infty}=\int_a^b f(x)dx \]

实践中，考虑取 \(m\le 3\)，算完之后，令 \(n\gets n+1\) 来提高精度。

下面是 P4526 的代码。

// Problem: P4526 【模板】自适应辛普森法 2
// Platform: Luogu
// URL: https://www.luogu.com.cn/problem/P4526
// Memory Limit: 62 MB
// Time Limit: 500000 ms
// Author:British Union
// Long live UOB and koala
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
using namespace std;
#define db double 
class Romberg{
	protected:
	function<db(db)> f;
	db calc(db a,db b,db h){
		db ans=0.0;
		for(db i=a+h/2.0;i<b;i+=h)ans+=f(i);
		return ans;
	}
	public:
	Romberg(function<db(db)> F=[](db x){return x;}){f=F;}
	db operator ()(db a,db b,db eps){
		db T1,T2,S1,S2,C1,C2,R1,R2,h,S;
		int cnt,k;
		h=b-a,T1=(f(a)+f(b))*(b-a)/2.0,cnt=0,k=1;
		while(1){
			// cout<<cnt<<endl;
			S=calc(a,b,h),T2=(T1+h*S)/2.0;
			S2=(4*T2-T1)/3.0,h=h/2.0;
			T1=T2,S1=S2,C1=C2,R1=R2,cnt++;
			if(k==1){k++;continue;}
			C2=(16*S2-S1)/15.0;
			if(k==2){k++;continue;}
			R2=(64*C2-C1)/63.0;
			if(k==3){k++;continue;}
			if(fabs(R1-R2)<eps||cnt>40)break;
		}
		return R2;
	}
};
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	db a;cin>>a;
	if(a<0){
		cout<<"orz"<<endl;
		return 0;
	}
	printf("%.5lf",Romberg([a](db x){return pow(x,a/x-x);})(1e-100,20,1e-12));
	return 0;
}

P3779 应用一下

令 \(n=Y,\hat X=\sum_{i=1}^n X\)。

希望求的就是 \(P(A\le X\le B)\)。

首先有比较 simple 的 FFT，但是 \(n,X\) 很大时会寄。

考虑中心极限定理。

正态分布函数为

\[P(X\le t)=\frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\infty}^t \exp(-\dfrac{(x-\mu)^2}{2\sigma^2})dx \]

称其满足 \(N(\mu,\sigma^2)\) 的正态分布。

而设满足 \(X\sim N(0,1)\) 的 \(X\) 满足分布 \(P(X\le t)=\Phi(t)\)。

中心极限定理：设 \(n\) 个满足同一分布的随机变量 \(X_1\sim X_n\)。

设 \(\mu= E[X],\sigma^2=V[X]\)。

\[Y=\frac{\sum_{i=1}^n X_i-n\mu}{\sqrt{n\sigma^2}}\\ Y\sim N(0,1) \]

依靠中心极限定理，在 \(X,n\) 较大的时候通过正态分布近似计算。而 \(\sum X_i\) 的范围可以变成 \(Y\) 的范围。

// Problem: P3779 [SDOI2017] 龙与地下城
// Platform: Luogu
// URL: https://www.luogu.com.cn/problem/P3779
// Memory Limit: 500 MB
// Time Limit: 1000 ms
// Author:British Union
// Long live UOB and koala
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define db double
class Romberg{
	protected:
	function<db(db)> f;
	db calc(db a,db b,db h){
		db ans=0;
		for(db i=a+h/2;i<b;i+=h)ans+=f(i);
		return ans;
	}
	public:
	Romberg(function<db(db)> F=[](db x){return x;}){f=F;}
	db operator ()(db a,db b,db eps){
		if(b<a)return 0.0;
		db T1,T2,S1,S2,C1,C2,R1,R2,h,S;
		int cnt,k;
		h=b-a,T1=(f(a)+f(b))*(b-a)/2,cnt=0,k=1;
		while(1){
			S=calc(a,b,h),T2=(T1+h*S)/2;
			S2=(4*T1-T1)/3,h=h/2;
			T1=T2,S1=S2,C1=C2,R1=R2,cnt++;
			if(k==1){k++;continue;}
			C2=(16*S2-S1)/15.0;
			if(k==2){k++;continue;}
			R2=(64*C2-C1)/63.0;
			if(k==3){k++;continue;}
			if(fabs(R1-R2)<eps||cnt>20)break;
		}
		return R2;
	}
};
const int maxn=4e5+5;
const db pi=acos(-1.0);
int lim=1,L=0;
int r[maxn];
#define qsy complex<db>
void predo(int n){
	lim=1,L=0;
	while(lim<=n)lim<<=1,L++;
	for(int i=1;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<L-1);
}
void fft(qsy a[],int tp){
	for(int i=0;i<lim;i++)if(r[i]>i)swap(a[r[i]],a[i]);
	for(int mid=1;mid<lim;mid<<=1){
		qsy w(cos(pi/mid),tp*sin(pi/mid));
		for(int j=0;j<lim;j+=(mid<<1)){
			qsy W(1,0);
			for(int k=0;k<mid;k++,W=W*w){
				qsy x=a[j+k],y=a[j+k+mid]*W;
				a[j+k]=(x+y),a[j+k+mid]=(x-y);
			}
		}
	}
	if(tp==-1)for(int i=0;i<lim;i++)a[i]=a[i]/(1.0*lim);
}
int T,x,n,N=2e5;
db calc(db A){
	A=min(A,50.0);
	return (1/sqrt(2*pi))*Romberg([](db x){return exp(-x*x/2);})(-50,A,0.001);
}
qsy t[maxn],tmp[maxn];
signed main(){
	// ios::sync_with_stdio(0);
	// cin.tie(0);cout.tie(0);
	cin>>T;
	while(T--){
		cin>>x>>n;
		db E=(x-1.0)/2.0,D=sqrt(1.0*n*(x*x-1.0)/12.0);
		if(n<=800){
			N=x*n;predo(N+1);
			for(int i=0;i<lim;i++)t[i]=tmp[i]=qsy(0,0);
			for(int i=0;i<x;i++)t[i]=tmp[i]=qsy(1.0/x,0);
			--n;
			// fft(tmp,1);
			for(int i=0;(1<<i)<=n;i++){
				if(n&(1<<i)){
					fft(tmp,1);fft(t,1);
					for(int j=0;j<lim;j++)t[j]=t[j]*tmp[j];
					fft(t,-1);fft(tmp,-1);
					for(int j=N+1;j<lim;j++)t[j]=tmp[j]=0;
				}
				fft(tmp,1);
				for(int j=0;j<lim;j++)tmp[j]=tmp[j]*tmp[j];
				fft(tmp,-1);
				for(int j=N+1;j<lim;j++)tmp[j]=0;
			}
			for(int i=1;i<=10;i++){
				int A,B;cin>>A>>B;
				db ans=0;
				for(int j=A;j<=B;j++)ans+=t[j].real();
				printf("%.6lf\n",ans);
			}
		}else{
			for(int i=1;i<=10;i++){
				int A,B;cin>>A>>B;
				printf("%.6lf\n",calc((0.0+B-n*E)/D)-calc((0.0+A-n*E)/D));
			}
		}
	}
	return 0;
}