|
//數(shù)值計算實驗 數(shù)值積分
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <gsl/gsl_integration.h>
using namespace std;
//被積函數(shù)
double f(double x)
{
//為便于調(diào)試，先弄個有原函數(shù)的 y = x^2 + x^3 - 2*x^4
return x*x + x*x*x - 2*x*x*x*x;
}
//被積函數(shù) 給gsl用的
double g(double x, void * params)
{
return f(x);
}
//原函數(shù) 用于調(diào)試算法
double F(double x)
{
return (x*x*x)/3.0 + (x*x*x*x)/4.0 - 0.4*x*x*x*x*x;
}
//輸出被積函數(shù)的精確解
double Jinque(const double a, const double b)
{
return (F(b) - F(a));
}
//梯形法 求函數(shù)在[a,b]上的定積分，積分區(qū)間分為n部分
double Tixing(const double & a, const double & b, const int & n)
{
double sum = 0.0;
double gaps = (b-a)/double(n); //每個間隔的長度
for (int i = 0; i < n; i++)
{
sum += (gaps/2.0) * (f(a + i*gaps) + f(a + (i+1)*gaps));
}
return sum;
}
//拋物線法
double Paowuxian(const double & a, const double & b, const int & n)
{
double sum = 0.0;
double gaps = (b-a)/double(n); //每個間隔的長度
double h = gaps/2.0;
for (int i = 0; i < n; i++)
{
sum += (h/3.0) * (f(a + i*gaps) + f(a + (i+1)*gaps) + 4.0*f((2*a + (2*i+1)*gaps)/2.0));
}
return sum;
}
//柯特斯公式
double Cotes(const double & a, const double & b, const int & n)
{
double sum = 0.0;
double gaps = (b-a)/double(n); //每個間隔的長度
double h = gaps/2.0;
for (int i = 0; i < n; i++)
{
sum += (h/45.0) * (7.0*f(a + i*gaps) +
32.0*f(a + i*gaps + 0.25*gaps) +
12.0*f(a + i*gaps + 0.5*gaps) +
32.0*f(a + i*gaps + 0.75*gaps) +
7.0*f(a + (i+1)*gaps));
}
return sum;
}
//gsl解法，參考gsl文檔
double gslIntegration(double & a, double & b)
{
gsl_function gf;
gf.function = g;
double r, er;
unsigned int n;
gsl_integration_qng(&gf, a, b, 1e-10, 1e-10, &r, &er, &n);
return r;
}
int main()
{
double a, b;
int n;
cout<<"請輸入積分區(qū)間:"<<endl;
cout<<"a = ";
cin>>a;
cout<<"b = ";
cin>>b;
cout<<"請輸入分割被積區(qū)間的數(shù)量:";
cin>>n;
if (a > b || n <= 1)
{
cout<<"輸入錯誤！"<<endl;
exit(1);
}
//設置輸出精度
cout.precision(10);
//輸出精確解
double result = Jinque(a, b);
cout<<"函數(shù)在["<<a<<","<<b<<"]上的定積分為:"<<result<<endl;
//梯形法
double result1 = Tixing(a, b, n);
cout<<"梯形法:"<<endl;
cout<<"函數(shù)在["<<a<<","<<b<<"]上的定積分為:"<<result1<<" 相對誤差為:"
<<abs((result1 - result)/result)*100<<"%"<<endl;
//拋物線法
double result2 = Paowuxian(a, b, n);
cout<<"拋物線法:"<<endl;
cout<<"函數(shù)在["<<a<<","<<b<<"]上的定積分為:"<<result2<<" 相對誤差為:"
<<abs((result2 - result)/result)*100<<"%"<<endl;
//柯特斯公式法
double result3 = Cotes(a, b, n);
cout<<"柯特斯法:"<<endl;
cout<<"函數(shù)在["<<a<<","<<b<<"]上的定積分為:"<<result3<<" 相對誤差為:"
<<abs((result3 - result)/result)*100<<"%"<<endl;
//調(diào)用gsl函數(shù)
double result4 = gslIntegration(a, b);
cout<<"gsl函數(shù)結(jié)果:"<<endl;
cout<<"函數(shù)在["<<a<<","<<b<<"]上的定積分為:"<<result4<<" 相對誤差為:"
<<abs((result4 - result)/result)*100<<"%"<<endl;
return 0;
}
|
