Commit ef2f3139 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Implement a new expression for Lagrange polynomials in order to reduce the risk of roundoff errors.

parent 6f4f0e0a
......@@ -21,6 +21,7 @@
#include <cassert>
#include "../../Utils/FNoCopyable.hpp"
#include "../../Utils/FMath.hpp"
/**
......@@ -56,19 +57,36 @@ struct FUnifRoots : FNoCopyable
*/
static FReal L(const unsigned int n, FReal x)
{
//std::cout << x << std::endl;
assert(std::fabs(x)-1.<10.*std::numeric_limits<FReal>::epsilon());
if (std::fabs(x)>1.) {
//std::cout << "x=" << x << " out of bounds!" << std::endl;
x = (x > FReal( 1.) ? FReal( 1.) : x);
x = (x < FReal(-1.) ? FReal(-1.) : x);
}
// Specific precomputation of scale factor
// in order to minimize round-off errors
// NB: scale factor could be hardcoded (just as the roots)
FReal scale;
int omn = order-n-1;
if(omn%2) scale=-1.; // (-1)^(n-1-(k+1)+1)=(-1)^(omn-1)
else scale=1.;
scale*=FMath::pow(2.,order-1)*FMath::factorial(n)*FMath::factorial(omn);
// compute L
FReal L=FReal(1.);
for(unsigned int m=0;m<order;++m){
if(m!=n)
L *= (x-FUnifRoots<order>::roots[m])/(FUnifRoots<order>::roots[n]-FUnifRoots<order>::roots[m]);
if(m!=n){
// previous version with risks of round-off error
//L *= (x-FUnifRoots<order>::roots[m])/(FUnifRoots<order>::roots[n]-FUnifRoots<order>::roots[m]);
// new version (reducing round-off)
// regular grid on [-1,1] (h simplifies, only the size of the domain and a remains i.e. 2. and -1.)
L *= ((order-1)*(x+1.)-2.*m);
}
}
L*=scale;
return FReal(L);
}
......
......@@ -103,6 +103,16 @@ struct FMath{
return (1 << power);
}
/** To get factorial */
static double factorial(int inValue){
if(inValue==0) return 1.;
else {
double result = inValue;
while(--inValue > 1) result *= inValue;
return result;
}
}
/** To know if a value is between two others */
template <class NumType>
static bool Between(const NumType inValue, const NumType inMin, const NumType inMax){
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment