1421 lines
47 KiB
C++
1421 lines
47 KiB
C++
/****************************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
Integration.cc
|
|
|
|
Defines all functions to perform integration of functions.
|
|
|
|
******************************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
|
|
namespace ABACUS {
|
|
|
|
// ********************* class I_table ******************************
|
|
|
|
I_table::I_table (DP (*function_ref) (DP, DP), DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec) {
|
|
|
|
function = function_ref;
|
|
|
|
if (!(*this).Load(rhomin_ref, rhomax_ref, Nvals_ref, req_prec)) {
|
|
|
|
Nvals = Nvals_ref;
|
|
rhomin = rhomin_ref;
|
|
rhomax = rhomax_ref;
|
|
prec = req_prec;
|
|
rho_tbl = new DP[Nvals + 1];
|
|
I_tbl = new DP[Nvals + 1];
|
|
|
|
logalpha = log(rhomax/rhomin)/Nvals;
|
|
alpha = exp(logalpha);
|
|
|
|
for (int i = 0; i <= Nvals ; ++i) {
|
|
rho_tbl[i] = rhomin * pow(alpha, DP(i));
|
|
I_tbl[i] = function (rho_tbl[i], req_prec);
|
|
}
|
|
|
|
(*this).Save();
|
|
}
|
|
}
|
|
|
|
|
|
DP I_table::Return_val (DP req_rho) {
|
|
|
|
DP used_rho = fabs(req_rho);
|
|
|
|
if (used_rho < rhomin || used_rho >= rhomax)
|
|
{
|
|
cerr << "requesting I of " << setprecision(16) << used_rho << " with rho_min = " << rhomin << endl;
|
|
return (function (used_rho, prec));
|
|
}
|
|
|
|
else {
|
|
int index = int(log(used_rho/rhomin)/logalpha); // index between 0 and Nvals - 1
|
|
|
|
// Do linear interpolation between values at index and index + 1
|
|
|
|
return((rho_tbl[index + 1] - used_rho) * I_tbl[index] + (used_rho - rho_tbl[index]) * I_tbl[index + 1])
|
|
/(rho_tbl[index + 1] - rho_tbl[index]);
|
|
|
|
}
|
|
}
|
|
|
|
void I_table::Save () {
|
|
|
|
stringstream outfile_strstream;
|
|
outfile_strstream << "I_rhomin_" << rhomin << "_rhomax_" << rhomax
|
|
<< "_Nvals_" << Nvals << "_prec_" << prec << ".dat";
|
|
string outfile_str = outfile_strstream.str();
|
|
const char* outfilename = outfile_str.c_str();
|
|
|
|
ofstream outfile;
|
|
outfile.open(outfilename);
|
|
outfile.precision(16);
|
|
|
|
outfile << rho_tbl[0] << "\t" << I_tbl[0];
|
|
for (int i = 1; i <= Nvals; ++i) outfile << endl << rho_tbl[i] << "\t" << I_tbl[i];
|
|
|
|
outfile.close();
|
|
}
|
|
|
|
bool I_table::Load (DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec) {
|
|
|
|
stringstream infile_strstream;
|
|
infile_strstream << "I_rhomin_" << rhomin_ref << "_rhomax_" << rhomax_ref
|
|
<< "_Nvals_" << Nvals_ref << "_prec_" << req_prec << ".dat";
|
|
string infile_str = infile_strstream.str();
|
|
const char* infilename = infile_str.c_str();
|
|
|
|
ifstream infile;
|
|
infile.open(infilename);
|
|
|
|
if (infile.fail()) return (false);
|
|
|
|
Nvals = Nvals_ref;
|
|
rhomin = rhomin_ref;
|
|
rhomax = rhomax_ref;
|
|
prec = req_prec;
|
|
rho_tbl = new DP[Nvals + 1];
|
|
I_tbl = new DP[Nvals + 1];
|
|
|
|
logalpha = log(rhomax/rhomin)/Nvals;
|
|
alpha = exp(logalpha);
|
|
|
|
int i = 0;
|
|
while (infile.peek() != EOF) {
|
|
infile >> rho_tbl[i] >> I_tbl[i];
|
|
i++;
|
|
}
|
|
|
|
if (i != Nvals + 1) ABACUSerror("Bad input file in I.Load.");
|
|
|
|
return(true);
|
|
}
|
|
|
|
|
|
// ********************* class Integral_table ******************************
|
|
|
|
Integral_table::Integral_table (DP (*function_ref) (DP, DP, int), const char* filenameprefix,
|
|
DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec, int max_nr_pts)
|
|
{
|
|
function = function_ref;
|
|
|
|
if (!(*this).Load(filenameprefix, rhomin_ref, rhomax_ref, Nvals_ref, req_prec, max_nr_pts)) {
|
|
|
|
Nvals = Nvals_ref;
|
|
rhomin = rhomin_ref;
|
|
rhomax = rhomax_ref;
|
|
prec = req_prec;
|
|
maxnrpts = max_nr_pts;
|
|
|
|
rho_tbl = new DP[Nvals + 1];
|
|
I_tbl = new DP[Nvals + 1];
|
|
|
|
logalpha = log(rhomax/rhomin)/Nvals;
|
|
alpha = exp(logalpha);
|
|
|
|
for (int i = 0; i <= Nvals ; ++i) {
|
|
rho_tbl[i] = rhomin * pow(alpha, DP(i));
|
|
I_tbl[i] = function (rho_tbl[i], req_prec, max_nr_pts);
|
|
}
|
|
|
|
(*this).Save(filenameprefix);
|
|
}
|
|
}
|
|
|
|
DP Integral_table::Return_val (DP req_rho) {
|
|
|
|
if (req_rho < rhomin || req_rho >= rhomax)
|
|
{
|
|
cerr << "Requesting I of " << setprecision(16) << req_rho << " with rho_min = " << rhomin << endl;
|
|
return (function (req_rho, prec, maxnrpts));
|
|
}
|
|
|
|
else {
|
|
int index = int(log(req_rho/rhomin)/logalpha); // index between 0 and Nvals - 1
|
|
|
|
// Do linear interpolation between values at index and index + 1
|
|
|
|
return((rho_tbl[index + 1] - req_rho) * I_tbl[index] + (req_rho - rho_tbl[index]) * I_tbl[index + 1])/(rho_tbl[index + 1] - rho_tbl[index]);
|
|
}
|
|
}
|
|
|
|
void Integral_table::Save (const char* filenameprefix) {
|
|
|
|
stringstream outfile_strstream;
|
|
outfile_strstream << "Integral_table_" << filenameprefix << "_rhomin_" << rhomin << "_rhomax_" << rhomax
|
|
<< "_Nvals_" << Nvals << "_prec_" << prec << "_maxnrpts_" << maxnrpts << ".dat";
|
|
string outfile_str = outfile_strstream.str();
|
|
const char* outfilename = outfile_str.c_str();
|
|
|
|
ofstream outfile;
|
|
outfile.open(outfilename);
|
|
outfile.precision(16);
|
|
|
|
outfile << rho_tbl[0] << "\t" << I_tbl[0];
|
|
for (int i = 1; i <= Nvals; ++i) outfile << endl << rho_tbl[i] << "\t" << I_tbl[i];
|
|
|
|
outfile.close();
|
|
}
|
|
|
|
bool Integral_table::Load (const char* filenameprefix, DP rhomin_ref, DP rhomax_ref,
|
|
int Nvals_ref, DP req_prec, int max_nr_pts) {
|
|
|
|
stringstream infile_strstream;
|
|
infile_strstream << "Integral_table_" << filenameprefix << "_rhomin_" << rhomin_ref << "_rhomax_" << rhomax_ref
|
|
<< "_Nvals_" << Nvals_ref << "_prec_" << req_prec << "_maxnrpts_" << max_nr_pts << ".dat";
|
|
string infile_str = infile_strstream.str();
|
|
const char* infilename = infile_str.c_str();
|
|
|
|
ifstream infile;
|
|
infile.open(infilename);
|
|
|
|
if (infile.fail()) {
|
|
return (false);
|
|
}
|
|
|
|
Nvals = Nvals_ref;
|
|
rhomin = rhomin_ref;
|
|
rhomax = rhomax_ref;
|
|
prec = req_prec;
|
|
maxnrpts = max_nr_pts;
|
|
|
|
rho_tbl = new DP[Nvals + 1];
|
|
I_tbl = new DP[Nvals + 1];
|
|
|
|
logalpha = log(rhomax/rhomin)/Nvals;
|
|
alpha = exp(logalpha);
|
|
|
|
int i = 0;
|
|
while (infile.peek() != EOF) {
|
|
infile >> rho_tbl[i] >> I_tbl[i];
|
|
i++;
|
|
}
|
|
|
|
if (i != Nvals + 1) ABACUSerror("Bad input file in Integral_table.Load.");
|
|
|
|
return(true);
|
|
}
|
|
|
|
|
|
|
|
// *********************************** Simple Riemann sum *********************************************
|
|
|
|
DP Integrate_Riemann (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax, int Npts)
|
|
{
|
|
if (xmax < xmin) return(-Integrate_Riemann (function, args, arg_to_integ, xmax, xmin, Npts));
|
|
|
|
DP dx = (xmax - xmin)/Npts;
|
|
|
|
DP result = 0.0;
|
|
|
|
for (int i = 0; i < Npts; ++i) {
|
|
args[arg_to_integ] = xmin + (i + 0.5) * dx;
|
|
result += function (args);
|
|
}
|
|
result *= dx;
|
|
|
|
return(result);
|
|
}
|
|
|
|
DP Integrate_Riemann_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
|
|
DP xmin, DP xmax, int Npts)
|
|
{
|
|
if (xmax < xmin) return(-Integrate_Riemann_using_table (function, args, arg_to_integ, Itable, xmax, xmin, Npts));
|
|
|
|
DP dx = (xmax - xmin)/Npts;
|
|
|
|
DP result = 0.0;
|
|
|
|
for (int i = 0; i < Npts; ++i) {
|
|
args[arg_to_integ] = xmin + (i + 0.5) * dx;
|
|
result += function (args, Itable);
|
|
}
|
|
result *= dx;
|
|
|
|
return(result);
|
|
}
|
|
|
|
|
|
|
|
// *********************************** Adaptive Riemann sum (straight recursive) *******************************************
|
|
|
|
DP Integrate_rec_main (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax, DP req_prec, int rec_level, int max_rec_level, DP* f)
|
|
{
|
|
if (rec_level > max_rec_level) {
|
|
return((xmax - xmin) * (f[0] + f[1] + f[2])/3.0);
|
|
}
|
|
|
|
DP* f1 = new DP[9];
|
|
|
|
DP dx = (xmax - xmin)/9.0;
|
|
|
|
args[arg_to_integ] = xmin + 0.5 * dx;
|
|
f1[0] = function (args);
|
|
f1[1] = f[0];
|
|
args[arg_to_integ] = xmin + 2.5 * dx;
|
|
f1[2] = function (args);
|
|
args[arg_to_integ] = xmin + 3.5 * dx;
|
|
f1[3] = function (args);
|
|
f1[4] = f[1];
|
|
args[arg_to_integ] = xmin + 5.5 * dx;
|
|
f1[5] = function (args);
|
|
args[arg_to_integ] = xmin + 6.5 * dx;
|
|
f1[6] = function (args);
|
|
f1[7] = f[2];
|
|
args[arg_to_integ] = xmin + 8.5 * dx;
|
|
f1[8] = function (args);
|
|
|
|
DP I1_pre = 3.0 * dx * f[0];
|
|
DP I1 = dx * (f1[0] + f1[1] + f1[2]);
|
|
|
|
if (fabs(I1 - I1_pre) > req_prec || rec_level < 5)
|
|
I1 = Integrate_rec_main (function, args, arg_to_integ, xmin, xmin + 3.0 * dx,
|
|
req_prec, rec_level + 1, max_rec_level, f1);
|
|
|
|
DP I2_pre = 3.0 * dx * f[1];
|
|
DP I2 = dx * (f1[3] + f1[4] + f1[5]);
|
|
|
|
if (fabs(I2 - I2_pre) > req_prec || rec_level < 5)
|
|
I2 = Integrate_rec_main (function, args, arg_to_integ, xmin + 3.0 * dx, xmin + 6.0 * dx,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[3]);
|
|
|
|
DP I3_pre = 3.0 * dx * f[2];
|
|
DP I3 = dx * (f1[6] + f1[7] + f1[8]);
|
|
|
|
if (fabs(I3 - I3_pre) > req_prec || rec_level < 5)
|
|
I3 = Integrate_rec_main (function, args, arg_to_integ, xmin + 6.0 * dx, xmax,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[6]);
|
|
|
|
delete[] f1;
|
|
|
|
return(I1 + I2 + I3);
|
|
}
|
|
|
|
DP Integrate_rec (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax,
|
|
DP req_prec, int max_rec_level)
|
|
{
|
|
if (xmax < xmin) return(-Integrate_rec(function, args, arg_to_integ, xmax, xmin, req_prec, max_rec_level));
|
|
|
|
DP* f = new DP[3];
|
|
DP dx = (xmax - xmin)/3.0;
|
|
|
|
DP sum_fabs_f = 0.0;
|
|
|
|
for (int i = 0; i < 3; ++i) {
|
|
args[arg_to_integ] = xmin + (i + 0.5) * dx;
|
|
f[i] = function (args);
|
|
sum_fabs_f += fabs(f[i]);
|
|
}
|
|
|
|
DP req_prec_rec = sum_fabs_f > 1.0 ? sum_fabs_f * req_prec : req_prec;
|
|
|
|
DP answer = Integrate_rec_main (function, args, arg_to_integ, xmin, xmax, req_prec_rec, 0, max_rec_level, f);
|
|
|
|
delete[] f;
|
|
|
|
return(answer);
|
|
}
|
|
|
|
DP Integrate_rec_main_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
|
|
DP xmin, DP xmax, DP req_prec, int rec_level, int max_rec_level, DP* f)
|
|
{
|
|
if (rec_level > max_rec_level) {
|
|
return((xmax - xmin) * (f[0] + f[1] + f[2])/3.0);
|
|
}
|
|
|
|
DP* f1 = new DP[9];
|
|
|
|
DP dx = (xmax - xmin)/9.0;
|
|
|
|
args[arg_to_integ] = xmin + 0.5 * dx;
|
|
f1[0] = function (args, Itable);
|
|
f1[1] = f[0];
|
|
args[arg_to_integ] = xmin + 2.5 * dx;
|
|
f1[2] = function (args, Itable);
|
|
args[arg_to_integ] = xmin + 3.5 * dx;
|
|
f1[3] = function (args, Itable);
|
|
f1[4] = f[1];
|
|
args[arg_to_integ] = xmin + 5.5 * dx;
|
|
f1[5] = function (args, Itable);
|
|
args[arg_to_integ] = xmin + 6.5 * dx;
|
|
f1[6] = function (args, Itable);
|
|
f1[7] = f[2];
|
|
args[arg_to_integ] = xmin + 8.5 * dx;
|
|
f1[8] = function (args, Itable);
|
|
|
|
DP I1_pre = 3.0 * dx * f[0];
|
|
DP I1 = dx * (f1[0] + f1[1] + f1[2]);
|
|
|
|
if (fabs(I1 - I1_pre) > req_prec || rec_level < 5)
|
|
I1 = Integrate_rec_main_using_table (function, args, arg_to_integ, Itable, xmin, xmin + 3.0 * dx,
|
|
req_prec, rec_level + 1, max_rec_level, f1);
|
|
|
|
DP I2_pre = 3.0 * dx * f[1];
|
|
DP I2 = dx * (f1[3] + f1[4] + f1[5]);
|
|
|
|
if (fabs(I2 - I2_pre) > req_prec || rec_level < 5)
|
|
I2 = Integrate_rec_main_using_table (function, args, arg_to_integ, Itable, xmin + 3.0 * dx, xmin + 6.0 * dx,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[3]);
|
|
|
|
DP I3_pre = 3.0 * dx * f[2];
|
|
DP I3 = dx * (f1[6] + f1[7] + f1[8]);
|
|
|
|
if (fabs(I3 - I3_pre) > req_prec || rec_level < 5)
|
|
I3 = Integrate_rec_main_using_table (function, args, arg_to_integ, Itable, xmin + 6.0 * dx, xmax,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[6]);
|
|
|
|
delete[] f1;
|
|
|
|
return(I1 + I2 + I3);
|
|
}
|
|
|
|
DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
|
|
DP xmin, DP xmax, DP req_prec, int max_rec_level)
|
|
{
|
|
if (xmax < xmin) return(-Integrate_rec_using_table (function, args, arg_to_integ, Itable, xmax, xmin,
|
|
req_prec, max_rec_level));
|
|
|
|
DP* f = new DP[3];
|
|
DP dx = (xmax - xmin)/3.0;
|
|
|
|
DP sum_fabs_f = 0.0;
|
|
|
|
for (int i = 0; i < 3; ++i) {
|
|
args[arg_to_integ] = xmin + (i + 0.5) * dx;
|
|
f[i] = function (args, Itable);
|
|
sum_fabs_f += fabs(f[i]);
|
|
}
|
|
|
|
DP req_prec_rec = sum_fabs_f > 1.0 ? sum_fabs_f * req_prec : req_prec;
|
|
|
|
DP answer = Integrate_rec_main_using_table (function, args, arg_to_integ, Itable, xmin, xmax,
|
|
req_prec_rec, 0, max_rec_level, f);
|
|
|
|
delete[] f;
|
|
|
|
return(answer);
|
|
}
|
|
|
|
DP Integrate_rec_main_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
|
|
DP xmin, DP xmax, DP req_prec, int rec_level, int max_rec_level, DP* f, ofstream& outfile)
|
|
{
|
|
if (rec_level > max_rec_level) {
|
|
return((xmax - xmin) * (f[0] + f[1] + f[2])/3.0);
|
|
}
|
|
|
|
DP* f1 = new DP[9];
|
|
|
|
DP dx = (xmax - xmin)/9.0;
|
|
|
|
args[arg_to_integ] = xmin + 0.5 * dx;
|
|
f1[0] = function (args, Itable);
|
|
outfile << args << "\t" << f1[0] << endl;
|
|
f1[1] = f[0];
|
|
args[arg_to_integ] = xmin + 2.5 * dx;
|
|
f1[2] = function (args, Itable);
|
|
outfile << args << "\t" << f1[2] << endl;
|
|
args[arg_to_integ] = xmin + 3.5 * dx;
|
|
f1[3] = function (args, Itable);
|
|
outfile << args << "\t" << f1[3] << endl;
|
|
f1[4] = f[1];
|
|
args[arg_to_integ] = xmin + 5.5 * dx;
|
|
f1[5] = function (args, Itable);
|
|
outfile << args << "\t" << f1[5] << endl;
|
|
args[arg_to_integ] = xmin + 6.5 * dx;
|
|
f1[6] = function (args, Itable);
|
|
outfile << args << "\t" << f1[6] << endl;
|
|
f1[7] = f[2];
|
|
args[arg_to_integ] = xmin + 8.5 * dx;
|
|
f1[8] = function (args, Itable);
|
|
outfile << args << "\t" << f1[8] << endl;
|
|
|
|
DP I1_pre = 3.0 * dx * f[0];
|
|
DP I1 = dx * (f1[0] + f1[1] + f1[2]);
|
|
|
|
if (fabs(I1 - I1_pre) > req_prec || rec_level < 5)
|
|
I1 = Integrate_rec_main_using_table (function, args, arg_to_integ, Itable, xmin, xmin + 3.0 * dx,
|
|
req_prec, rec_level + 1, max_rec_level, f1, outfile);
|
|
|
|
DP I2_pre = 3.0 * dx * f[1];
|
|
DP I2 = dx * (f1[3] + f1[4] + f1[5]);
|
|
|
|
if (fabs(I2 - I2_pre) > req_prec || rec_level < 5)
|
|
I2 = Integrate_rec_main_using_table (function, args, arg_to_integ, Itable, xmin + 3.0 * dx, xmin + 6.0 * dx,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[3], outfile);
|
|
|
|
DP I3_pre = 3.0 * dx * f[2];
|
|
DP I3 = dx * (f1[6] + f1[7] + f1[8]);
|
|
|
|
if (fabs(I3 - I3_pre) > req_prec || rec_level < 5)
|
|
I3 = Integrate_rec_main_using_table (function, args, arg_to_integ, Itable, xmin + 6.0 * dx, xmax,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[6], outfile);
|
|
|
|
delete[] f1;
|
|
|
|
return(I1 + I2 + I3);
|
|
}
|
|
|
|
DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
|
|
DP xmin, DP xmax, DP req_prec, int max_rec_level, ofstream& outfile)
|
|
{
|
|
if (xmax < xmin) return(-Integrate_rec_using_table (function, args, arg_to_integ, Itable, xmax, xmin,
|
|
req_prec, max_rec_level, outfile));
|
|
|
|
DP* f = new DP[3];
|
|
DP dx = (xmax - xmin)/3.0;
|
|
|
|
DP sum_fabs_f = 0.0;
|
|
|
|
for (int i = 0; i < 3; ++i) {
|
|
args[arg_to_integ] = xmin + (i + 0.5) * dx;
|
|
f[i] = function (args, Itable);
|
|
outfile << args << "\t" << f[i] << endl;
|
|
sum_fabs_f += fabs(f[i]);
|
|
}
|
|
|
|
DP req_prec_rec = sum_fabs_f > 1.0 ? sum_fabs_f * req_prec : req_prec;
|
|
|
|
DP answer = Integrate_rec_main_using_table (function, args, arg_to_integ, Itable, xmin, xmax,
|
|
req_prec_rec, 0, max_rec_level, f, outfile);
|
|
|
|
delete[] f;
|
|
|
|
return(answer);
|
|
}
|
|
|
|
DP Integrate_rec_main_using_table_and_file (DP (*function) (Vect_DP, I_table, ofstream&),
|
|
Vect_DP& args, int arg_to_integ, I_table Itable,
|
|
DP xmin, DP xmax, DP req_prec, int rec_level, int max_rec_level,
|
|
DP* f, ofstream& outfile)
|
|
{
|
|
if (rec_level > max_rec_level) {
|
|
return((xmax - xmin) * (f[0] + f[1] + f[2])/3.0);
|
|
}
|
|
|
|
DP* f1 = new DP[9];
|
|
|
|
DP dx = (xmax - xmin)/9.0;
|
|
|
|
args[arg_to_integ] = xmin + 0.5 * dx;
|
|
f1[0] = function (args, Itable, outfile);
|
|
outfile << args << "\t" << f1[0] << endl;
|
|
f1[1] = f[0];
|
|
args[arg_to_integ] = xmin + 2.5 * dx;
|
|
f1[2] = function (args, Itable, outfile);
|
|
outfile << args << "\t" << f1[2] << endl;
|
|
args[arg_to_integ] = xmin + 3.5 * dx;
|
|
f1[3] = function (args, Itable, outfile);
|
|
outfile << args << "\t" << f1[3] << endl;
|
|
f1[4] = f[1];
|
|
args[arg_to_integ] = xmin + 5.5 * dx;
|
|
f1[5] = function (args, Itable, outfile);
|
|
outfile << args << "\t" << f1[5] << endl;
|
|
args[arg_to_integ] = xmin + 6.5 * dx;
|
|
f1[6] = function (args, Itable, outfile);
|
|
outfile << args << "\t" << f1[6] << endl;
|
|
f1[7] = f[2];
|
|
args[arg_to_integ] = xmin + 8.5 * dx;
|
|
f1[8] = function (args, Itable, outfile);
|
|
outfile << args << "\t" << f1[8] << endl;
|
|
|
|
DP I1_pre = 3.0 * dx * f[0];
|
|
DP I1 = dx * (f1[0] + f1[1] + f1[2]);
|
|
|
|
if (fabs(I1 - I1_pre) > req_prec || rec_level < 5)
|
|
I1 = Integrate_rec_main_using_table_and_file (function, args, arg_to_integ, Itable, xmin, xmin + 3.0 * dx,
|
|
req_prec, rec_level + 1, max_rec_level, f1, outfile);
|
|
|
|
DP I2_pre = 3.0 * dx * f[1];
|
|
DP I2 = dx * (f1[3] + f1[4] + f1[5]);
|
|
|
|
if (fabs(I2 - I2_pre) > req_prec || rec_level < 5)
|
|
I2 = Integrate_rec_main_using_table_and_file (function, args, arg_to_integ, Itable, xmin + 3.0 * dx, xmin + 6.0 * dx,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[3], outfile);
|
|
|
|
DP I3_pre = 3.0 * dx * f[2];
|
|
DP I3 = dx * (f1[6] + f1[7] + f1[8]);
|
|
|
|
if (fabs(I3 - I3_pre) > req_prec || rec_level < 5)
|
|
I3 = Integrate_rec_main_using_table_and_file (function, args, arg_to_integ, Itable, xmin + 6.0 * dx, xmax,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[6], outfile);
|
|
|
|
delete[] f1;
|
|
|
|
return(I1 + I2 + I3);
|
|
}
|
|
|
|
DP Integrate_rec_using_table_and_file (DP (*function) (Vect_DP, I_table, ofstream&), Vect_DP& args, int arg_to_integ,
|
|
I_table Itable, DP xmin, DP xmax, DP req_prec, int max_rec_level, ofstream& outfile)
|
|
{
|
|
if (xmax < xmin)
|
|
return(-Integrate_rec_using_table_and_file (function, args, arg_to_integ, Itable, xmax, xmin,
|
|
req_prec, max_rec_level, outfile));
|
|
|
|
DP* f = new DP[3];
|
|
DP dx = (xmax - xmin)/3.0;
|
|
|
|
DP sum_fabs_f = 0.0;
|
|
|
|
for (int i = 0; i < 3; ++i) {
|
|
args[arg_to_integ] = xmin + (i + 0.5) * dx;
|
|
f[i] = function (args, Itable, outfile);
|
|
outfile << args << "\t" << f[i] << endl;
|
|
sum_fabs_f += fabs(f[i]);
|
|
}
|
|
|
|
DP req_prec_rec = sum_fabs_f > 1.0 ? sum_fabs_f * req_prec : req_prec;
|
|
|
|
DP answer = Integrate_rec_main_using_table_and_file (function, args, arg_to_integ, Itable, xmin, xmax,
|
|
req_prec_rec, 0, max_rec_level, f, outfile);
|
|
|
|
delete[] f;
|
|
|
|
return(answer);
|
|
}
|
|
|
|
// THESE TWO FUNCTIONS HAVE NOT BEEN TESTED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
DP Integrate_exp_rec_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ,
|
|
I_table Itable, DP xmin, DP xmax, DP eta_old, DP req_prec, int rec_level,
|
|
int max_rec_level, DP* f, ofstream& outfile)
|
|
{
|
|
DP dx_over_x_old = 0.5 * (eta_old - 1.0/eta_old);
|
|
if (rec_level > max_rec_level) {
|
|
return(dx_over_x_old * xmin * eta_old * (f[0] + eta_old * (f[1] + eta_old * f[2])));
|
|
}
|
|
|
|
DP* f1 = new DP[9];
|
|
DP* x_pts = new DP[9];
|
|
DP eta = pow(xmax/xmin, 1.0/18.0);
|
|
DP eta_sq = eta * eta;
|
|
DP dx_over_x = 0.5 * (eta - 1.0/eta);
|
|
|
|
x_pts[0] = xmin * eta;
|
|
for (int i = 1; i < 9; ++i) x_pts[i] = x_pts[i-1] * eta_sq;
|
|
|
|
args[arg_to_integ] = x_pts[0];
|
|
f1[0] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f1[0] << endl;
|
|
f1[1] = f[0];
|
|
args[arg_to_integ] = x_pts[2];
|
|
f1[2] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f1[2] << endl;
|
|
args[arg_to_integ] = x_pts[3];
|
|
f1[3] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f1[3] << endl;
|
|
f1[4] = f[1];
|
|
args[arg_to_integ] = x_pts[5];
|
|
f1[5] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f1[5] << endl;
|
|
args[arg_to_integ] = x_pts[6];
|
|
f1[6] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f1[6] << endl;
|
|
f1[7] = f[2];
|
|
args[arg_to_integ] = x_pts[8];
|
|
f1[8] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f1[8] << endl;
|
|
|
|
DP I1_pre = dx_over_x_old * x_pts[1] * f1[1];
|
|
DP I1 = dx_over_x * (x_pts[0] * f1[0] + x_pts[1] * f1[1] + x_pts[2] * f1[2]);
|
|
|
|
if (fabs(I1 - I1_pre) > req_prec || rec_level < 5)
|
|
I1 = Integrate_exp_rec_using_table (function, args, arg_to_integ, Itable, xmin, x_pts[2] * eta, eta,
|
|
req_prec, rec_level + 1, max_rec_level, f1, outfile);
|
|
|
|
DP I2_pre = dx_over_x_old * x_pts[4] * f1[4];
|
|
DP I2 = dx_over_x * (x_pts[3] * f1[3] + x_pts[4] * f1[4] + x_pts[5] * f1[5]);
|
|
|
|
if (fabs(I2 - I2_pre) > req_prec || rec_level < 5)
|
|
I2 = Integrate_exp_rec_using_table (function, args, arg_to_integ, Itable, x_pts[2] * eta, x_pts[5] * eta, eta,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[3], outfile);
|
|
|
|
DP I3_pre = dx_over_x_old * f1[7];
|
|
DP I3 = dx_over_x * (x_pts[6] * f1[6] + x_pts[7] * f1[7] + x_pts[8] * f1[8]);
|
|
|
|
if (fabs(I3 - I3_pre) > req_prec || rec_level < 5)
|
|
I3 = Integrate_exp_rec_using_table (function, args, arg_to_integ, Itable, x_pts[5] * eta, xmax, eta,
|
|
req_prec, rec_level + 1, max_rec_level, &f1[6], outfile);
|
|
|
|
delete[] f1;
|
|
|
|
return(I1 + I2 + I3);
|
|
}
|
|
|
|
DP Integrate_exp_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
|
|
DP xmin, DP xmax, DP req_prec, int max_rec_level, ofstream& outfile)
|
|
{
|
|
// This uses an exponential progression of sampling points.
|
|
|
|
if (xmin <= 0.0) ABACUSerror("Use xmin > 0.0 in Integrate_exp.");
|
|
if (xmax <= 0.0) ABACUSerror("Use xmax > 0.0 in Integrate_exp.");
|
|
|
|
if (xmax < xmin) return(-Integrate_exp_using_table (function, args, arg_to_integ, Itable, xmax, xmin,
|
|
req_prec, max_rec_level, outfile));
|
|
|
|
DP* f = new DP[3];
|
|
DP eta = pow(xmax/xmin, 1.0/6.0);
|
|
|
|
DP sum_fabs_f = 0.0;
|
|
|
|
args[arg_to_integ] = xmin * eta;
|
|
f[0] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f[0] << endl;
|
|
sum_fabs_f += fabs(f[0]);
|
|
args[arg_to_integ] *= eta * eta;
|
|
f[1] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f[1] << endl;
|
|
sum_fabs_f += fabs(f[1]);
|
|
args[arg_to_integ] *= eta * eta;
|
|
f[2] = function (args, Itable);
|
|
outfile << args[arg_to_integ] << "\t" << f[2] << endl;
|
|
sum_fabs_f += fabs(f[2]);
|
|
|
|
DP req_prec_rec = sum_fabs_f > 1.0 ? sum_fabs_f * req_prec : req_prec;
|
|
|
|
DP answer = Integrate_exp_rec_using_table (function, args, arg_to_integ, Itable, xmin, xmax, eta,
|
|
req_prec_rec, 0, max_rec_level, f, outfile);
|
|
|
|
delete[] f;
|
|
|
|
return(answer);
|
|
}
|
|
|
|
|
|
|
|
// *********************************** Adaptive Riemann sum (better implementation) *******************************************
|
|
|
|
std::ostream& operator<< (std::ostream& s, const Integral_result& res)
|
|
{
|
|
s << res.integ_est << "\t" << res.abs_prec << "\t" << res.rel_prec << "\t" << res.n_vals;
|
|
|
|
return(s);
|
|
}
|
|
|
|
Integral_data::~Integral_data()
|
|
{
|
|
if (data != 0) delete[] data;
|
|
if (abs_d2f_dx != 0) delete[] abs_d2f_dx;
|
|
}
|
|
|
|
Integral_data::Integral_data (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin_ref, DP xmax_ref)
|
|
{
|
|
xmin = xmin_ref;
|
|
xmax = xmax_ref;
|
|
integ_res.n_vals = 9;
|
|
|
|
data = new data_pt[integ_res.n_vals];
|
|
abs_d2f_dx = new DP[integ_res.n_vals/3];
|
|
|
|
integ_res.integ_est = 0.0;
|
|
|
|
// Initialize x, f and dx:
|
|
DP dx = (xmax - xmin)/9.0;
|
|
for (int i = 0; i < 9; ++i) {
|
|
data[i].x = xmin + (i + 0.5) * dx;
|
|
args[arg_to_integ] = data[i].x;
|
|
data[i].f = function (args);
|
|
data[i].dx = dx;
|
|
integ_res.integ_est += data[i].dx * data[i].f;
|
|
}
|
|
|
|
max_abs_d2f_dx = 0.0;
|
|
integ_res.abs_prec = 0.0;
|
|
for (int j = 0; j < 3; ++j) {
|
|
abs_d2f_dx[j] = dx * fabs(data[3*j].f - 2.0 * data[3*j + 1].f + data[3*j + 2].f);
|
|
max_abs_d2f_dx = ABACUS::max(max_abs_d2f_dx, abs_d2f_dx[j]);
|
|
integ_res.abs_prec += abs_d2f_dx[j];
|
|
}
|
|
|
|
integ_res.rel_prec = integ_res.abs_prec/integ_res.integ_est;
|
|
}
|
|
|
|
Integral_data::Integral_data (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ,
|
|
I_table Itable, DP xmin_ref, DP xmax_ref)
|
|
{
|
|
xmin = xmin_ref;
|
|
xmax = xmax_ref;
|
|
integ_res.n_vals = 9;
|
|
|
|
data = new data_pt[integ_res.n_vals];
|
|
abs_d2f_dx = new DP[integ_res.n_vals/3];
|
|
|
|
integ_res.integ_est = 0.0;
|
|
|
|
// Initialize x, f and dx:
|
|
DP dx = (xmax - xmin)/9.0;
|
|
for (int i = 0; i < 9; ++i) {
|
|
data[i].x = xmin + (i + 0.5) * dx;
|
|
args[arg_to_integ] = data[i].x;
|
|
data[i].f = function (args, Itable);
|
|
data[i].dx = dx;
|
|
integ_res.integ_est += data[i].dx * data[i].f;
|
|
}
|
|
|
|
max_abs_d2f_dx = 0.0;
|
|
integ_res.abs_prec = 0.0;
|
|
for (int j = 0; j < 3; ++j) {
|
|
abs_d2f_dx[j] = dx * fabs(data[3*j].f - 2.0 * data[3*j + 1].f + data[3*j + 2].f);
|
|
max_abs_d2f_dx = ABACUS::max(max_abs_d2f_dx, abs_d2f_dx[j]);
|
|
integ_res.abs_prec += abs_d2f_dx[j];
|
|
}
|
|
|
|
integ_res.rel_prec = integ_res.abs_prec/integ_res.integ_est;
|
|
}
|
|
|
|
Integral_data::Integral_data (DP (*function) (Vect_DP, Integral_table), Vect_DP& args, int arg_to_integ,
|
|
Integral_table Itable, DP xmin_ref, DP xmax_ref)
|
|
{
|
|
xmin = xmin_ref;
|
|
xmax = xmax_ref;
|
|
integ_res.n_vals = 9;
|
|
|
|
data = new data_pt[integ_res.n_vals];
|
|
abs_d2f_dx = new DP[integ_res.n_vals/3];
|
|
|
|
integ_res.integ_est = 0.0;
|
|
|
|
// Initialize x, f and dx:
|
|
DP dx = (xmax - xmin)/9.0;
|
|
for (int i = 0; i < 9; ++i) {
|
|
data[i].x = xmin + (i + 0.5) * dx;
|
|
args[arg_to_integ] = data[i].x;
|
|
data[i].f = function (args, Itable);
|
|
data[i].dx = dx;
|
|
integ_res.integ_est += data[i].dx * data[i].f;
|
|
}
|
|
|
|
max_abs_d2f_dx = 0.0;
|
|
integ_res.abs_prec = 0.0;
|
|
for (int j = 0; j < 3; ++j) {
|
|
abs_d2f_dx[j] = dx * fabs(data[3*j].f - 2.0 * data[3*j + 1].f + data[3*j + 2].f);
|
|
max_abs_d2f_dx = ABACUS::max(max_abs_d2f_dx, abs_d2f_dx[j]);
|
|
integ_res.abs_prec += abs_d2f_dx[j];
|
|
}
|
|
|
|
integ_res.rel_prec = integ_res.abs_prec/integ_res.integ_est;
|
|
}
|
|
|
|
void Integral_data::Save (ofstream& outfile)
|
|
{
|
|
// Reset file writing position
|
|
outfile.seekp(0);
|
|
int i = 0;
|
|
while (i < integ_res.n_vals && !(data[i].x == 0.0 && data[i].f == 0.0)) {
|
|
outfile << setprecision(16) << data[i].x << "\t" << data[i].f << endl;
|
|
i++;
|
|
}
|
|
}
|
|
|
|
void Integral_data::Improve_estimate (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, int max_nr_pts)
|
|
{
|
|
// We generate a max of 3* n_vals points
|
|
data_pt* new_data = new data_pt[3* integ_res.n_vals];
|
|
DP* new_abs_d2f_dx = new DP[integ_res.n_vals];
|
|
|
|
// Check points in batches of 3; if needed, improve
|
|
int threei = 0;
|
|
int index_new = 0;
|
|
int i, j;
|
|
|
|
integ_res.abs_prec = 0.0;
|
|
integ_res.integ_est = 0.0;
|
|
|
|
DP new_max_abs_d2f_dx = 0.0;
|
|
|
|
for (i = 0; i < integ_res.n_vals/3; ++i) {
|
|
|
|
threei = 3 * i;
|
|
|
|
if (abs_d2f_dx[i] <= 0.1 * max_abs_d2f_dx || index_new + integ_res.n_vals - threei > max_nr_pts) {
|
|
|
|
// simply transfer the data points into new_data
|
|
new_abs_d2f_dx[index_new/3] = abs_d2f_dx[i];
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3]);
|
|
for (j = 0; j < 3; ++j) {
|
|
new_data[index_new].x = data[threei + j].x;
|
|
new_data[index_new].f = data[threei + j].f;
|
|
new_data[index_new].dx = data[threei + j].dx;
|
|
integ_res.integ_est += new_data[index_new].dx * new_data[index_new].f;
|
|
index_new++;
|
|
}
|
|
integ_res.abs_prec += abs_d2f_dx[i];
|
|
|
|
}
|
|
|
|
else { // create six new entries and transfer the three existing ones
|
|
|
|
new_data[index_new].dx = data[threei].dx/3.0;
|
|
for (j = 1; j < 9; ++j) new_data[index_new + j].dx = new_data[index_new].dx;
|
|
|
|
new_data[index_new].x = data[threei].x - new_data[index_new].dx;
|
|
new_data[index_new + 1].x = data[threei].x;
|
|
new_data[index_new + 2].x = data[threei].x + new_data[index_new].dx;
|
|
|
|
new_data[index_new + 3].x = data[threei + 1].x - new_data[index_new].dx;
|
|
new_data[index_new + 4].x = data[threei + 1].x;
|
|
new_data[index_new + 5].x = data[threei + 1].x + new_data[index_new].dx;
|
|
|
|
new_data[index_new + 6].x = data[threei + 2].x - new_data[index_new].dx;
|
|
new_data[index_new + 7].x = data[threei + 2].x;
|
|
new_data[index_new + 8].x = data[threei + 2].x + new_data[index_new].dx;
|
|
|
|
args[arg_to_integ] = new_data[index_new].x;
|
|
new_data[index_new].f = function(args);
|
|
new_data[index_new + 1].f = data[threei].f;
|
|
args[arg_to_integ] = new_data[index_new + 2].x;
|
|
new_data[index_new + 2].f = function(args);
|
|
|
|
args[arg_to_integ] = new_data[index_new + 3].x;
|
|
new_data[index_new + 3].f = function(args);
|
|
new_data[index_new + 4].f = data[threei + 1].f;
|
|
args[arg_to_integ] = new_data[index_new + 5].x;
|
|
new_data[index_new + 5].f = function(args);
|
|
|
|
args[arg_to_integ] = new_data[index_new + 6].x;
|
|
new_data[index_new + 6].f = function(args);
|
|
new_data[index_new + 7].f = data[threei + 2].f;
|
|
args[arg_to_integ] = new_data[index_new + 8].x;
|
|
new_data[index_new + 8].f = function(args);
|
|
|
|
new_abs_d2f_dx[index_new/3] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new].f - 2.0 * new_data[index_new + 1].f + new_data[index_new + 2].f));
|
|
new_abs_d2f_dx[index_new/3 + 1] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new + 3].f - 2.0 * new_data[index_new + 4].f + new_data[index_new + 5].f));
|
|
new_abs_d2f_dx[index_new/3 + 2] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new + 6].f - 2.0 * new_data[index_new + 7].f + new_data[index_new + 8].f));
|
|
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3]);
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3 + 1]);
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3 + 2]);
|
|
|
|
integ_res.integ_est += new_data[index_new].dx * (new_data[index_new].f + new_data[index_new + 1].f + new_data[index_new + 2].f
|
|
+ new_data[index_new + 3].f + new_data[index_new + 4].f + new_data[index_new + 5].f
|
|
+ new_data[index_new + 6].f + new_data[index_new + 7].f + new_data[index_new + 8].f);
|
|
|
|
integ_res.abs_prec += new_abs_d2f_dx[index_new/3] + new_abs_d2f_dx[index_new/3 + 1] + new_abs_d2f_dx[index_new/3 + 2];
|
|
|
|
index_new += 9;
|
|
|
|
} // else
|
|
|
|
} // for (i < nvals/3)
|
|
|
|
integ_res.rel_prec = integ_res.abs_prec/integ_res.integ_est;
|
|
|
|
integ_res.n_vals = index_new;
|
|
|
|
delete[] data;
|
|
data = new_data;
|
|
|
|
max_abs_d2f_dx = new_max_abs_d2f_dx;
|
|
|
|
delete[] abs_d2f_dx;
|
|
abs_d2f_dx = new_abs_d2f_dx;
|
|
|
|
return;
|
|
}
|
|
|
|
void Integral_data::Improve_estimate (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ,
|
|
I_table Itable, int max_nr_pts)
|
|
{
|
|
// We generate a max of 3* n_vals points
|
|
data_pt* new_data = new data_pt[3* integ_res.n_vals];
|
|
DP* new_abs_d2f_dx = new DP[integ_res.n_vals];
|
|
|
|
// Check points in batches of 3; if needed, improve
|
|
int threei = 0;
|
|
int index_new = 0;
|
|
int i, j;
|
|
|
|
integ_res.abs_prec = 0.0;
|
|
integ_res.integ_est = 0.0;
|
|
|
|
DP new_max_abs_d2f_dx = 0.0;
|
|
|
|
for (i = 0; i < integ_res.n_vals/3; ++i) {
|
|
|
|
threei = 3 * i;
|
|
|
|
if (abs_d2f_dx[i] <= 0.1 * max_abs_d2f_dx || index_new + integ_res.n_vals - threei > max_nr_pts) {
|
|
|
|
// simply transfer the data points into new_data
|
|
new_abs_d2f_dx[index_new/3] = abs_d2f_dx[i];
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3]);
|
|
for (j = 0; j < 3; ++j) {
|
|
new_data[index_new].x = data[threei + j].x;
|
|
new_data[index_new].f = data[threei + j].f;
|
|
new_data[index_new].dx = data[threei + j].dx;
|
|
integ_res.integ_est += new_data[index_new].dx * new_data[index_new].f;
|
|
index_new++;
|
|
}
|
|
integ_res.abs_prec += abs_d2f_dx[i];
|
|
|
|
}
|
|
|
|
else { // create six new entries and transfer the three existing ones
|
|
|
|
new_data[index_new].dx = data[threei].dx/3.0;
|
|
for (j = 1; j < 9; ++j) new_data[index_new + j].dx = new_data[index_new].dx;
|
|
|
|
new_data[index_new].x = data[threei].x - new_data[index_new].dx;
|
|
new_data[index_new + 1].x = data[threei].x;
|
|
new_data[index_new + 2].x = data[threei].x + new_data[index_new].dx;
|
|
|
|
new_data[index_new + 3].x = data[threei + 1].x - new_data[index_new].dx;
|
|
new_data[index_new + 4].x = data[threei + 1].x;
|
|
new_data[index_new + 5].x = data[threei + 1].x + new_data[index_new].dx;
|
|
|
|
new_data[index_new + 6].x = data[threei + 2].x - new_data[index_new].dx;
|
|
new_data[index_new + 7].x = data[threei + 2].x;
|
|
new_data[index_new + 8].x = data[threei + 2].x + new_data[index_new].dx;
|
|
|
|
args[arg_to_integ] = new_data[index_new].x;
|
|
new_data[index_new].f = function(args, Itable);
|
|
new_data[index_new + 1].f = data[threei].f;
|
|
args[arg_to_integ] = new_data[index_new + 2].x;
|
|
new_data[index_new + 2].f = function(args, Itable);
|
|
|
|
args[arg_to_integ] = new_data[index_new + 3].x;
|
|
new_data[index_new + 3].f = function(args, Itable);
|
|
new_data[index_new + 4].f = data[threei + 1].f;
|
|
args[arg_to_integ] = new_data[index_new + 5].x;
|
|
new_data[index_new + 5].f = function(args, Itable);
|
|
|
|
args[arg_to_integ] = new_data[index_new + 6].x;
|
|
new_data[index_new + 6].f = function(args, Itable);
|
|
new_data[index_new + 7].f = data[threei + 2].f;
|
|
args[arg_to_integ] = new_data[index_new + 8].x;
|
|
new_data[index_new + 8].f = function(args, Itable);
|
|
|
|
new_abs_d2f_dx[index_new/3] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new].f - 2.0 * new_data[index_new + 1].f + new_data[index_new + 2].f));
|
|
new_abs_d2f_dx[index_new/3 + 1] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new + 3].f - 2.0 * new_data[index_new + 4].f + new_data[index_new + 5].f));
|
|
new_abs_d2f_dx[index_new/3 + 2] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new + 6].f - 2.0 * new_data[index_new + 7].f + new_data[index_new + 8].f));
|
|
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3]);
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3 + 1]);
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3 + 2]);
|
|
|
|
integ_res.integ_est += new_data[index_new].dx * (new_data[index_new].f + new_data[index_new + 1].f + new_data[index_new + 2].f
|
|
+ new_data[index_new + 3].f + new_data[index_new + 4].f + new_data[index_new + 5].f
|
|
+ new_data[index_new + 6].f + new_data[index_new + 7].f + new_data[index_new + 8].f);
|
|
|
|
integ_res.abs_prec += new_abs_d2f_dx[index_new/3] + new_abs_d2f_dx[index_new/3 + 1] + new_abs_d2f_dx[index_new/3 + 2];
|
|
|
|
index_new += 9;
|
|
|
|
} // else
|
|
|
|
} // for (i < nvals/3)
|
|
|
|
integ_res.rel_prec = integ_res.abs_prec/integ_res.integ_est;
|
|
|
|
integ_res.n_vals = index_new;
|
|
|
|
delete[] data;
|
|
data = new_data;
|
|
|
|
max_abs_d2f_dx = new_max_abs_d2f_dx;
|
|
|
|
delete[] abs_d2f_dx;
|
|
abs_d2f_dx = new_abs_d2f_dx;
|
|
|
|
return;
|
|
}
|
|
|
|
void Integral_data::Improve_estimate (DP (*function) (Vect_DP, Integral_table), Vect_DP& args, int arg_to_integ,
|
|
Integral_table Itable, int max_nr_pts)
|
|
{
|
|
// We generate a max of 3* n_vals points
|
|
data_pt* new_data = new data_pt[3* integ_res.n_vals];
|
|
DP* new_abs_d2f_dx = new DP[integ_res.n_vals];
|
|
|
|
// Check points in batches of 3; if needed, improve
|
|
int threei = 0;
|
|
int index_new = 0;
|
|
int i, j;
|
|
|
|
integ_res.abs_prec = 0.0;
|
|
integ_res.integ_est = 0.0;
|
|
|
|
DP new_max_abs_d2f_dx = 0.0;
|
|
|
|
for (i = 0; i < integ_res.n_vals/3; ++i) {
|
|
|
|
threei = 3 * i;
|
|
|
|
if (abs_d2f_dx[i] <= 0.1 * max_abs_d2f_dx || index_new + integ_res.n_vals - threei > max_nr_pts) {
|
|
|
|
// simply transfer the data points into new_data
|
|
new_abs_d2f_dx[index_new/3] = abs_d2f_dx[i];
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3]);
|
|
for (j = 0; j < 3; ++j) {
|
|
new_data[index_new].x = data[threei + j].x;
|
|
new_data[index_new].f = data[threei + j].f;
|
|
new_data[index_new].dx = data[threei + j].dx;
|
|
integ_res.integ_est += new_data[index_new].dx * new_data[index_new].f;
|
|
index_new++;
|
|
}
|
|
integ_res.abs_prec += abs_d2f_dx[i];
|
|
|
|
}
|
|
|
|
else { // create six new entries and transfer the three existing ones
|
|
|
|
new_data[index_new].dx = data[threei].dx/3.0;
|
|
for (j = 1; j < 9; ++j) new_data[index_new + j].dx = new_data[index_new].dx;
|
|
|
|
new_data[index_new].x = data[threei].x - new_data[index_new].dx;
|
|
new_data[index_new + 1].x = data[threei].x;
|
|
new_data[index_new + 2].x = data[threei].x + new_data[index_new].dx;
|
|
|
|
new_data[index_new + 3].x = data[threei + 1].x - new_data[index_new].dx;
|
|
new_data[index_new + 4].x = data[threei + 1].x;
|
|
new_data[index_new + 5].x = data[threei + 1].x + new_data[index_new].dx;
|
|
|
|
new_data[index_new + 6].x = data[threei + 2].x - new_data[index_new].dx;
|
|
new_data[index_new + 7].x = data[threei + 2].x;
|
|
new_data[index_new + 8].x = data[threei + 2].x + new_data[index_new].dx;
|
|
|
|
args[arg_to_integ] = new_data[index_new].x;
|
|
new_data[index_new].f = function(args, Itable);
|
|
new_data[index_new + 1].f = data[threei].f;
|
|
args[arg_to_integ] = new_data[index_new + 2].x;
|
|
new_data[index_new + 2].f = function(args, Itable);
|
|
|
|
args[arg_to_integ] = new_data[index_new + 3].x;
|
|
new_data[index_new + 3].f = function(args, Itable);
|
|
new_data[index_new + 4].f = data[threei + 1].f;
|
|
args[arg_to_integ] = new_data[index_new + 5].x;
|
|
new_data[index_new + 5].f = function(args, Itable);
|
|
|
|
args[arg_to_integ] = new_data[index_new + 6].x;
|
|
new_data[index_new + 6].f = function(args, Itable);
|
|
new_data[index_new + 7].f = data[threei + 2].f;
|
|
args[arg_to_integ] = new_data[index_new + 8].x;
|
|
new_data[index_new + 8].f = function(args, Itable);
|
|
|
|
new_abs_d2f_dx[index_new/3] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new].f - 2.0 * new_data[index_new + 1].f + new_data[index_new + 2].f));
|
|
new_abs_d2f_dx[index_new/3 + 1] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new + 3].f - 2.0 * new_data[index_new + 4].f + new_data[index_new + 5].f));
|
|
new_abs_d2f_dx[index_new/3 + 2] =
|
|
fabs(new_data[index_new].dx * (new_data[index_new + 6].f - 2.0 * new_data[index_new + 7].f + new_data[index_new + 8].f));
|
|
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3]);
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3 + 1]);
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3 + 2]);
|
|
|
|
integ_res.integ_est += new_data[index_new].dx * (new_data[index_new].f + new_data[index_new + 1].f + new_data[index_new + 2].f
|
|
+ new_data[index_new + 3].f + new_data[index_new + 4].f + new_data[index_new + 5].f
|
|
+ new_data[index_new + 6].f + new_data[index_new + 7].f + new_data[index_new + 8].f);
|
|
|
|
integ_res.abs_prec += new_abs_d2f_dx[index_new/3] + new_abs_d2f_dx[index_new/3 + 1] + new_abs_d2f_dx[index_new/3 + 2];
|
|
|
|
index_new += 9;
|
|
|
|
} // else
|
|
|
|
} // for (i < nvals/3)
|
|
|
|
integ_res.rel_prec = integ_res.abs_prec/integ_res.integ_est;
|
|
|
|
integ_res.n_vals = index_new;
|
|
|
|
delete[] data;
|
|
data = new_data;
|
|
|
|
max_abs_d2f_dx = new_max_abs_d2f_dx;
|
|
|
|
delete[] abs_d2f_dx;
|
|
abs_d2f_dx = new_abs_d2f_dx;
|
|
|
|
return;
|
|
}
|
|
|
|
Integral_result Integrate_optimal (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax,
|
|
DP req_rel_prec, DP req_abs_prec, int max_nr_pts)
|
|
{
|
|
if (xmax < xmin) ABACUSerror("Use xmax > xmin in Integrate.");
|
|
|
|
Integral_data integ_dat (function, args, arg_to_integ, xmin, xmax);
|
|
|
|
while ((integ_dat.integ_res.rel_prec > req_rel_prec || integ_dat.integ_res.abs_prec > req_abs_prec)
|
|
&& integ_dat.integ_res.n_vals < max_nr_pts) {
|
|
integ_dat.Improve_estimate (function, args, arg_to_integ, max_nr_pts);
|
|
}
|
|
|
|
return(integ_dat.integ_res);
|
|
}
|
|
|
|
Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ,
|
|
I_table Itable, DP xmin, DP xmax,
|
|
DP req_rel_prec, DP req_abs_prec, int max_nr_pts)
|
|
{
|
|
if (xmax < xmin) ABACUSerror("Use xmax > xmin in Integrate.");
|
|
|
|
Integral_data integ_dat (function, args, arg_to_integ, Itable, xmin, xmax);
|
|
|
|
while ((integ_dat.integ_res.rel_prec > req_rel_prec || integ_dat.integ_res.abs_prec > req_abs_prec)
|
|
&& integ_dat.integ_res.n_vals < max_nr_pts) {
|
|
integ_dat.Improve_estimate (function, args, arg_to_integ, Itable, max_nr_pts);
|
|
}
|
|
|
|
return(integ_dat.integ_res);
|
|
}
|
|
|
|
Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, Integral_table), Vect_DP& args, int arg_to_integ,
|
|
Integral_table Itable, DP xmin, DP xmax,
|
|
DP req_rel_prec, DP req_abs_prec, int max_nr_pts)
|
|
{
|
|
if (xmax < xmin) ABACUSerror("Use xmax > xmin in Integrate.");
|
|
|
|
Integral_data integ_dat (function, args, arg_to_integ, Itable, xmin, xmax);
|
|
|
|
while ((integ_dat.integ_res.rel_prec > req_rel_prec || integ_dat.integ_res.abs_prec > req_abs_prec)
|
|
&& integ_dat.integ_res.n_vals < max_nr_pts) {
|
|
|
|
integ_dat.Improve_estimate (function, args, arg_to_integ, Itable, max_nr_pts);
|
|
}
|
|
|
|
return(integ_dat.integ_res);
|
|
}
|
|
|
|
Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ,
|
|
I_table Itable, DP xmin, DP xmax,
|
|
DP req_rel_prec, DP req_abs_prec, int max_nr_pts, ofstream& outfile)
|
|
{
|
|
if (xmax < xmin) ABACUSerror("Use xmax > xmin in Integrate.");
|
|
|
|
Integral_data integ_dat (function, args, arg_to_integ, Itable, xmin, xmax);
|
|
|
|
while ((integ_dat.integ_res.rel_prec > req_rel_prec || integ_dat.integ_res.abs_prec > req_abs_prec)
|
|
&& integ_dat.integ_res.n_vals < max_nr_pts) {
|
|
|
|
integ_dat.Improve_estimate (function, args, arg_to_integ, Itable, max_nr_pts);
|
|
|
|
integ_dat.Save (outfile);
|
|
|
|
}
|
|
|
|
return(integ_dat.integ_res);
|
|
}
|
|
|
|
|
|
// FOR COMPLEX FUNCTIONS:
|
|
// *********************************** Adaptive Riemann sum (better implementation) *******************************************
|
|
|
|
std::ostream& operator<< (std::ostream& s, const Integral_result_CX& res)
|
|
{
|
|
s << real(res.integ_est) << "\t" << imag(res.integ_est) << "\t"
|
|
<< res.abs_prec << "\t" << res.rel_prec << "\t" << res.n_vals;
|
|
|
|
return(s);
|
|
}
|
|
|
|
Integral_data_CX::~Integral_data_CX()
|
|
{
|
|
if (data != 0) delete[] data;
|
|
if (abs_d2f_dx != 0) delete[] abs_d2f_dx;
|
|
}
|
|
|
|
Integral_data_CX::Integral_data_CX (complex<DP> (*function) (Vect_DP), Vect_DP& args, int arg_to_integ,
|
|
DP xmin_ref, DP xmax_ref)
|
|
{
|
|
xmin = xmin_ref;
|
|
xmax = xmax_ref;
|
|
integ_res.n_vals = 9;
|
|
|
|
data = new data_pt_CX[integ_res.n_vals];
|
|
abs_d2f_dx = new DP[integ_res.n_vals/3];
|
|
|
|
integ_res.integ_est = 0.0;
|
|
|
|
// Initialize x, f and dx:
|
|
DP dx = (xmax - xmin)/9.0;
|
|
for (int i = 0; i < 9; ++i) {
|
|
data[i].x = xmin + (i + 0.5) * dx;
|
|
args[arg_to_integ] = data[i].x;
|
|
data[i].f = function (args);
|
|
data[i].dx = dx;
|
|
integ_res.integ_est += data[i].dx * data[i].f;
|
|
}
|
|
|
|
max_abs_d2f_dx = 0.0;
|
|
integ_res.abs_prec = 0.0;
|
|
for (int j = 0; j < 3; ++j) {
|
|
abs_d2f_dx[j] = dx * abs(data[3*j].f - 2.0 * data[3*j + 1].f + data[3*j + 2].f);
|
|
max_abs_d2f_dx = ABACUS::max(max_abs_d2f_dx, abs_d2f_dx[j]);
|
|
integ_res.abs_prec += abs_d2f_dx[j];
|
|
}
|
|
|
|
integ_res.rel_prec = integ_res.abs_prec/abs(integ_res.integ_est);
|
|
}
|
|
|
|
// No implementation with I_table yet for complex...
|
|
void Integral_data_CX::Save (ofstream& outfile)
|
|
{
|
|
// Reset file writing position
|
|
outfile.seekp(0);
|
|
int i = 0;
|
|
while (i < integ_res.n_vals && !(data[i].x == 0.0 && data[i].f == 0.0)) {
|
|
outfile << setprecision(16) << data[i].x << "\t" << data[i].f << endl;
|
|
i++;
|
|
}
|
|
}
|
|
|
|
void Integral_data_CX::Improve_estimate (complex<DP> (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, int max_nr_pts)
|
|
{
|
|
// We generate a max of 3* n_vals points
|
|
data_pt_CX* new_data = new data_pt_CX[3* integ_res.n_vals];
|
|
DP* new_abs_d2f_dx = new DP[integ_res.n_vals];
|
|
|
|
// Check points in batches of 3; if needed, improve
|
|
int threei = 0;
|
|
int index_new = 0;
|
|
int i, j;
|
|
|
|
integ_res.abs_prec = 0.0;
|
|
integ_res.integ_est = 0.0;
|
|
|
|
DP new_max_abs_d2f_dx = 0.0;
|
|
|
|
for (i = 0; i < integ_res.n_vals/3; ++i) {
|
|
|
|
threei = 3 * i;
|
|
|
|
if (abs_d2f_dx[i] <= 0.1 * max_abs_d2f_dx || index_new + integ_res.n_vals - threei > max_nr_pts) {
|
|
|
|
// simply transfer the data points into new_data
|
|
new_abs_d2f_dx[index_new/3] = abs_d2f_dx[i];
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3]);
|
|
for (j = 0; j < 3; ++j) {
|
|
new_data[index_new].x = data[threei + j].x;
|
|
new_data[index_new].f = data[threei + j].f;
|
|
new_data[index_new].dx = data[threei + j].dx;
|
|
integ_res.integ_est += new_data[index_new].dx * new_data[index_new].f;
|
|
index_new++;
|
|
}
|
|
integ_res.abs_prec += abs_d2f_dx[i];
|
|
|
|
}
|
|
|
|
else { // create six new entries and transfer the three existing ones
|
|
|
|
new_data[index_new].dx = data[threei].dx/3.0;
|
|
for (j = 1; j < 9; ++j) new_data[index_new + j].dx = new_data[index_new].dx;
|
|
|
|
new_data[index_new].x = data[threei].x - new_data[index_new].dx;
|
|
new_data[index_new + 1].x = data[threei].x;
|
|
new_data[index_new + 2].x = data[threei].x + new_data[index_new].dx;
|
|
|
|
new_data[index_new + 3].x = data[threei + 1].x - new_data[index_new].dx;
|
|
new_data[index_new + 4].x = data[threei + 1].x;
|
|
new_data[index_new + 5].x = data[threei + 1].x + new_data[index_new].dx;
|
|
|
|
new_data[index_new + 6].x = data[threei + 2].x - new_data[index_new].dx;
|
|
new_data[index_new + 7].x = data[threei + 2].x;
|
|
new_data[index_new + 8].x = data[threei + 2].x + new_data[index_new].dx;
|
|
|
|
args[arg_to_integ] = new_data[index_new].x;
|
|
new_data[index_new].f = function(args);
|
|
new_data[index_new + 1].f = data[threei].f;
|
|
args[arg_to_integ] = new_data[index_new + 2].x;
|
|
new_data[index_new + 2].f = function(args);
|
|
|
|
args[arg_to_integ] = new_data[index_new + 3].x;
|
|
new_data[index_new + 3].f = function(args);
|
|
new_data[index_new + 4].f = data[threei + 1].f;
|
|
args[arg_to_integ] = new_data[index_new + 5].x;
|
|
new_data[index_new + 5].f = function(args);
|
|
|
|
args[arg_to_integ] = new_data[index_new + 6].x;
|
|
new_data[index_new + 6].f = function(args);
|
|
new_data[index_new + 7].f = data[threei + 2].f;
|
|
args[arg_to_integ] = new_data[index_new + 8].x;
|
|
new_data[index_new + 8].f = function(args);
|
|
|
|
new_abs_d2f_dx[index_new/3] = abs(new_data[index_new].dx * (new_data[index_new].f - 2.0 * new_data[index_new + 1].f + new_data[index_new + 2].f));
|
|
new_abs_d2f_dx[index_new/3 + 1] = abs(new_data[index_new].dx * (new_data[index_new + 3].f - 2.0 * new_data[index_new + 4].f + new_data[index_new + 5].f));
|
|
new_abs_d2f_dx[index_new/3 + 2] = abs(new_data[index_new].dx * (new_data[index_new + 6].f - 2.0 * new_data[index_new + 7].f + new_data[index_new + 8].f));
|
|
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3]);
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3 + 1]);
|
|
new_max_abs_d2f_dx = ABACUS::max(new_max_abs_d2f_dx, new_abs_d2f_dx[index_new/3 + 2]);
|
|
|
|
integ_res.integ_est += new_data[index_new].dx * (new_data[index_new].f + new_data[index_new + 1].f + new_data[index_new + 2].f
|
|
+ new_data[index_new + 3].f + new_data[index_new + 4].f + new_data[index_new + 5].f
|
|
+ new_data[index_new + 6].f + new_data[index_new + 7].f + new_data[index_new + 8].f);
|
|
|
|
integ_res.abs_prec += new_abs_d2f_dx[index_new/3] + new_abs_d2f_dx[index_new/3 + 1] + new_abs_d2f_dx[index_new/3 + 2];
|
|
|
|
index_new += 9;
|
|
|
|
} // else
|
|
|
|
} // for (i < nvals/3)
|
|
|
|
integ_res.rel_prec = integ_res.abs_prec/abs(integ_res.integ_est);
|
|
|
|
integ_res.n_vals = index_new;
|
|
|
|
delete[] data;
|
|
data = new_data;
|
|
|
|
max_abs_d2f_dx = new_max_abs_d2f_dx;
|
|
|
|
delete[] abs_d2f_dx;
|
|
abs_d2f_dx = new_abs_d2f_dx;
|
|
|
|
return;
|
|
}
|
|
|
|
// No implementation with I_table yet for complex...
|
|
Integral_result_CX Integrate_optimal (complex<DP> (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax,
|
|
DP req_rel_prec, DP req_abs_prec, int max_nr_pts)
|
|
{
|
|
if (xmax < xmin) ABACUSerror("Use xmax > xmin in Integrate.");
|
|
|
|
Integral_data_CX integ_dat (function, args, arg_to_integ, xmin, xmax);
|
|
|
|
while ((integ_dat.integ_res.rel_prec > req_rel_prec || integ_dat.integ_res.abs_prec > req_abs_prec) && integ_dat.integ_res.n_vals < max_nr_pts) {
|
|
|
|
integ_dat.Improve_estimate (function, args, arg_to_integ, max_nr_pts);
|
|
}
|
|
|
|
return(integ_dat.integ_res);
|
|
}
|
|
|
|
|
|
} // namespace ABACUS
|