#include <iostream> #include <iomanip> #include <cmath> using namespace std; Quadrature: trapezoidal formula double func_integral_trapezoid (double lo, double Hi, double (*func) (double), int n = +) {//lo: Lower limit; hi: Upper limit; Func: letter Number N: Equal fraction if (n <= 0) n = 100; Double X; Double step = (Hi-lo)/n; Double result = 0.0; x = lo; for (int i = 1; i < n; i++) {x + = step; result + = Func (x);} result + = (Func (LO) + Func (HI))/2; Result *= step; return result; }//quadrature: Simpson formula double Func_integral_simpson (double lo, double Hi, double (*func) (double), int n = +) {//lo: Lower limit; hi: Upper limit; F UNC: function N: Equal fraction if (n <= 0) n = 100; Double X; Double step = (Hi-lo)/n; Double result1 = 0.0; x = lo; for (int i = 1; i < n; i++) {x + = step; Result1 + = Func (x);} result1 *= 2; Double result2 = 0.0; x = lo + step/2; for (int i = 0; i < n; i++) {result2 + = Func (x); x + = step;} result2 *= 4; Double result = result1 + result2 + func (LO) + func (HI); Result *= STEP/6; return result; }//quadrature: Cotes formula double Func_integRal_cotes (double lo, double Hi, double (*func) (double), int n = +) {//lo: Lower bound; hi: Upper limit; Func: function; N: Equal fraction if (n <= 0) n = 100; Double X; Double step = (Hi-lo)/n; Double result1 = 0.0; x = lo; for (int i = 1; i < n; i++) {x + = step; Result1 + = Func (x);} result1 *= 14; Double result2 = 0.0; x = lo + step/2; Double result3 = 0.0; Double x1 = lo + step/4; Double x2 = lo + STEP/4 * 3; for (int i = 0; i < n; i++) {result2 + = func (x), RESULT3 + = func (x1) + func (x2), x + = step, x1 + = step, x2 + = step; RESULT2 *= 12; RESULT3 *= 32; Double RESULT4 = (func (LO) + Func (HI)) * 7; Double result = result1 + result2 + result3 + result4; Result *= step/90; return result; }//quadrature: Romberg formula double Func_integral_romberg (double lo, double Hi, double (*func) (double), int k = 4) {//lo: Lower bound; hi: Upper limit; Func : function; k: equal exponent (interval divided into 2^k parts) int size = k + 1; Double *matrix = new Double[size*size]; for (int i = 0; i < size*size; i++) matrix[i] = 0.0; Double step = Hi-lo; Matrix[0] = Func_integral_trapezoid (lo, Hi, Func, 1); for (int i = 1; i < size; i++) {int n = 1 << (i-1), for (int k = 0; k < n, k++) {matrix[i*size + 0] + = Fun C (lo + (k + 0.5) *step); } matrix[i*size + 0] *= step; Matrix[i*size + 0] + = matrix[(i-1) *size]; Matrix[i*size + 0]/= 2.0; Step/= 2.0; } Double temp = 1.0; Double Factor1, Factor2; for (int j = 1; J < Size, J + +) {temp *= 4.0; Factor1 = temp/(temp-1); Factor2 = 1/(temp-1); for (int i = j; i < size; i++) {matrix[i*size + j] = Factor1*matrix[i*size + j-1]-factor2*matrix[(i-1) *size + j-1];}} Double result = Matrix[k*size + K]; Delete[] Matrix; return result; }//Test integrand, 0 to 1 points for Pi double func_test1 (double x) {return 4/(1 + x*x);} int main () {cout << trapezoid numerica L integration "<< Endl; cout << setprecision << func_integral_trapezoid (0, 1, func_test1, 1024x768) << Endl << Endl; cout << "Simpson numerical integration" << Endl; cout << setprecision () << Func_integral_simpson (0, 1, func_test1, 1024x768) << Endl << Endl; cout << "Cotes numerical integration" << Endl; cout << setprecision << func_integral_cotes (0, 1, func_test1, 1024x768) << Endl << Endl; cout << "Romberg numerical integraion" << Endl; cout << setprecision () << Func_integral_romberg (0, 1, Func_test1, ten) << Endl << Endl; GetChar (); return 0; }
C + + implementation of a single integral