變步長梯形積分演算法求解函式定積分
阿新 • • 發佈:2019-01-01
演算法基本原理:把原區間分為一系列小期間(n份),在每個小區間上都用小的梯形面積來近似代替原函式的積分,當小區間足夠小時,就可以得到原來積分的近似值。但是這個過程中有一個問題是步長的取值,步長太大精度難以保證,步長太小會導致計算量的增加,所以,實際計算中常常採用變步長的方法,使得在步長逐次減小的過程中,求得的積分結果滿足要求的精度為止。
首先,給出兩個計算公式
(1)
(2) //將步長h二分,計算以h/2為步長的積分
步驟:
- 取n=1,利用公式(1)計算積分值
- 進行二分,利用新的公式(2)計算新的積分值
- 進行判斷,如果兩次計算積分值的差在給定的誤差範圍之內,二分後的積分值就是所需的結果,計算結束,否則返回第二步繼續執行
因此,主要計算的問題有兩個
一. 被積函式值的運算(設定Function類)
二. 變步長梯形積分的實現(設定Trapz類)
#ifndef _TRAPZINT_H_ #define _TRAPZINT_H_ class Function{ public: virtual double operator()(double x) const = 0;//純虛擬函式過載運算子() virtual ~Function(){}//解構函式 }; class MyFunction :public Function{ public: virtual double operator()(double x) const;//覆蓋虛擬函式 }; class Integration{//抽象類Integration的定義 public: virtual double operator()(double a, double b, double eps) const = 0; virtual ~Integration(){} }; class Trapz :public Integration{//公有派生類Trapz的定義 public: Trapz(const Function &f) :f(f){}//建構函式 virtual double operator()(double a, double b, double eps) const;//覆蓋虛擬函式 private: const Function &f;//私有成員,Function類物件的指標 }; #endif
#include"Trapzint.h"
#include<iostream>
#include<cmath>
using namespace std;
//計算被積函式的函式值
double MyFunction::operator()(double x) const{
return log(1.0 + x) / (1.0 + x*x);
}
//積分運算過程,過載為運算子()
double Trapz::operator()(double a, double b, double eps) const{
//計算區間[a,b]之間誤差在eps之內的積分值
bool done = false;//Trapz類的虛擬函式成員
int n = 1;
double h = b - a;
double tn = h*(f(a) + f(b)) / 2;//計算n=1時的積分
double t2n;//步長為一半的時候的積分值
do{
double sum = 0;
for (int k = 0; k < n; k++){
double x = a + (k + 0.5)*h;
sum += f(x);
}
t2n = (tn + h*sum) / 2.0;//利用公式(2)計算積分
if (fabs(tn - t2n) < eps)
done = true;//判斷積分誤差
else{
tn = t2n;
n *= 2;
h /= 2;
}
} while (!done);
return t2n;
}
#include"Trapzint.h"
#include<iostream>
#include<iomanip>
using namespace std;
int main(){
MyFunction f;//定義MyFunction類的物件
Trapz trapz(f);//定義Trapz類的物件
//計算並輸出結果
cout << "Trapz Int:" << setprecision(7) << trapz(0, 2, 1e-7) << endl;
return 0;
}
結果