Coverage for src/methodsnm/intrule_1d.py: 65%
78 statements
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-27 13:22 +0000
1"""
2This module provides classes for 1D numerical integration rules.
3"""
5from abc import ABC, abstractmethod
6import numpy as np
7from numpy import array
8from methodsnm.intrule import IntRule
10class IntRule1D(IntRule):
11 """
12 Abstract base class for 1D numerical integration rules.
13 """
14 interval = None
15 def __init__(self, interval=(0,1)):
16 """
17 Initializes the integration rule with the given interval.
18 """
19 self.interval = interval
22def evaluate_exactness_degree(rule, max_order=None):
23 """
24 Evaluates the exactness degree of the given integration rule.
26 Parameters:
27 rule (IntRule1D): The integration rule to evaluate.
28 max_order (int or None): maximum order to check for exactness degree. If None, the exactness degree is checked until it is found.
30 Returns:
31 int: The exactness degree of the integration rule.
32 """
33 a,b = rule.interval
34 i=0
35 while True:
36 if not np.isclose(rule.integrate(lambda x: x**i), (b**(i+1)-a**(i+1))/(i+1), rtol=1e-12):
37 if i == 0:
38 print("Warning: Exactness degree is below 0")
39 return i-1
40 i += 1
41 if max_order is not None and i > max_order:
42 raise ValueError("Could not determine exactness degree")
44class MidPointRule(IntRule1D):
45 """
46 Class for the midpoint rule for 1D numerical integration.
47 """
48 def __init__(self, interval = (0,1)):
49 """
50 Initializes the midpoint rule with the given interval.
52 Parameters:
53 interval (tuple): The interval to integrate over.
54 """
55 self.interval = interval
56 a,b = interval
57 self.nodes = array([[0.5*a+0.5*b]])
58 self.weights = array([1.0*(b-a)])
59 self.exactness_degree = 1
62class NewtonCotesRule(IntRule1D):
63 """
64 Class for the Newton-Cotes rule for 1D numerical integration.
65 """
66 def __init__(self, n =None, nodes=None, interval=(0,1)):
67 """
68 Initializes the Newton-Cotes rule with the given interval and number of nodes.
70 Parameters:
71 n (int or list): The number of nodes or a list of nodes to use for integration.
72 nodes (list): A list of nodes to use for integration.
73 interval (tuple): The interval to integrate over.
74 """
75 self.interval = interval
76 a,b = interval
77 if nodes is None and n is None:
78 raise ValueError("Either n or nodes must be specified")
79 if isinstance(n, list):
80 nodes = n
81 else:
82 nodes = np.linspace(a,b,n)
83 n = len(nodes)
85 self.nodes = array(nodes)
87 h = b-a
88 # Compute weights on [0,1]
89 x = [(nodes[i]-a)/h for i in range(n)]
90 # sum w[i] * x[i]**k = 1/(k+1)
91 A = np.array([[x[i]**k for i in range(n)] for k in range(n)])
92 b = np.array([1.0/(k+1) for k in range(n)])
93 w = np.linalg.solve(A,b)
94 self.weights = w*h
95 self.exactness_degree = evaluate_exactness_degree(self)
97class NP_GaussLegendreRule(IntRule1D):
98 """
99 Wrapper class for the Gauss-Legendre rule for 1D numerical integration of numpy.
100 """
101 def __init__(self, n, interval=(0,1)):
102 """
103 Initializes the Gauss-Legendre rule with the given interval and number of nodes.
105 Parameters:
106 n (int): The number of nodes to use for integration.
107 interval (tuple): The interval to integrate over.
108 """
109 self.interval = interval
110 a,b = interval
111 self.nodes = np.empty((n,1))
112 nodes = np.polynomial.legendre.leggauss(n)[0]
113 self.nodes[:,0] = 0.5*(a+b) + 0.5*(b-a)*nodes
114 self.weights = np.polynomial.legendre.leggauss(n)[1]
115 self.weights = 0.5*(b-a)*self.weights
116 self.exactness_degree = 2*n-1
118class GaussLegendreRule(IntRule1D):
119 """
120 Class for the Gauss-Legendre rule for 1D numerical integration.
121 """
122 def __init__(self, n, interval=(0,1)):
123 """
124 Initializes the Gauss-Legendre rule with the given interval and number of nodes.
126 Parameters:
127 n (int): The number of nodes to use for integration.
128 interval (tuple): The interval to integrate over.
129 """
130 self.interval = interval
131 a,b = interval
133 raise NotImplementedError("Not implemented")
135import scipy
136class SP_GaussJacobiRule(IntRule1D):
137 """
138 Wrapper class for the Gauss-Jacobi rule for 1D numerical integration of numpy.
139 """
140 def __init__(self, n, alpha, beta, interval=(0,1)):
141 """
142 Initializes the Gauss-Legendre rule with the given interval and number of nodes.
144 Parameters:
145 n (int): The number of nodes to use for integration.
146 interval (tuple): The interval to integrate over.
147 alpha, beta (float): The parameters of the Jacobi polynomial.
148 """
149 self.interval = interval
150 a,b = interval
151 self.alpha = alpha
152 self.beta = beta
153 nodes, weights = scipy.special.roots_jacobi(n,alpha,beta)
154 self.nodes = np.empty((n,1))
155 self.nodes[:,0] = 0.5*(a+b) + 0.5*(b-a)*nodes
156 self.weights = 0.5**(alpha+beta+1)*(b-a)*weights
157 self.exactness_degree = None
159class GaussJacobiRule(IntRule1D):
160 """
161 Class for the Gauss-Jacobi rule for 1D numerical integration.
162 """
163 def __init__(self, n, alpha, beta, interval=(0,1)):
164 """
165 Initializes the Gauss-Jacobi rule with the given interval and number of nodes.
167 Parameters:
168 n (int): The number of nodes to use for integration.
169 interval (tuple): The interval to integrate over.
170 alpha, beta (float): The parameters of the Jacobi polynomial.
171 """
172 self.interval = interval
173 a,b = interval
174 self.alpha = alpha
175 self.beta = beta
177 raise NotImplementedError("Not implemented")