C++ integration

Bonjour,

Ma connaissance insuffisante du C++ me bloque, probablemment a cause des
templates.
Je voulais utiliser une routine que j'ai trouve sur githib, mais sans
commentaire.

Voici
Ceci fonctionne pour integrer integrand2 (une fonction qui n'a pas de parametre),
mais comment puisse integrer func (qui necessite des parameters) ?
Je peux faire
func fcn (20, 4.) ;
mais apres ?

Merci d'avance

#include <iostream>
#include <functional>
#include "clenshaw_curtis_quadrature.h"

double integrand2 (const double x) {
return (x+1)*(x-2) ;
} ;

struct func {
  double p1;
  double p2;
  func (double p1_, double p2_) : p1(p1_), p2(p2_) {}
  double operator()(const double x)
  {
    return x*x*x + p1*x*x + p2*x;
  }
};

int main()
{
    double result = ejmahler_integration::integrateClenshawCurtis<8>(integrand2, -1.0, 2.0) ;
    std::cout << result << std::endl;
}

Disponible sur
https://github.com/ejmahler/constexpr-case-study/tree/master
-------------------------------------
clenshaw_curtis_quadrature.h
------------------------
#include <array>

#ifndef M_PI
#define M_PI           3.14159265358979323846
#endif

#include "discrete_cosine_transform.h"
#include "square_matrix.h"

namespace ejmahler_integration {
    template<size_t N, typename floating_t, class Function>
    floating_t integrateClenshawCurtis(Function integrand, floating_t from, floating_t to);

    template<size_t N, typename floating_t>
    constexpr std::array<floating_t, N/2+1> clenshawCurtisPoints(void);

    template<size_t N, typename floating_t>
    constexpr std::array<floating_t, N/2+1> clenshawCurtisWeights(void);
}

template<size_t N, typename floating_t, class Function>
floating_t ejmahler_integration::integrateClenshawCurtis(Function integrand, floating_t a, floating_t b)
{
    constexpr std::array<floating_t, N/2+1> quadraturePoints = clenshawCurtisPoints<N, floating_t>();
    constexpr std::array<floating_t, N/2+1> quadratureWeights = clenshawCurtisWeights<N, floating_t>();

    floating_t halfDiff = (b - a) / 2;
    floating_t halfSum = (a + b) / 2;

    floating_t sum(0);
    for(size_t i = 0; i < N/2; i++)
    {
        sum += quadratureWeights[i] * (integrand(halfSum + halfDiff * quadraturePoints[i]) + integrand(halfSum - halfDiff * quadraturePoints[i]));
    }

    //the final element in quadraturePoints will always be zero. so if we compute it inside the loop it'll be double-evaluated
    sum += quadratureWeights[N/2] * integrand(halfSum);
    return halfDiff * sum;
}

template<size_t N, typename floating_t>
constexpr std::array<floating_t, N/2+1> ejmahler_integration::clenshawCurtisPoints(void)
{
    std::array<floating_t, N/2+1> result{};

    for(size_t i = 0; i <= N/2; i++)
    {
        result[i] = std::cos(i * floating_t(M_PI) / N);
    }

    return result;
}

template<size_t N, typename floating_t>
constexpr std::array<floating_t, N/2+1> ejmahler_integration::clenshawCurtisWeights(void)
{
    float scale = floating_t(2) / N;

    std::array<floating_t, N/2+1> input{};
    for(size_t i = 0; i < N/2+1; i++)
    {
        input[i] = scale / (1 - floating_t(4 * i * i));
    }

    auto result = ejmahler_integration::dctType1(input);

    result[0] *= .5;

    return result;
}
------------------------------------
===========================================================================
 Patrick DUPRÉ                                 | | email: pdupre@???
===========================================================================