LA3485二分+求解积分方程+辛普森算法计算积分

 1 /*LA3485:
 2 求解积分方程
 3 关于这道题的数学模型：
 4 给定抛物线长度L，抛物线函数f(x)=a(x-d)(x+d),求解 |a*d*d|的值，a>0,曲线积分函数lf(x)=sqrt(1+f'(x)^2)
 5 a越大，|a*d*d|越大，L越长，所以可以二分求解
 6 
 7 辛普森算法解函数f(x)在区间(a,b)上的积分：
 8 模板如下：
 9 double simpsonF(double a,double b);//返回测试值
10 double simpsonM(double a,double b,eps,double A);//自适应simpson递归
11 //答案是：simpsonM(a,b,eps,simpsonF(a,b));
12 double simpsonF(double a,double b)
13 {
14 double c=a+(b-a)/2;
15 return (F(a)+4*F(c)+F(b))*(b-a)/6;
16 }
17 double simpsonM(double a,double b,double eps,double A)
18 {
19 double c=a+(b-a)/2;
20 double L=simpsonF(a,c),R=simpsonF(c,b);
21 if (fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15.0;
22 return simpsonM(a,c,eps/2,L)+simpsonM(c,b,eps/2,R);
23 }
24 简要原理分析：细节部分自己以后需要学习！
25 传统求微积分的方法，是手动计算，但是很多难以写出被积函数，这里用计算机帮助计算。
26 任然使用先微分，再求和的方法。微分就是等分区间，去重点函数值，近似成矩形，最终将面积求和。
27 如果直接使用这个方法，越细分，近似程度越高，计算量越大，不细分，精度不够。以前遇到题目，用朴素算法写错了。
28 我们在函数曲线平缓的地方，少细分一些，在陡峭的地方（容易损失精度）多细分。折中的办法满足了时间和精度的要求。
29 
30 */
31 
32 
33 #include<iostream>
34 #include<stdio.h>
35 #include<string.h>
36 #include<algorithm>
37 #include<stdlib.h>
38 #include<math.h>
39 #include<queue>
40 #include<vector>
41 #include<map>
42 #define LL long long
43 using namespace std;
44 double d,l;
45 double F(double k,double x)//k是枚举的系数
46 {
47     k=k/d/d;
48     return sqrt(1+(2*k*(x))*(2*k*(x)));
49 }
50 double simpsonF(double a,double b,double k)
51 {
52 double c=a+(b-a)/2;
53 return (F(a,k)+4*F(c,k)+F(b,k))*(b-a)/6;
54 }
55 double simpsonM(double a,double b,double eps,double A,double k)
56 {
57 double c=a+(b-a)/2;
58 double L=simpsonF(a,c,k),R=simpsonF(c,b,k);
59 if (fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15.0;
60 return simpsonM(a,c,eps/2,L,k)+simpsonM(c,b,eps/2,R,k);
61 }
62 double Q(double k)
63 {
64     double ans=simpsonM(-d,d,1e-7,simpsonF(-d,d,k),k);
65     return ans-l;
66 }
67 int T,D,H,B,L;
68 int main()
69 {
70 //    freopen("out.txt","w",stdout);
71     cin>>T;
72     for(int cas=1;cas<=T;cas++)
73     {
74         cin>>D>>H>>B>>L;
75         d=(B+0.0)/ceil((B+0.0)/D)/2.0;
76         l=(L+0.0)/ceil((B+0.0)/D);
77         double l1=0,r1=H;
78         while(r1-l1>1e-5)
79         {
80             double m=(l1+r1)/2;
81             if (Q(m)>0) r1=m;else l1=m;
82         }
83         if(cas>1) cout<<endl;
84         printf("Case %d:
",cas);
85         printf("%.2lf
",H-l1);
86 
87     }
88     return 0;
89 }