P4525 【模板】自适应辛普森法 1 题解

HYC0907

于 2024-08-22 20:35:43 发布

阅读量968

点赞数 19

文章标签：算法 c++ 线性代数

本文链接：https://blog.csdn.net/HYC0907/article/details/141401463

版权

题目描述

试计算积分

$\int_{L}^{R}\frac{cx+d}{ax+b}dx$

结果保留至小数点后 6 位。

数据保证计算过程中分母不为 0 且积分能够收敛。

输入格式

一行，包含 6 个实数 a,b,c,d,L,R。

输出格式

一行，积分值，保留至小数点后 6 位。

输入输出样例

输入 #1

1 2 3 4 5 6

输出 #1

2.732937

说明/提示

a,b,c,d∈[−10,10]，−100≤L<R≤100 且 R−L≥1。

一、纯数学推导

①a 不等于 0

我们需要先判断 a 是否等于 0 ，如果不等于 0 ，那么：

1.变形

首先，我们可以把 $\frac{cx+d}{ax+b}$ 变成 $\frac{c}{a}+\frac{d-\frac{bc}{a}}{ax+b}$ 。为什么？

我们观察被积函数 $\frac{cx+d}{ax+b}$ ，为了更容易地积分，我们尝试将分子重写为与分母线性相关的形式加上一个常数项。具体来说，我们希望将分子 $cx+d$ 改写为 $k(ax+b)+m$ 的形式，其中 k 和 m 是常数。

通过比较系数，我们可以得到：

$cx+d=k(ax+b)+m$

展开并整理得：

$cx+d=kax+kb+m$

由于两边的 x 和常数项的系数必须相等，我们可以列出方程组：

