luoguP4526 【模板】自适应辛普森法2

guapi2333

于 2021-01-13 21:43:11 发布

阅读量163

点赞数

CC 4.0 BY-SA版权

分类专栏：计算几何

本文链接：https://blog.csdn.net/guapi2333/article/details/112593993

计算几何专栏收录该内容

4 篇文章

订阅专栏

博客围绕积分∫0∞xxa−xdx(∣a∣≤50)展开，先分析给定积分，因无法用常规方法求原函数，故考虑审敛法，介绍了无穷界积分和无穷间断点的审敛法，判断积分在不同情况下的收敛性，还给出a≥0时的代码实现及总结。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Description

试计算积分 $∫0∞xax−xdx(∣a∣≤50)\int_{0}^{\infty}x^{\frac{a}{x}-x}dx(|a|\le50)$ ，保留至小数点后 $5$ 位。若积分发散，请输出orz。

Solution

如何分析题目中给定积分

设题目给定的无界（广义）积分中被积函数为 $f (x)$ ；首先，上述积分中被积函数自变量 $x$ 取不到 $0$ 且无上界；其次，由于 $f (x)$ 是“幂指函数”，故先“指对互化”：

$f(x)=xax−x=e(ax−x)lnxf(x)=x^{\frac{a}{x}-x}=e^{(\frac{a}{x}-x)lnx}$ 设 $g(x)=(ax−x)lnxg(x)=(\frac{a}{x}-x)lnx$ ，则 $f(x)=e^{g(x)}$

通常情况下，对于此题，我们判断 $x→0x\to0$ 这种情况下积分是否收敛可试着求 $lim⁡ε→0+∫0+ε∞xax−xdx\lim_{\varepsilon\to0^+}\int_{0+\varepsilon}^{\infty}x^{\frac{a}{x}-x}dx$ 看其是否存在（即是否等于一常数）；对于 $x→+∞x\to+\infty$ 的情况，则看 $f (x)$ 原函数在 $x→+∞x\to+\infty$ 时极限是否存在；但此积分着实没法用常规方法求原函数，故考虑审敛法。

（注意：被积函数在某处趋近于正无穷（常数）并不一定代表其积分发散（收敛）！即：不能简单根据被积函数的取值判断收敛发散。

Eg：1、 $∫01lnxdx\int_{0}^{1}lnxdx$ ：此积分可用上面设 $ε\varepsilon$ 的方法分部积分后得到其 $= - 1$ ，收敛；2、虽然 $x→+∞x\to+\infty$ 时 $1x\frac{1}{x}$ 趋近于 $0$ ，但 $∫1+∞1xdx\int_{1}^{+\infty}\frac{1}{x}dx$ 显然发散。）

关于审敛法

一、无穷界积分的审敛法

设 $f (x)$ 在积分区间 $[a,+∞)[a,+\infty)$ 上恒有 $f(x)≥0f(x)\ge0$ ：

（1）、若存在常数 $p > 1$ 使得 $lim⁡x→+∞xpf(x)=c\lim_{x\to+\infty}x^pf(x)=c$ ，那么 $∫a+∞f(x)dx\int_{a}^{+\infty}f(x)dx$ 收敛。

（2）、若 $lim⁡x→+∞xf(x)=d>0\lim_{x\to+\infty}xf(x)=d>0$ ，那么 $∫a+∞f(x)dx\int_{a}^{+\infty}f(x)dx$ 发散。

将积分区间换为 $(−∞,a](-\infty,a]$ 时定理仍成立。

二、无穷间断点的审敛法

设 $a$ 是 $f (x)$ 的无穷间断点， $f (x)$ 在 $(a, b]$ 上有 $f(x)≥0f(x)\ge0$ ：

1、若存在常数 $0 < p < 1$ 使得 $lim⁡x→a+(x−a)pf(x)=c\lim_{x\to a^+}(x-a)^pf(x)=c$ ，那么 $∫abf(x)dx\int_{a}^{b}f(x)dx$ 收敛。

2、若 $lim⁡x→a+(x−a)f(x)=d>0\lim_{x\to a^+}(x-a)f(x)=d>0$ ，那么 $∫abf(x)dx\int_{a}^{b}f(x)dx$ 发散。

类似的，设 $b$ 是 $f (x)$ 的无穷间断点， $f (x)$ 在 $[a, b)$ 上有 $f(x)≥0f(x)\ge0$ ：

1、若存在常数 $0 < p < 1$ 使得 $lim⁡x→b−(b−x)pf(x)=c\lim_{x\to b^-}(b-x)^pf(x)=c$ ，那么 $∫abf(x)dx\int_{a}^{b}f(x)dx$ 收敛。

2、若 $lim⁡x→b−(b−x)f(x)=d>0\lim_{x\to b^-}(b-x)f(x)=d>0$ ，那么 $∫abf(x)dx\int_{a}^{b}f(x)dx$ 发散。

（注意：以上 $c$ 均为常数， $d$ 可为常数也可取无穷）

上述两种审敛法的证明回头拿到同济版高数课本再补。

