Rosetta Code

Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

Page Discussion

Numerical integration/Gauss-Legendre Quadrature (view source)

Revision as of 22:01, 27 January 2012

238 bytes added ,  12 years ago
Updated D version
(D code: the task asked for N=5)
(Updated D version)
Line 231:
<lang d>import std.stdio, std.math;
 
final /*const*/ final class GaussLegendreQuadrature(size_t N, FP=double,
size_t NBITS=50) {
/*const*/ public double[N] lroots, weight;
//FP[N + 1][N + 1] lcoef = 0.0;
FP[N + 1][N + 1] lcoef;
 
this() pure nothrow {
foreachFP[N (ref+ row;1][N + 1] lcoef) = 0.0;
legendreCoef(/*ref*/ row[] = 0.0lcoef);
legendreCoeflegendreRoots(lcoef);
legendreRoots();
}
 
static private void legendreCoef()ref pureFP[N nothrow+ {1][N + 1] lcoef)
pure nothrow {
lcoef[0][0] = lcoef[1][1] = 1;
foreach (int n; 2 .. N + 1) { // n must be signed
Line 254 ⟶ 252:
}
 
private FPvoid legendreEvallegendreRoots(inconst intref n,FP[N in+ FP1][N + 1] xlcoef)
const pure nothrow {
static FP slegendreEval(const =ref lcoefFP[nN + 1][nN + 1]; lcoef,
in int n, in FP x) pure nothrow {
foreach_reverse (i; 1 .. n+1)
FP s = s * x + lcoef[n][i - 1n];
return s foreach_reverse (i; 1 .. n+1)
s = s * x + lcoef[n][i - 1];
}▼
return s;
▲ }
 
private static FP legendreDiff(inconst intref n,FP[N in+ FP1][N x)+ 1] lcoef,
const in int n, in FP x) pure nothrow {
return n * (x * legendreEval(lcoef, n, x) - legendreEval(n - 1, x))/
(x ^^ 2 legendreEval(lcoef, n - 1, x);) /
(x ^^ 2 - 1);
}
 
private void legendreRoots() pure nothrow {
foreach (i; 1 .. N + 1) {
FP x = cos(PI * (i - 0.25) / (N + 0.5));
Line 274:
do {
x1 = x;
x -= legendreEval(lcoef, N, x) / legendreDiff(N, x);
legendreDiff(lcoef, N, x);
} while (feqrel(x, x1) < NBITS);
lroots[i - 1] = x;
 
x1 = legendreDiff(lcoef, N, x);
weight[i - 1] = 2 / ((1 - x ^^ 2) * (x1 ^^ 2));
}
}
 
public FP integrate(in FP function(FP x) f,
in FP a, in FP b) const {
immutable FP c1 = (b - a) / 2,
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/183689"