Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 17:07, 23 August 2018
, 5 years ago→{{header|Sidef}}: updated code
(Adjust epsilon) |
m (→{{header|Sidef}}: updated code) |
||
Line 2,010:
<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)
}
Line 2,022:
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))) * cos(Num.pi * t
}
func newton_raphson(f, f_prime, r, eps = 2e-16) {
loop {
dr.abs >= eps || break
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))
}
Line 2,048 ⟶ 2,050:
func nodes(n) {
gather {
take(Pair(0, weight(n, 0))) if n.is_odd
var r = legendre_root(n, i)
var w = weight(n, r)
take(Pair(r, w), Pair(-r, w))
}.each(1 .. (n >> 1))
}
}
Line 2,059 ⟶ 2,061:
func quadrature(n, f, a, b, nds = nodes(n)) {
func scale(x) { (x*(b - a) + a + b) / 2 }
(b - a) / 2 * nds.
}
| |||