$\left\{\begin{matrix} ka=c\\ kb+m=d\end{matrix}\right.$

解这个方程组，我们得到：

$k=\frac{c}{a},m=d-\frac{bc}{a}$

因此，分子可以重写为：

$cx+d=\frac{c}{a}(ax+b)+(d-\frac{bc}{a})$

即：

$\frac{cx+d}{ax+b}=\frac{c}{a}+\frac{d-\frac{bc}{a}}{ax+b}$

2.积分

我们再在上面的基础上进行积分：

$\int \frac{cx+d}{ax+b}dx=\int (\frac{c}{a}+\frac{d-\frac{bc}{a}}{ax+b}dx)$

$=\int \frac{c}{a}dx+\int \frac{d-\frac{bc}{a}}{ax+b}dx$

$=\frac{c}{a}x+\frac{d-\frac{bc}{a}}{a}\int \frac{1}{u}du$ （令 u = ax+b , du = a dx ）

$=\frac{c}{a}x+\frac{d-\frac{bc}{a}}{a}\ln \left | u \right |+C$

$=\frac{c}{a}x+\frac{d-\frac{bc}{a}}{a}\ln \left | ax+b \right |+C$ （代回 u = ax+b ）

3.转化

我们最后来转化一下，使这个式子变成一个更加简洁的式子。

$\frac{c}{a}x=\frac{cx}{a}$

$\frac{d-\frac{bc}{a}}{a}=\frac{ad-bc}{a^{2}}$

所以

$\frac{c}{a}x+\frac{d-\frac{bc}{a}}{a}\ln \left | ax+b \right |+C=\frac{cx}{a}+\frac{ad-bc}{a^{2}}\ln \left | ax+b \right |+C$

②a 等于 0

如果 a 等于 0 的话， $ax+b$ 就变成了 $0x+b=b$ ，积分就可以简化为：
$\int \frac{cx+d}{b}dx$

$=\frac{1}{b}\int (cx+d)dx$

$=\frac{1}{b}(\frac{c}{2}x^{2}+dx)+C$

$=\frac{c}{2b}x^{2}+\frac{d}{b}x+C$

③代码

所以，代码是：

#include<bits/stdc++.h>
using namespace std;
double a,b,c,d,l,r;
double f1(double x)
{
	return (c*x)/a+((a*d-b*c)/(a*a))*log(fabs(a*x+b));
}
double f2(double x)
{
	return (c/(2*b))*x*x+(d/b)*x;
}
int main() {
	cin>>a>>b>>c>>d>>l>>r;
	if(a==0.0)
	{
		printf("%.6f",f2(r)-f2(l));
	}
	else
	{
		printf("%.6f",f1(r)-f1(l));
	}
	return 0;
}

二、自适应辛普森法

当然，这道题可以用自适应辛普森法去做。

1.思路

自适应辛普森算法是一类近似算法，主要用于求较难反导的函数的积分。其思想是利用二次函数曲线来不断拟合所求曲线。而自适应则是用于优化时间复杂度的方法。

什么是拟合？

如果给定了一个函数，我们可以用一个图像与它近似的初等函数来代替它，这样的过程称之为“拟合”。

辛普森公式

$\int_{l}^{r}f(x)dx\approx \frac{r-l}{6}\begin{bmatrix} f(l)+f(r)+4f(\frac{l+r}{2})\end{bmatrix}$

辛普森公式推导

对于一个二次函数 $f(x)=ax^{2}+bx+c$ ，求积分可得：

$F(x)=\int_{0}^{x}f(x)dx=\frac{a}{3}x^{3}+\frac{b}{2}x^{2}+cx+D$

那么：

$\int_{l}^{r}f(x)dx=F(r)-F(l)$

$=\frac{a}{3}(r^{3}-l^{3})+\frac{b}{2}(r^{2}-l^{2})+c(r-l)$

$=(r-l)(\frac{a}{3}(l^{2}+r^{2}+lr)+\frac{b}{2}(l+r)+c)$

$=\frac{r-l}{6}(2al^{2}+2ar^{2}+2alr+3bl+3br+6c)$

$=\frac{r-l}{6}((al^{2}+bl+c)+(ar^{2}+br+c)+4(a(\frac{l+r}{2})^{2}+b(\frac{l+r}{2})+c))$

$=\frac{r-l}{6}(f(l)+f(r)+4f(\frac{l+r}{2}))$

我们可以用辛普森公式去求定积分。我们把积分区间分成很多的非常小的区间，每一个区间，我们都拟合一下。

关键是我们要分成多少段？

如果我们分得太少了，答案就会不精确；如果我们分得太多了，时间会太长。

所以，关键就是“自适应”。精度不够低，我们就再继续分；精度够了，我们就不再分下去了。核心过程就是，递归二分划分区间，精度达到标准，终止递归。否则继续递归。

模板：

const double eps=精度要求;
double f(double x)
{
	需要积分的函数
}
double simpson(double l,double r)
{
	return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6;
}
double answer(double l,double r,double ans)
{
	double mid=(l+r)*0.5;
	double fl=simpson(l,mid);
	double fr=simpson(mid,r);
	if(fabs(fl+fr-ans)<=eps)
	{
		return fl+fr;
	}
	else
	{
		return answer(l,mid,fl)+answer(mid,r,fr);
	}
}

2.代码

#include<bits/stdc++.h>
using namespace std;
double a,b,c,d,l,r;
double eps=1e-12;
double f(double x)
{
	return (c*x+d)/(a*x+b);
}
double simpson(double l,double r)
{
	return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6;
}
double answer(double l,double r,double ans)
{
	double mid=(l+r)*0.5;
	double fl=simpson(l,mid);
	double fr=simpson(mid,r);
	if(fabs(fl+fr-ans)<=eps)
	{
		return fl+fr;
	}
	else
	{
		return answer(l,mid,fl)+answer(mid,r,fr);
	}
}
int main()
{
	cin>>a>>b>>c>>d>>l>>r;
	printf("%.6f",answer(l,r,simpson(l,r)));
	return 0;
}