forked from hbayraktaroglu/quadpy
-
Notifications
You must be signed in to change notification settings - Fork 0
/
_gauss_radau.py
41 lines (31 loc) · 1.16 KB
/
_gauss_radau.py
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
import numpy as np
import orthopy
from ..tools import scheme_from_rc
from ._helpers import C1Scheme
def gauss_radau(n, a=0.0, b=0.0):
assert n >= 2
degree = 2 * n - 1
rc = orthopy.c1.jacobi.RecurrenceCoefficients("monic", a, b, symbolic=False)
_, alpha, beta = np.array([rc[k] for k in range(n)]).T
points, weights = _radau(alpha, beta, rc.int_1, -1.0)
return C1Scheme("Gauss-Radau", degree, weights, points)
def _radau(alpha, beta, int_1, xr):
"""From <http://www.scientificpython.net/pyblog/radau-quadrature>:
Compute the Radau nodes and weights with the preassigned node xr.
Based on the section 7 of the paper
Some modified matrix eigenvalue problems,
Gene Golub,
SIAM Review Vol 15, No. 2, April 1973, pp.318--334.
"""
from scipy.linalg import solve_banded
beta[0] = int_1
n = len(alpha) - 1
f = np.zeros(n)
f[-1] = beta[-1]
A = np.vstack((np.sqrt(beta), alpha - xr))
J = np.vstack((A[:, 0:-1], A[0, 1:]))
delta = solve_banded((1, 1), J, f)
alphar = alpha.copy()
alphar[-1] = xr + delta[-1]
x, w = scheme_from_rc(alphar, beta, int_1, mode="numpy")
return x, w