Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 07:49, 16 December 2015
, 8 years agoAdded Sidef
m (→version 2: added/changed whitespace and comments, made other changes.) |
m (Added Sidef) |
||
Line 1,480:
────── ───────────────────────────────────────────────────────────────────────────────── ─────────────
20.0357498548198037979491872389316561203562082463657269288113065020927852103607 {exact value}
</pre>
=={{header|Sidef}}==
{{trans|Perl 6}}
<lang ruby>func legendre_pair((1), x) { (x, 1) }
func legendre_pair( n, x) {
var (m1, m2) = legendre_pair(n - 1, x);
var u = (1 - 1/n);
((1 + u)*x*m1 - u*m2, m1);
}
func legendre((0), _) { 1 }
func legendre( n, x) { [legendre_pair(n, x)][0] }
func legendre_prime({ .is_zero }, _) { 0 }
func legendre_prime({ .is_one }, _) { 1 }
func legendre_prime(n, x) {
var (m0, m1) = legendre_pair(n, x);
(m1 - x*m0) * n / (1 - x**2);
}
func approximate_legendre_root(n, k) {
# Approximation due to Francesco Tricomi
var t = ((4*k - 1) / (4*n + 2));
(1 - ((n - 1)/(8 * n**3))) * (Math.pi * t -> cos);
}
func newton_raphson(f, f_prime, r, eps = 2e-16) {
while (var dr = (-f(r) / f_prime(r)) -> abs >= eps) {
r += dr;
}
return r;
}
func legendre_root(n, k) {
newton_raphson(legendre.method(:call, n), legendre_prime.method(:call, n),
approximate_legendre_root(n, k));
}
func weight(n, r) { 2 / ((1 - r**2) * legendre_prime(n, r)**2) }
func nodes(n) {
gather {
take(Pair(0, weight(n, 0))) if n.is_odd;
(n/2 -> int).times { |i|
var r = legendre_root(n, i);
var w = weight(n, r);
take(Pair(r, w), Pair(-r, w));
}
}
}
func quadrature(n, f, a, b, nds = nodes(n)) {
func scale(x) { (x*(b - a) + a + b) / 2 }
(b - a) / 2 * nds.map{ .second * f(scale(.first)) }.sum
}
[@(5..10), 20].each { |i|
printf("Gauss-Legendre %2d-point quadrature ∫₋₃⁺³ exp(x) dx ≈ %.15f\n",
i, quadrature(i, {.exp}, -3, +3))
}</lang>
{{out}}
<pre>
Gauss-Legendre 5-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035577718385561
Gauss-Legendre 6-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035746975092344
Gauss-Legendre 7-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749819726600
Gauss-Legendre 8-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749854494515
Gauss-Legendre 9-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749854817432
Gauss-Legendre 10-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749854819791
Gauss-Legendre 20-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.035749854819805
</pre>
| |||