Anonymous user
Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 14:26, 10 July 2013
, 10 years agoFortran added
(Fortran added) |
|||
Line 305:
Integrating exp(x) over [-3, 3]: 20.035577718386
Compred to actual: 20.035749854820</pre>
=={{header|Fortran}}==
<lang 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, z1, z2
do n = 1,20
a = -3
b = 3
allocate (r(2,n))
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
deallocate(r)
end do
contains
function polyval(coeff, x) result (y)
real(kind=p) :: coeff(:), x, y
integer :: k
y = coeff(1)
do k = 2, size(coeff)
y = y*x+coeff(k)
end do
end function
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(:)
allocate(p0(1), p1(2))
p0 = [1._p]
p1 = [1._p, 0._p]
do k = 2, n
allocate(tmp(k+1))
tmp = ((2*k-1)*[p1,0._p]-(k-1)*[0._p, 0._p,p0])/k
deallocate(p0)
call move_alloc(p1, p0)
call move_alloc(tmp, p1)
end do
do i = 1, n
x = cos(pi*(i-0.25_p)/(n+0.5_p))
iter = 0
do iter = 1, 10
f = polyval(p1, x)
df = n/(x*x-1)*(x*f-polyval(p0,x))
dx = f / df
x = x - dx
if (abs(dx)<10*epsilon(dx)) exit
end do
r(1,i) = x
r(2,i) = 2/((1-x**2)*df**2)
end do
end function
end program
</lang>
<pre>
n numerical integral error
--------------------------------------------------
1 6.00000000000000000000000000000000000 -14.
2 17.4874646410555689643606840462449514 -2.5
3 19.8536919968055821921309108927158374 -0.18
4 20.0286883952907008527738054439857600 -0.71E-02
5 20.0355777183855621539285357252750835 -0.17E-03
6 20.0357469750923438830654575585499463 -0.29E-05
7 20.0357498197266007755718729372892887 -0.35E-07
8 20.0357498544945172882260918041685114 -0.33E-09
9 20.0357498548174338368864419454853160 -0.24E-11
10 20.0357498548197898711175766908545118 -0.14E-13
11 20.0357498548198037305529147159728024 -0.67E-16
12 20.0357498548198037976759531014515043 -0.27E-18
13 20.0357498548198037979482458118938813 -0.94E-21
14 20.0357498548198037979491844483850590 -0.28E-23
15 20.0357498548198037979491872316559285 -0.73E-26
16 20.0357498548198037979491872389497857 0.18E-28
17 20.0357498548198037979491872391086672 0.18E-27
18 20.0357498548198037979491872397203703 0.79E-27
19 20.0357498548198037979491872404290602 0.15E-26
20 20.0357498548198037979491872397574252 0.83E-27
</pre>
=={{header|Go}}==
| |||