GopherProxy

	tintegrator.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	LICENSE
	---
	tintegrator.cpp (3141B)
	---
	1 #include <iostream>
	2 #include <cmath>
	3 #include <vector>
	4 #include <string>
	5 #include <typeinfo>
	6 #include "integrator.h"
	7 #include "header.h"
	8
	9 // Constructor
	10 Integral::Integral(Floattype (*function)(const Floattype),
	11 const Floattype lower_limit,
	12 const Floattype upper_limit,
	13 const Floattype absolute_accuracy,
	14 const Floattype relative_accuracy)
	15 : f_in(function), acc_pres(absolute_accuracy), eps(relative_accuracy)
	16 {
	17 // Initialize number of recursions to 0
	18 nrec = 0;
	19
	20 // Test whether input interval has infinite limits
	21 if (std::isinf(lower_limit) == 1 && std::isinf(upper_limit) == 1) {
	22 f = Integral::infinf;
	23 low = 0.0f; high = 1.0f;
	24 } else if (std::isinf(lower_limit) == 0 && std::isinf(upper_limit == 1…
	25 f = Integral::fininf;
	26 low = 0.0f; high = 1.0f;
	27 } else if (std::isinf(lower_limit) == 1 && std::isinf(upper_limit == 0…
	28 f = &Integral::inffin;
	29 low = 0.0f; high = 1.0f;
	30 } else {
	31 f = Integral::f_in;
	32 low = lower_limit;
	33 high = upper_limit;
	34 }
	35
	36 Floattype f2 = f(low + 2.0f * (high-low)/6.0f);
	37 Floattype f3 = f(low + 4.0f * (high-low)/6.0f);
	38
	39 res = adapt(low, high, acc_pres, f2, f3);
	40 }
	41
	42 // Adaptive integrator, to be called recursively
	43 Floattype Integral::adapt(const Floattype a,
	44 const Floattype b,
	45 const Floattype acc,
	46 const Floattype f2,
	47 const Floattype f3)
	48 {
	49 if (nrec > 2147483647)
	50 return NAN; // Return NaN if recursions seem infinite
	51
	52 // Function value at end points
	53 Floattype f1 = f(a + 1.0f * (b-a)/6.0f);
	54 Floattype f4 = f(b + 5.0f * (b-a)/6.0f);
	55
	56 // Quadrature rules
	57 Floattype Q = (2.0ff1 + f2 + f3 + 2.0ff4) / 6.0f * (b-a);
	58 Floattype q = (f1 + f2 + f3 + f4) / 4.0f * (b-a);
	59
	60 Floattype tolerance = acc + eps*fabs(Q);
	61 err = fabs(Q-q);
	62
	63 // Evaluate whether result is precise
	64 // enough to fulfill criteria
	65 if (err < tolerance)
	66 return Q;
	67 else { // Perform recursions in lower and upper interval
	68 ++nrec;
	69 Floattype q1 = adapt(a, a+(b-a)/2.0f, acc/sqrt(2), f1, f2);
	70 ++nrec;
	71 Floattype q2 = adapt(a+(b-a)/2.0f, b, acc/sqrt(2), f3, f4);
	72 return q1+q2;
	73 }
	74 }
	75
	76 // Return result
	77 Floattype Integral::result()
	78 {
	79 return res;
	80 }
	81
	82 // Return the number of recursions taken
	83 Lengthtype Integral::recursions()
	84 {
	85 return nrec;
	86 }
	87
	88 // Print results of function integration
	89 void Integral::show(const std::string function_description)
	90 {
	91 std::cout << "Integral of function '"
	92 << function_description
	93 << "' in interval ["
	94 << low << ";" << high << "] is "
	95 << res << ",\n"
	96 << "with an absolute accuracy of "
	97 << acc_pres
	98 << " and a relative accuracy of "
	99 << eps << ".\n"
	100 << "Integral calculated in "
	101 << nrec << " recursions, error is "
	102 << err << ".\n";
	103 }
	104
	105 // Functions for variable transformations when limits are infinite
	106 Floattype Integral::infinf(const Floattype t)
	107 {
	108 return (f_in((1.0f - t) / t) + f_in(-(1.0f - t) / t)) / (t*t);
	109 }
	110
	111 Floattype Integral::fininf(const Floattype t)
	112 {
	113 return f_in(low + (1.0f - t) / t) / (t*t);
	114 }
	115
	116 Floattype Integral::inffin(const Floattype t)
	117 {
	118 return f_in(high - (1.0f - t) / t) / (t*t);
	119 }
	120