HDU - 1071 - The area - 高斯约旦消元法 - 自适应辛普森法积分
http://acm.hdu.edu.cn/showproblem.php?pid=1071
解一个给定三个点的坐标二次函数某区域的积分值。
设出方程之后高斯消元得到二次函数。然后再消元得到直线。
两次积分然后相减就可以了。
把自适应辛普森改成了传入函数指针的形式,有点多此一举。
可以这样做的原因,是因为这道题保证要求的区域都是在第一象限,否则不能直接定积分。
语言:G++
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
namespace Gauss_Jordan_Elimination {
const int MAXN=10;
double a[MAXN][MAXN],del;
double ans[MAXN];
const double eps=1e-8;
bool gauss_jordan(int n) {
for(int i=1; i<=n; i++) {
int mx=i;
for(int j=i+1; j<=n; j++)
if(fabs(a[j][i])>fabs(a[mx][i]))
mx=j;
for(int j=1; j<=n+1; j++)
swap(a[i][j],a[mx][j]);
if(fabs(a[i][i])<eps)
return false;
for(int j=1; j<=n+1; j++) {
if(j!=i) {
double tmp=a[j][i]/a[i][i];
for(int k=i+1; k<=n+1; k++) {
a[j][k]-=a[i][k]*tmp;
}
}
}
}
for(int i=1; i<=n; i++)
ans[i]=a[i][n+1]/a[i][i];
return true;
}
}
using namespace Gauss_Jordan_Elimination;
double A,B,C,D,E;
namespace Adaptive_Simpson_Integral {
/* 备注:
1.直接往上面模板的f函数中输入本题给定函数,然后调用一次asr(l,r)即可。
2.可能需要修改精度要求,测试过至少取需求精度的100倍或1000倍就挺好。
*/
double seps=1e-8;
double f(double x) {
//需要积分的函数
return A*x*x+B*x+C;
}
double g(double x) {
//需要积分的函数
return D*x+E;
}
double simpson(double l,double r,double (*F)(double)) {
double mid=(l+r)/2;
return (F(l)+4*F(mid)+F(r))*(r-l)/6;
}
double asr(double l,double r,double A,double (*F)(double)) {
double mid=(l+r)/2;
double L=simpson(l,mid,F),R=simpson(mid,r,F);
if(fabs(L+R-A)<=15*seps)
return L+R+(L+R-A)/15.0;
return asr(l,mid,L,F)+asr(mid,r,R,F);
}
double asr(double l,double r,double (*F)(double)) {
return asr(l,r,simpson(l,r,F),F);
}
/* 备注:
1.直接往上面模板的F函数中输入本题给定函数,然后调用一次asr(l,r)即可。
2.可能需要修改精度要求,测试过至少取需求精度的100倍或1000倍就挺好。
*/
}
using namespace Adaptive_Simpson_Integral;
double x[3],y[3];
int main() {
int t;
scanf("%d",&t);
while(t--) {
scanf("%lf%lf%lf%lf%lf%lf",&x[0],&y[0],&x[1],&y[1],&x[2],&y[2]);
for(int i=1; i<=3; i++) {
a[i][1]=x[i-1]*x[i-1];
a[i][2]=x[i-1];
a[i][3]=1;
a[i][4]=y[i-1];
}
gauss_jordan(3);
A=ans[1],B=ans[2],C=ans[3];
//cout<<"A="<<A<<" B="<<B<<" C="<<C<<endl;
for(int i=1; i<=2; i++) {
a[i][1]=x[i];
a[i][2]=1;
a[i][3]=y[i];
}
gauss_jordan(2);
D=ans[1],E=ans[2];
//cout<<"D="<<D<<" E="<<E<<endl;
double ans1=asr(x[1],x[2],f);
double ans2=asr(x[1],x[2],g);
printf("%.2f\n",ans1-ans2);
}
}