跟着PoPoQQQ学算法

最新推荐文章于 2023-04-12 11:19:05 发布

mrcrack

最新推荐文章于 2023-04-12 11:19:05 发布

阅读量512

点赞数 1

分类专栏：跟着大佬学算法

本文链接：https://blog.csdn.net/mrcrack/article/details/88068382

版权

跟着大佬学算法专栏收录该内容

154 篇文章 8 订阅

订阅专栏

跟着PoPoQQQ学算法

PoPoQQQ博客https://blog.csdn.net/popoqqq

此文不包括与https://blog.csdn.net/mrcrack/article/details/84471041雷同的算法

测评地址:

Simpson 积分辛普森积分

1.【模板】自适应辛普森法1 //在线测评地址https://www.luogu.org/problemnew/show/P4525

2.【模板】自适应辛普森法2 //在线测评地址https://www.luogu.org/problemnew/show/P4526

题解:

1.【模板】自适应辛普森法1

//P4525 【模板】自适应辛普森法1
//在线测评地址https://www.luogu.org/problemnew/show/P4525
//该题不定积分展开形式c*x/a+(d-b*c/a)/a*ln(ax+b)+常数
//通过此文https://blog.csdn.net/u014399084/article/details/70249206大致了解C++头文件位置
//想看看,头文件中的ln函数.没找到,求助网络https://bbs.csdn.net/topics/390620134?page=1发现ln就是log
//测试了log(2)=0.693147确实如此,用linux自带计算器检验了ln2=0.693147181
//样例通过,提交40分,测试点5-10WA,根据不定积分算出的结果也会出问题?
//以下为40分代码.2019-3-1 22:21
//翻看题解https://www.luogu.org/problemnew/solution/P4525比较偏好作者: iki9 更新时间: 2018-12-23 20:08写法
//发现c*x/a+(d-b*c/a)/a*ln(ax+b)+常数应写成c*x/a+(d-b*c/a)/a*ln|ax+b|+常数漏了绝对值
//还需讨论a=0的情况.不定积分展开形式c*x*x/2/b+d*x/b+常数
//以下为采用不定积分计算的AC代码.2019-3-1 22:43
//提一点,C中fabs用在浮点,abs用在int;而C++中abs可以用在浮点,提一点,这个地方若有疑问,还是要及时测试.
//有个疑问,c*x/a+(d-b*c/a)/a*ln|ax+b|+常数绝对值怎么来的,当然ln里的数要大于0,有机会要看看微积分解决这个问题
//辛普森公式原理可看此文https://www.cnblogs.com/widsom/p/7694882.html配图说明
//辛普森公式推导,可以下3篇文章结合起来看
//https://blog.csdn.net/VictoryCzt/article/details/80660113
//https://blog.csdn.net/fengzhizi76506/article/details/54848920
//https://blog.csdn.net/xyz32768/article/details/81392369
//a^3-b^3=(a-b)(a^2+ab+b^2)可这样得到
//(a-b)(Aa^2+Bab+Cb^2)=Aa^3-Cb^3+Ba^2b+Cab^2-Aa^2b-Bab^2=Aa^3-Cb^3+(B-A)a^2b+(C-B)ab^2
//可得A=B,B=C,A=1,C=1,B=1
//喜欢此文代码的写法https://www.cnblogs.com/widsom/p/7694882.html
//测试样例时,又遇到段错误这次不慌了
//排查在递归,发现double mid=(left+right)/2;//此处写成double mid=(left+right);造成递归溢出,即一直递归
//修改,还是遇到段错误
//排查,return(right-left)/6*( f(right)+4*f(mid)+f(left));//此处写成return f(right)+4*f(mid)+f(left);
//修改,样例通过,提交AC.2019-3-2 13:39
//新的东西,再简单,还是容易手忙脚乱,熟练确实需要经历一定的操练.
//以下为Simpson 积分辛普森积分 AC代码.2019-3-2- 13:41
#include <stdio.h>
#include <math.h>
#define eps 1e-10
double a,b,c,d;
double f(double x){
   return (c*x+d)/(a*x+b);
}
double Simpson(double left,double right){
   double mid=(left+right)/2;
   return(right-left)/6*( f(right)+4*f(mid)+f(left));//此处写成return f(right)+4*f(mid)+f(left);
}
double Integral(double left,double right){//积分
   double mid=(left+right)/2;//此处写成double mid=(left+right);造成递归溢出,即一直递归
   double left_ans,right_ans,ans;
   left_ans=Simpson(left,mid),right_ans=Simpson(mid,right),ans=Simpson(left,right);
   if(fabs(left_ans+right_ans-ans)<=eps)return ans;
   return Integral(left,mid)+Integral(mid,right);
}
int main(){
   double left,right;
   scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&left,&right);
   printf("%lf\n",Integral(left,right));
   return 0;
}

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

//P4526 【模板】自适应辛普森法2
//在线测评地址https://www.luogu.org/problemnew/show/P4526
//用不了积分公式,那么 Simpson 积分辛普森积分
//困扰就是,积分发散,如何判定
//先编能有输出的情况
//还有,上界无穷,如何处理
//样例输出一直为0,
//仔细一想,f(0)=0,f(INF)=0,有问题
//应分成两段,0-sqrt(a),sqrt(a)-INF
//a/x-x => x=sqrt(a)应是极大值
//printf("%.5lf\n",Integral(eps,sqrt(a))+Integral(sqrt(a),INF));//此处写成printf("%.5lf\n",Integral(eps,INF));
//进行上述修改,样例通过,
//因,sqrt(a),需分成,a=0,a<0,a>0,3种情况进行讨论
//测试了,发现a=0,a<0 段错误,果断将这两种情况归为积分发散
//提交90分,测试点6WA.总体满意.2019-3-2 21:57
//翻看讨论版https://www.luogu.org/discuss/lists?forumname=P4526
//大致有数,测试点6,数据a=0的情况,从该题角度,有解
//else if(a==0)printf("%.5lf\n",Integral(eps,INF));//添加此种情况
//修改,提交AC.2019-3-2 22:06
//在处理该题中,还是对x^(a/x-x)进行了讨论,看来,数学的函数知识还是蛮要紧的.2019-3-2 22:12
#include <stdio.h>
#include <math.h>
#define eps 1e-10
#define INF 1e10
double a;
double f(double x){
   return pow(x,a/x-x);
}
double Simpson(double left,double right){
   return (right-left)/6*(f(left)+4*f((left+right)/2)+f(right));//此处写成return (right-left)/6*(f(left)+f((left+right)/2)+f(right));
}
double Integral(double left,double right){
   double mid=(left+right)/2;
   double left_ans,right_ans,ans;
   left_ans=Simpson(left,mid),right_ans=Simpson(mid,right),ans=Simpson(left,right);
   if(fabs(left_ans+right_ans-ans)<=eps)return ans;//此处写成if(fabs(left_ans+right_ans-ans)<=ans)return ans;昏招
   return Integral(left,mid)+Integral(mid,right);
}
int main(){
   scanf("%lf",&a);
   if(a>0)printf("%.5lf\n",Integral(0,sqrt(a))+Integral(sqrt(a),INF));//此处写成printf("%.5lf\n",Integral(eps,INF));
   else if(a==0)printf("%.5lf\n",Integral(eps,INF));//添加此种情况
   else printf("orz\n");
   return 0;
}