Anonymous user
Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 08:26, 11 July 2013
, 10 years ago→{{header|Fortran}}: use of automatic re-allocation, more accurate evaluation of P' using Horner's method
m (→{{header|Fortran}}: use of automatic re-allocation, more accurate evaluation of P' using Horner's method) |
|||
Line 307:
=={{header|Fortran}}==
<lang Fortran>
! -assume realloc_lhs
! when compiled with Intel Fortran.
program gauss
implicit none
integer, parameter :: p = 16 ! quadruple precision
integer :: n = 10, k
real(kind=p), allocatable :: r(:,:)
real(kind=p) :: z, a, b, exact
do n = 1,20
a = -3; b = 3
b = 3▼
r = gaussquad(n)
z = (b-a)/2*dot_product(r(2,:),exp((a+b)/2+r(1,:)*(b-a)/2))
exact = exp(3.0_p)-exp(-3.0_p)
print "(i0,1x,g0,1x,g10.2)",n, z, z-exact
end do
contains
function gaussquad(n) result(r)
integer :: n
real(kind=p), parameter :: pi = 4*atan(1._p)
real(kind=p) :: r(2, n), x, f, df, dx
integer :: i, iter
real(kind = p), allocatable :: p0(:), p1(:), tmp(:)
p0 = [1._p]
Line 348 ⟶ 338:
do k = 2, n
tmp = ((2*k-1)*[p1,0._p]-(k-1)*[0._p, 0._p,p0])/k
end do
do i = 1, n
Line 358 ⟶ 345:
iter = 0
do iter = 1, 10
f =
df = f + x*df
f = p1(k) + x * f
dx = f / df
x = x - dx
Line 374 ⟶ 364:
n numerical integral error
--------------------------------------------------
1 6.
2 17.
3 19.
4 20.
5 20.
6 20.
7 20.
8 20.
9 20.
10 20.
11 20.
12 20.
13 20.
14 20.
15 20.
16 20.
17 20.
18 20.
19 20.
20 20.
</pre>
| |||