Adaptive Simpson's Method (Adaptive Integral)

출처: https://discuss.codechef.com/t/a-tutorial-on-adaptive-simpsons-method/19991

A tutorial on Adaptive Simpson's Method

위 출처를 보면 쉽게 adaptive integral 구현이 가능함.

아래의 예제는 f1 이라는 복소수 값을 가지는 함수를 0.0~tf 구간에서 32 포인트에 대한 uniform 적분과 Adaptive 적분 결과와 비교

> uniform이 확실히 빠르긴 하지만 계산의 정확도에 대한 판단이 어려움

> f1 함수가 변수 n과 z를 입력받아야 해서 calc 와 simpson에도 변수 n과 z를 같이 넣어줌

#include <iostream>
#include <cmath>
#include <fstream>
#include <string>
#include <iomanip>
#include <time.h>
#include <ostream>
#include <complex>
#include <list>
 
using namespace std;
 
#define pi 3.1415926535897932384626433832795028841971693993751058209749445923078164
const double eps = 1e-10; // represents how much precision you need
 
complex<double> sinhs(complex<double> x)
{
    return (abs(x) <= 1E-10) ? 1. : sinh(x) / x;
}
 
complex<double> f1(double n, complex<double> z, double t)
{
    if (fabs(t) <= 1E-4)
    {
        complex<double> s = sinhs(0.5 * z * t);
        return (z * z * s * s) / (2 * n * pi * cosh(0.5 * pi * t) * sinhs(n * pi * t));
    }    
    else
    {
        double ep = t * (1 - exp(-2 * n * pi * t));
        complex<double> ez = exp(-z * t);
        complex<double> temp1 = 2. * exp((z - 0.5 * pi - n * pi) * t) * (1. + ez * ez - 2. * ez);
        complex<double> temp2 = ((1. + exp(-pi * t)) * ep);
        return temp1 / temp2;
    }
}
 
complex<double> simpson(double l, double r, double n, complex<double> z)
{
    return (r - l) / 6.0 * (f1(n, z, l) + 4.0 * f1(n, z, (l + r) / 2) + f1(n, z, r));
}
complex<double> calc(double l, double r, complex<double> ans, double n, complex<double> z)
{
    double m = (l + r) / 2;
    complex<double> x = simpson(l, m, n, z);
    complex<double> y = simpson(m, r, n, z);
    if (abs(x + y - ans) < eps)
        return x + y;
    return calc(l, m, x, n, z) + calc(m, r, y, n, z);
}
 
int main()
{
    complex<double> jj = complex<double>(0.0, 1.0);
    complex<double> z = complex<double>(-4.5728, 1.3170);
    double n = 1.0;
 
    // move to the right half plane because psi(-z) = psi(z)
    if (z.real() < 0)
        z = -z;
    // move to the analytic strip z0 - 4Phi < real(z) < z0 < -- - (2Phi = n * pi)
    double z0 = 0.5 * pi + n * pi;
    int m = (int)ceil((real(z) - z0) / (2.0 * n * pi));
    z = z - 2.0 * n * pi * m;
 
    
    // original form
    double tf = 20. / (z0 - real(z));
    tf = (50 < tf) ? 50 : tf; // 0~0.5 정도로 적분해도 충분함 경험상 100 넘어가면 적분이 발산함
 
    // uniform integral
    int kk, subInterval;
    complex<double> integration, result, integration_cf;
    double lower, upper, stepSize, ll;
    lower = 0.0;
    upper = tf;
    subInterval = 32;
 
    double start, end;
    start = clock();
    integration = f1(n, z, lower) + f1(n, z, upper);
    stepSize = (upper - lower) / subInterval;
    
    for (kk = 1; kk <= subInterval - 1; kk++)
    {
        ll = lower + kk * stepSize;
        integration = integration + 2.0 * f1(n, z, ll);
    }
    integration = integration * stepSize/2.0;
 
    end = clock();
    cout << end - start << endl;
 
    // Adpative Simpson's Method
    start = clock();
    integration_cf = calc(lower, upper, simpson(lower, upper, n, z), n, z);
    end = clock();
    cout << end - start << endl;
 
    result = exp(integration * -0.5);
 
    // back to original z
    for (int k = 1; k <= m; k++)
    {
        result /= tan(0.5 * z + n * pi / 2 + 0.25 * pi);
        z += 2.0 * n * pi;
    }
 
    cout << "Result:" << result << endl;
    return 0;
}

---

Adapative Quadrature 예제

> 참고: https://stackoverflow.com/questions/12906618/adaptive-quadrature-c

Adaptive Quadrature (C++)

아래는 체크해보지는 않았지만... 조금 손보면 사용할 수 있을 듯

ouble trap_rule(double a, double b, double (*f)(double),double tolerance, int count)
{
    double coarse = coarse_helper(a,b, f); //getting the coarse and fine approximations from the helper functions
    double fine = fine_helper(a,b,f);
 
    if ((abs(coarse - fine) <= 3.0*tolerance) && (count >= minLevel))
        //return fine if |c-f| <3*tol, (fine is "good") and if count above
        //required minimum level
    {
        return fine;
    }
    else if (count >= maxLevel)
        //maxLevel is the maximum number of recursion we can go through
    {
        return fine;
 
    }
    else
    {
        //if none of these conditions are satisfied, split [a,b] into [a,c] and [c,b] performing trap_rule
        //on these intervals -- using recursion to adaptively approach a tolerable |coarse-fine| level
        //here, (a+b)/2 = c
        ++count;
        return  (trap_rule(a, (a+b)/2.0, f, tolerance/2.0, count) + trap_rule((a+b)/2.0, b, f, tolerance/2.0, count));
    }
}
 
 
EDIT: Helper and test functions:
 
//test function
double function_1(double a)
{
    return pow(a,2);
}
 
//"true" integral for comparison and tolerances
 
 
 
 
//helper functions
 
double coarse_helper(double a, double b, double (*f)(double))
{
    return 0.5*(b - a)*(f(a) + f(b)); //by definition of coarse approx
}
 
 
double fine_helper(double a, double b, double (*f)(double))
{
    double c = (a+b)/2.0;
    return 0.25*(b - a)*(f(a) + 2*f(c) + f(b)); //by definition of fine approx
}
 
double helper(double a, double b, double (*f)(double x), double tol)
{
    return trap_rule(a, b, f, tol, 1);
}
 
void main(void)
{
    std::cout << "First we approximate the integral of f(x) = x^2 on [0,2]" << std::endl;
 
    std::cout << "Enter a: ";
    std::cin >> a;
 
    std::cout << "Enter b: ";
    std::cin >> b;
 
    true_value1 = analytic_first(a,b);
 
    for (int i = 0; i<=8; i++)
    {
        result1 [i] = helper(a, b, function_1, tolerance[i]);
        error1 [i] = fabs(true_value1 - result1 [i]);
    }
 
    std::cout << "(Approximate integral of x^2, tolerance, error )" << std::endl;
    for (int i = 0; i<=8; i++)
    {
        std::cout << "(" << result1 [i] << "," << tolerance[i] << "," << error1[i] << ")" << std::endl;
    }
 
}