diff --git a/pySDC/helpers/spectral_helper.py b/pySDC/helpers/spectral_helper.py index 25777ec31e..469023c296 100644 --- a/pySDC/helpers/spectral_helper.py +++ b/pySDC/helpers/spectral_helper.py @@ -161,6 +161,10 @@ def get_differentiation_matrix(self): def get_integration_matrix(self): raise NotImplementedError() + def get_integration_weights(self): + """Weights for integration across entire domain""" + raise NotImplementedError() + def get_wavenumbers(self): """ Get the grid in spectral space @@ -379,6 +383,16 @@ def get_integration_matrix(self, lbnd=0): raise NotImplementedError(f'This function allows to integrate only from x=0, you attempted from x={lbnd}.') return S + def get_integration_weights(self): + """Weights for integration across entire domain""" + n = self.xp.arange(self.N, dtype=float) + + weights = (-1) ** n + 1 + weights[2:] /= 1 - (n**2)[2:] + + weights /= 2 / self.L + return weights + def get_differentiation_matrix(self, p=1): ''' Keep in mind that the T2T differentiation matrix is dense. @@ -808,6 +822,12 @@ def get_integration_matrix(self, p=1): k[0] = 1j * self.L return self.linalg.matrix_power(self.sparse_lib.diags(1 / (1j * k)), p) + def get_integration_weights(self): + """Weights for integration across entire domain""" + weights = self.xp.zeros(self.N) + weights[0] = self.L / self.N + return weights + def get_plan(self, u, forward, *args, **kwargs): if self.fft_lib.__name__ == 'mpi4py_fft.fftw': if 'axes' in kwargs.keys(): diff --git a/pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py b/pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py index 47f572c88d..ca286fe915 100644 --- a/pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py +++ b/pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py @@ -139,6 +139,35 @@ def test_integration_matrix(N): assert np.allclose(exact.coef[:-1], du) +@pytest.mark.base +@pytest.mark.parametrize('x0', [-1, 0]) +@pytest.mark.parametrize('x1', [0.789, 1]) +@pytest.mark.parametrize('N', [4, 7]) +def test_integral_whole_interval(x0, x1, N): + import numpy as np + from pySDC.helpers.spectral_helper import ChebychevHelper + from qmat.lagrange import LagrangeApproximation + + cheby = ChebychevHelper(N, x0=x0, x1=x1) + x = cheby.get_1dgrid() + + coeffs = np.random.random(N) + coeffs[-1] = 0 + + u_hat = coeffs + u = cheby.itransform(u_hat) + + weights = cheby.get_integration_weights() + integral = weights @ u_hat + + # generate a reference solution with qmat + lag = LagrangeApproximation(points=x) + Q = lag.getIntegrationMatrix(intervals=[(x0, x1)]) + integral_ref = (Q @ u)[0] + + assert np.isclose(integral, integral_ref) + + @pytest.mark.base @pytest.mark.parametrize('N', [4]) @pytest.mark.parametrize('d', [1, 2, 3]) diff --git a/pySDC/tests/test_helpers/test_spectral_helper_1d_fft.py b/pySDC/tests/test_helpers/test_spectral_helper_1d_fft.py index 4d456036c0..b1b7619aae 100644 --- a/pySDC/tests/test_helpers/test_spectral_helper_1d_fft.py +++ b/pySDC/tests/test_helpers/test_spectral_helper_1d_fft.py @@ -112,6 +112,36 @@ def test_integration_matrix(N, plot=False): assert np.allclose(expect, Du) +@pytest.mark.base +@pytest.mark.parametrize('x0', [-1, 0]) +@pytest.mark.parametrize('x1', [0.789, 1]) +@pytest.mark.parametrize('N', [32, 45]) +def test_integral_whole_interval(x0, x1, N): + import numpy as np + from pySDC.helpers.spectral_helper import FFTHelper + from qmat.lagrange import LagrangeApproximation + + helper = FFTHelper(N, x0=x0, x1=x1) + x = helper.get_1dgrid() + + u = np.zeros_like(x) + + num_coef = N // 2 - 1 + coeffs = np.random.random((2, N)) + u += coeffs[0, 0] + for i in range(1, num_coef + 1): + u += coeffs[0, i] * np.sin(2 * np.pi * i * x / helper.L) + u += coeffs[1, i] * np.cos(2 * np.pi * i * x / helper.L) + + u_hat = helper.transform(u) + + weights = helper.get_integration_weights() + integral = weights @ u_hat + integral_ref = coeffs[0, 0] * helper.L + + assert np.isclose(integral, integral_ref, atol=1e-7), abs(integral_ref - integral) + + @pytest.mark.base @pytest.mark.parametrize('N', [4, 32]) @pytest.mark.parametrize('v', [0, 4.78]) @@ -141,4 +171,5 @@ def test_tau_method(N, v): # test_tau_method(6, 1) # test_transform(True) # test_transform(False) - test_transform_cupy(4) + # test_transform_cupy(4) + test_integral_whole_interval(0, 2, 90)