Skip to content

Numerical integration

Gaussian Quadrature

beignet.gauss_laguerre_quadrature 

gauss_laguerre_quadrature(degree)
Source code in src/beignet/_gauss_laguerre_quadrature.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def gauss_laguerre_quadrature(degree):
    degree = int(degree)
    if degree <= 0:
        raise ValueError

    c = torch.zeros(degree + 1)
    c[-1] = 1.0

    m = laguerre_polynomial_companion(c)
    x = torch.linalg.eigvalsh(m)

    dy = evaluate_laguerre_polynomial(x, c)
    df = evaluate_laguerre_polynomial(x, differentiate_laguerre_polynomial(c))
    x = x - (dy / df)

    fm = evaluate_laguerre_polynomial(x, c[1:])
    fm = fm / torch.abs(fm).max()
    df = df / torch.abs(df).max()
    w = 1 / (fm * df)

    w = w / torch.sum(w)

    return x, w

beignet.gauss_legendre_quadrature 

gauss_legendre_quadrature(degree)
Source code in src/beignet/_gauss_legendre_quadrature.py 
 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
def gauss_legendre_quadrature(degree):
    degree = int(degree)

    if degree <= 0:
        raise ValueError

    c = torch.zeros(degree + 1)
    c[-1] = 1.0
    m = legendre_polynomial_companion(c)
    x = torch.linalg.eigvalsh(m)

    dy = evaluate_legendre_polynomial(x, c)
    df = evaluate_legendre_polynomial(x, differentiate_legendre_polynomial(c))
    x -= dy / df

    fm = evaluate_legendre_polynomial(x, c[1:])

    fm /= torch.abs(fm).max()
    df /= torch.abs(df).max()

    w = 1 / (fm * df)

    a = torch.flip(w, dims=[0])
    b = torch.flip(x, dims=[0])

    w = (w + a) / 2
    x = (x - b) / 2

    w = w * (2.0 / torch.sum(w))

    return x, w

beignet.gauss_physicists_hermite_polynomial_quadrature 

gauss_physicists_hermite_polynomial_quadrature(degree)
Source code in src/beignet/_gauss_physicists_hermite_polynomial_quadrature.py 
 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
def gauss_physicists_hermite_polynomial_quadrature(degree):
    degree = int(degree)
    if degree <= 0:
        raise ValueError

    c = torch.zeros(degree + 1)
    c[-1] = 1.0

    x = torch.linalg.eigvalsh(physicists_hermite_polynomial_companion(c))

    if degree == 0:
        output = torch.full(x.shape, 1 / math.sqrt(math.sqrt(math.pi)))
    else:
        a1 = torch.zeros_like(x)

        b1 = torch.ones_like(x) / math.sqrt(math.sqrt(math.pi))

        size = torch.tensor(degree)

        for _ in range(0, degree - 1):
            previous = a1

            a1 = -b1 * torch.sqrt((size - 1.0) / size)

            b1 = previous + b1 * x * torch.sqrt(2.0 / size)

            size = size - 1.0

        output = a1 + b1 * x * math.sqrt(2.0)

    dy = output

    n = degree - 1

    if n == 0:
        df = torch.full(x.shape, 1 / math.sqrt(math.sqrt(math.pi)))
    else:
        a = torch.zeros_like(x)

        b = torch.ones_like(x) / math.sqrt(math.sqrt(math.pi))

        size = torch.tensor(n)

        for _ in range(0, n - 1):
            previous = a

            a = -b * torch.sqrt((size - 1.0) / size)

            b = previous + b * x * torch.sqrt(2.0 / size)

            size = size - 1.0

        df = a + b * x * math.sqrt(2.0)

    df = df * math.sqrt(2 * degree)

    x = x - (dy / df)

    n = degree - 1

    if n == 0:
        fm = torch.full(x.shape, 1 / math.sqrt(math.sqrt(math.pi)))
    else:
        a = torch.zeros_like(x)

        b = torch.ones_like(x) / math.sqrt(math.sqrt(math.pi))

        size = torch.tensor(n)

        for _ in range(0, n - 1):
            previous = a

            a = -b * torch.sqrt((size - 1.0) / size)

            b = previous + b * x * torch.sqrt(2.0 / size)

            size = size - 1.0

        fm = a + b * x * math.sqrt(2.0)

    fm = fm / torch.abs(fm).max()

    w = 1 / (fm * fm)

    a = torch.flip(w, dims=[0])
    b = torch.flip(x, dims=[0])

    w = (w + a) / 2
    x = (x - b) / 2

    w = w * (math.sqrt(math.pi) / torch.sum(w))

    return x, w

beignet.gauss_probabilists_hermite_polynomial_quadrature 

gauss_probabilists_hermite_polynomial_quadrature(degree)
Source code in src/beignet/_gauss_probabilists_hermite_polynomial_quadrature.py 
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
def gauss_probabilists_hermite_polynomial_quadrature(degree):
    degree = int(degree)
    if degree <= 0:
        raise ValueError

    c = torch.zeros(degree + 1)
    c[-1] = 1.0
    m = probabilists_hermite_polynomial_companion(c)
    x = torch.linalg.eigvalsh(m)

    if degree == 0:
        dy = torch.full(x.shape, 1.0 / math.sqrt(math.sqrt(2.0 * math.pi)))
    else:
        a1 = torch.zeros_like(x)
        b1 = torch.ones_like(x) / math.sqrt(math.sqrt(2.0 * math.pi))

        size = torch.tensor(degree)

        for _ in range(0, degree - 1):
            previous = a1

            a1 = -b1 * torch.sqrt((size - 1.0) / size)

            b1 = previous + b1 * x * torch.sqrt(1.0 / size)

            size = size - 1.0

        dy = a1 + b1 * x

    n = degree - 1

    if n == 0:
        df = torch.full(x.shape, 1 / math.sqrt(math.sqrt(math.pi)))
    else:
        a = torch.zeros_like(x)

        b = torch.ones_like(x) / math.sqrt(math.sqrt(math.pi))

        size = torch.tensor(n)

        for _ in range(0, n - 1):
            previous = a

            a = -b * torch.sqrt((size - 1.0) / size)

            b = previous + b * x * torch.sqrt(2.0 / size)

            size = size - 1.0

        df = a + b * x * math.sqrt(2.0)

    df = df * math.sqrt(degree)

    x = x - (dy / df)

    n = degree - 1

    if n == 0:
        fm = torch.full(x.shape, 1.0 / math.sqrt(math.sqrt(2.0 * math.pi)))
    else:
        a = torch.zeros_like(x)
        b = torch.ones_like(x) / math.sqrt(math.sqrt(2.0 * math.pi))

        size = torch.tensor(n)

        for _ in range(0, n - 1):
            previous = a

            a = -b * torch.sqrt((size - 1.0) / size)

            b = previous + b * x * torch.sqrt(1.0 / size)

            size = size - 1.0

        fm = a + b * x

    fm = fm / torch.abs(fm).max()

    w = 1 / (fm * fm)

    a = torch.flip(w, dims=[0])
    b = torch.flip(x, dims=[0])

    w = (w + a) / 2
    x = (x - b) / 2

    w = w * (math.sqrt(2 * math.pi) / torch.sum(w))

    return x, w