-
Notifications
You must be signed in to change notification settings - Fork 0
/
quadrature.F90
189 lines (149 loc) · 4.77 KB
/
quadrature.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
module m_quadrature
implicit none
double precision, parameter :: PI = acos(-1.0d0)
public
contains
! Arfken p. 541
! Legendre polynomials
double precision pure recursive function Legendre(n,x) result(returnval)
integer, intent(in) :: n
double precision, intent(in) :: x
double precision :: Pn1, Pn2
if (n == 0) then
returnval = 1
else if (n == 1) then
returnval = x
else
Pn2 = Legendre(n-2,x)
Pn1 = Legendre(n-1,x)
returnval = 2*x*Pn1 - Pn2 - (x*Pn1-Pn2)/n
end if
end function
! Arfken p. 541
! Derivatives of Legendre polynomials
double precision pure recursive function DLegendre(n,x) result(returnval)
integer, intent(in) :: n
double precision, intent(in) :: x
double precision :: DPn1, Pn1
if (n == 0) then
returnval = 0
else if (n == 1) then
returnval = 1
else
DPn1 = DLegendre(n-1,x)
Pn1 = Legendre(n-1, x)
returnval = n * Pn1 + x * DPn1
end if
end function
! Second derivative of Legendre polynomials
double precision pure recursive function DDLegendre(n,x) result(returnval)
integer, intent(in) :: n
double precision, intent(in) :: x
double precision :: DPn1, DDPn1
if (n == 0) then
returnval = 0
else if (n == 1) then
returnval = 0
else
DPn1 = DLegendre(n-1,x)
DDPn1 = DDLegendre(n-1,x)
returnval = (n + 1.0) * DPn1 + x* DDPn1
end if
end function
! Gauss-Legendre quadrature points (Roots of Legendre polynomial of order N)
subroutine GL_Points(pts)
double precision, intent(inout) :: pts(:)
double precision, parameter :: accuracy = 1.0d-15
integer :: N, i, j
double precision :: a(size(pts)), linsp_start, linsp_end, magdiff, tmp
N = size(pts)
! Initial guess of node positions
! Newman p. 523
linsp_start = 3.0
linsp_end = 4.0*N-1
a = [(linsp_start + (linsp_end-linsp_start) * (i-1) / (N-1), i=1, N)]
a = a / (4*N+2)
pts = cos(PI*a+1/(8*N*N*tan(a)))
! Newton's method root finder
magdiff = 2.0 * accuracy
do while (magdiff > accuracy)
a = pts - [( Legendre(N, pts(i)), i=1, N )] / [( DLegendre(N, pts(j)), j=1, N )]
magdiff = maxval(abs(a - pts))
pts = a
end do
! Reverse the ordering so the negative ones come first
do i = 1, N/2
tmp = pts(N+1-i)
pts(N+1-i) = pts(i)
pts(i) = tmp
end do
end subroutine GL_Points
! Gauss-Legendre-Lobatto quadrature points (Roots of Legendre polynomial derivative of order N-1,
! plus 1 and -1).
subroutine GLL_Points(pts)
double precision, intent(inout) :: pts(:)
double precision, parameter :: accuracy = 1.0d-15
integer :: N, i, j
double precision :: a(size(pts)-2) , magdiff, tmp
double precision :: interiorPts(size(pts)-2)
N = size(pts)
call GL_Points(interiorPts)
! Newton's method root finder
magdiff = accuracy * 2.0
do while (magdiff > accuracy)
a = interiorPts - [( DLegendre(N-1, interiorPts(i)), i=1, N-2 )] / [( DDLegendre(N-1, interiorPts(j)), j=1, N-2 )]
magdiff = maxval(abs(a - interiorPts ))
interiorPts = a
end do
pts(N) = 1.0
pts(2:N-1) = interiorPts
pts(1) = -1.0
end subroutine GLL_Points
subroutine GL_Weights(pts, wghts)
double precision, intent(in) :: pts(:)
double precision, intent(inout) :: wghts(:)
double precision :: dleg(size(pts))
integer :: N , i
N = size(pts)
dleg = [( DLegendre(N, pts(i)), i = 1, N )]
wghts = 2.0 / ((1-pts**2) * (dleg**2))
end subroutine
subroutine GLL_Weights(pts, wghts)
double precision, intent(in) :: pts(:)
double precision, intent(inout) :: wghts(:)
double precision :: leg(size(pts))
integer :: N , i
N = size(pts)
leg = [( Legendre(N-1, pts(i)), i = 1, N )]
wghts = 2.0 / ( N*(N-1) * leg**2 )
end subroutine
! value of the i-th Lagrange polynomial on the [-1,1] domain at point x
! gx are the quadrature points
double precision pure function Lagrange(gx, i, x)
double precision, intent(in) :: gx(:), x
integer, intent(in) :: i
double precision :: prod
integer :: j
prod = 1.0
do j = 1, size(gx)
if(j /= i) prod = prod * (x - gx(j))/(gx(i)-gx(j))
end do
Lagrange = prod
end function
! value of the derivative of the i-th Lagrange polynomial on the [-1,1] domain at point x
! gx are the quadrature points
! Calculation takes advantage of log differentiation, f'/f = log(f)'
double precision pure function DLagrange(gx,i,x)
integer,intent(in) :: i
double precision,intent(in) :: gx(:), x
double precision :: L
integer :: m
L = 0.0
do m = 1,size(gx)
if (m /= i) then
L = L + 1.0/(x-gx(m))
end if
end do
DLagrange = L * Lagrange(gx,i,x)
end function
end module