Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

Content added Content deleted
(Added Perl example)
(→‎{{header|Julia}}: function eig was depreciated and replaced with eigen, the dot operator handling has also changed since this example was put up, added the needed module import, removed the return keyword as it is unnecessary in julia)
Line 928: Line 928:
=={{header|Julia}}==
=={{header|Julia}}==
This function computes the points and weights of an ''N''-point Gauss–Legendre quadrature rule on the interval (''a'',''b''). It uses the O(''N''<sup>2</sup>) algorithm described in Trefethen & Bau, ''Numerical Linear Algebra'', which finds the points and weights by computing the eigenvalues and eigenvectors of a real-symmetric tridiagonal matrix:
This function computes the points and weights of an ''N''-point Gauss–Legendre quadrature rule on the interval (''a'',''b''). It uses the O(''N''<sup>2</sup>) algorithm described in Trefethen & Bau, ''Numerical Linear Algebra'', which finds the points and weights by computing the eigenvalues and eigenvectors of a real-symmetric tridiagonal matrix:
<lang julia>function gauss(a,b, N)
<lang julia>using LinearAlgebra

λ, Q = eig(SymTridiagonal(zeros(N), [ n / sqrt(4n^2 - 1) for n = 1:N-1 ]))
function gauss(a, b, N)
return (λ + 1) * (b-a)/2 + a, [ 2*Q[1,i]^2 for i = 1:N ] * (b-a)/2
λ, Q = eigen(SymTridiagonal(zeros(N), [n / sqrt(4n^2 - 1) for n = 1:N-1]))
@. (λ + 1) * (b - a) / 2 + a, [2Q[1, i]^2 for i = 1:N] * (b - a) / 2
end</lang>
end</lang>
(This code is a simplified version of the <code>Base.gauss</code> subroutine in the Julia standard library.)
(This code is a simplified version of the <code>Base.gauss</code> subroutine in the Julia standard library.)
{{out}}
{{out}}
<pre>
<pre>
julia> x,w = gauss(-3,3,5)
julia> x, w = gauss(-3, 3, 5)
([-2.71854, -1.61541, 1.33227e-15, 1.61541, 2.71854], [0.710781, 1.43589, 1.70667, 1.43589, 0.710781])
([-2.718539537815986,-1.6154079303170459,1.3322676295501878e-15,1.6154079303170512,2.7185395378159916],[0.7107806551685667,1.4358860114981036,1.706666666666666,1.435886011498098,0.7107806551685655])


julia> sum(exp(x) .* w)
julia> sum(exp.(x) .* w)
20.03557771838554
20.03557771838554
</pre>
</pre>