对于此题： $x≥0x\ge0$ ，故 $f(x)≥0f(x)\ge0$ 恒成立。

1、 $x→0+x\to0^+$ ： $lim⁡x→0+f(x)=elim⁡x→0+(a−x2)lnxx=elim⁡x→0+(a−x2)lim⁡x→0+lnxx\lim_{x\to0^+}f(x)=e^{\lim_{x\to 0^+}\frac{(a-x^2)lnx}{x}}=e^{\lim_{x\to 0^+}(a-x^2)\lim_{x\to 0^+}\frac{lnx}{x}}$

$=ealim⁡x→0+lnxx==e^{a\lim_{x\to 0^+}\frac{lnx}{x}}=$ (1) $e+∞=+∞,a<0e^{+\infty}=+\infty,a<0$

(2) $e^{0}=1,a=0$

(3) $e−∞=0,a>0e^{-\infty}=0,a>0$

所以只有 $a < 0$ 时 $∫0∞f(x)dx\int_{0}^{\infty}f(x)dx$ 可能在 $x = 0$ 处发散。

又 $∵\because$ $lim⁡x→0+(x−0)f(x)=elim⁡x→0+(a−x2+x)lnxx=elim⁡x→0+lnx(1−x+ax)=e+∞=+∞\lim_{x\to 0^+}(x-0)f(x)=e^{\lim_{x\to 0^+}\frac{(a-x^2+x)lnx}{x}}=e^{\lim_{x\to 0^+}lnx(1-x+\frac{a}{x})}=e^{+\infty}=+\infty$ ，故 $∫0∞f(x)dx\int_{0}^{\infty}f(x)dx$ 在 $x = 0$ 处发散。

2、 $x→+∞x\to+\infty$ ：可利用洛必达法则求得：

$lim⁡x→+∞x2f(x)=elim⁡x→+∞(a−x2+2x)lnxx=elim⁡x→+∞[(−2x+2)lnx+ax−x+2]=e−∞=0\lim_{x\to+\infty}x^2f(x)=e^{\lim_{x\to+\infty}\frac{(a-x^2+2x)lnx}{x}}=e^{\lim_{x\to+\infty}[(-2x+2)lnx+\frac{a}{x}-x+2]}=e^{-\infty}=0$ ，故 $∫0∞f(x)dx\int_{0}^{\infty}f(x)dx$ 在 $x→+∞x\to +\infty$ 处收敛。

$a≥0a\ge0$ 情况代码实现

由于 $a > 0$ 时原积分在 $x→+∞x\to +\infty$ 处收敛且 $f(10)=10a10−10,f(20)=20a20−20f(10)=10^{\frac{a}{10}-10},f(20)=20^{\frac{a}{20}-20}$ 。即便 $a$ 取到题目给定上限 $50$ ， $f(10)=10^{-5}=0.0001,f(20)=20^{-15}$ 小的可以在保留 $5$ 位小数情况下忽略不记，故自适应 $S i m p s o n$ 法的左右端点可取 $0, 20$ （注意：右端点取值过大会导致初始状态用 $S i m p s o n$ 公式近似定积分值时得到 $f (m i d)$ 与 $f (r)$ 等于 $0$ 从而导致拟合误差较大）。同时，为了提高准确性，初始精度 $e p s$ 可适当调高。

总结

1、若写的自适应Simpson法实际运行的时候出了误差，可从修改初次调用Simpson函数时积分区间的端点 $l$ 和 $r$ 、长度 $e p s$ 的值入手；

2、求极限时：能求提则求提！！（再次重申：求提原理：若 $f (a) 、 f (b)$ 在 $x→ax\to a$ 时极限均存在，则 $limx→af(a)f(b)=limx→af(a)limx→af(b)lim_{x\to a}f(a)f(b)=lim_{x\to a}f(a)lim_{x\to a}f(b)$ ！）；

3、公式推导务求过程严谨。
Code

#include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;

double A;

inline double dcmp(double x,double y) { return fabs(x-y)<=1e-6; }
inline double f(double x) 
{ 
	if(dcmp(x,0.0))
	{
		if(dcmp(A,0.0))	return 1.0;
		return 0.0;		
	}
	return exp((A-x*x)*log(x)/x); 
}
inline double Simpson(double a,double b) { return (f(a)+4*f((a+b)/2.0)+f(b))*(b-a)/6.0; }

inline double queryans(double l,double r,double eps,double ans)
{
	register double mid=(l+r)/2.0;
	double ansl=Simpson(l,mid),ansr=Simpson(mid,r);
	if(fabs(ansl+ansr-ans)<=15*eps)	return ansl+ansr+(ansl+ansr-ans)/15.0;
	return queryans(l,mid,eps/2.0,ansl)+queryans(mid,r,eps/2.0,ansr);
}

int main()
{
	scanf("%lf",&A);
	if(A<0.0)	{ cout<<"orz"; return 0; }
	printf("%.5lf",queryans(0,20,1e-7,Simpson(0,20)));
	return 0;
}