diff --git a/wafo/integrate.py b/wafo/integrate.py index 5669a45..d1cb26d 100644 --- a/wafo/integrate.py +++ b/wafo/integrate.py @@ -8,7 +8,7 @@ from scipy import special as sp from wafo.plotbackend import plotbackend as plt from scipy.integrate import simps, trapz from wafo.demos import humps -from pychebfun import Chebfun +#from pychebfun import Chebfun _EPS = np.finfo(float).eps _POINTS_AND_WEIGHTS = {} @@ -1424,133 +1424,133 @@ def test_docstrings(): doctest.testmod() -def levin_integrate(): - ''' An oscillatory integral - Sheehan Olver, December 2010 - - - (Chebfun example quad/LevinIntegrate.m) - - This example computes the highly oscillatory integral of - - f * exp( 1i * w * g ), - - over (0,1) using the Levin method [1]. This method computes the integral - by rewriting it as an ODE - - u' + 1i * w * g' u = f, - - so that the indefinite integral of f * exp( 1i * w * g ) is - - u * exp( 1i * w * g ). - - - - We use as an example - - f = 1 / ( x + 2 ); - g = cos( x - 2 ); - w = 100000; - - # - References: - - [1] Levin, D., Procedures for computing one and two-dimensional integrals - of functions with rapid irregular oscillations, Maths Comp., 38 (1982) 531--538 - ''' - exp = np.exp - domain=[0, 1] - x = Chebfun.identity(domain=domain) - f = 1./(x+2) - g = np.cos(x-2) - D = np.diff(domain) - - - # Here is are plots of this integrand, with w = 100, in complex space - w = 100; - line_opts = dict(line_width=1.6) - font_opts = dict(font_size= 14) - # - - intg = f*exp(1j*w*g) - xs, ys, xi, yi, d = intg.plot_data(1000) - #intg.plot(with_interpolation_points=True) - #xi = np.linspace(0, 1, 1024) -# plt.plot(xs, ys) # , **line_opts) -# plt.plot(xi, yi, 'r.') +# def levin_integrate(): +# ''' An oscillatory integral +# Sheehan Olver, December 2010 +# +# +# (Chebfun example quad/LevinIntegrate.m) +# +# This example computes the highly oscillatory integral of +# +# f * exp( 1i * w * g ), +# +# over (0,1) using the Levin method [1]. This method computes the integral +# by rewriting it as an ODE +# +# u' + 1i * w * g' u = f, +# +# so that the indefinite integral of f * exp( 1i * w * g ) is +# +# u * exp( 1i * w * g ). +# +# +# +# We use as an example +# +# f = 1 / ( x + 2 ); +# g = cos( x - 2 ); +# w = 100000; +# +# # +# References: +# +# [1] Levin, D., Procedures for computing one and two-dimensional integrals +# of functions with rapid irregular oscillations, Maths Comp., 38 (1982) 531--538 +# ''' +# exp = np.exp +# domain=[0, 1] +# x = Chebfun.identity(domain=domain) +# f = 1./(x+2) +# g = np.cos(x-2) +# D = np.diff(domain) +# +# +# # Here is are plots of this integrand, with w = 100, in complex space +# w = 100; +# line_opts = dict(line_width=1.6) +# font_opts = dict(font_size= 14) +# # +# +# intg = f*exp(1j*w*g) +# xs, ys, xi, yi, d = intg.plot_data(1000) +# #intg.plot(with_interpolation_points=True) +# #xi = np.linspace(0, 1, 1024) +# # plt.plot(xs, ys) # , **line_opts) +# # plt.plot(xi, yi, 'r.') +# # #axis equal +# # plt.title('Complex plot of integrand') #,**font_opts) +# # plt.show('hold') +# ## +# # and of just the real part +# # intgr = np.real(intg) +# # xs, ys, xi, yi, d = intgr.plot_data(1000) +# #intgr.plot() +# # plt.plot(xs, np.real(ys)) # , **line_opts) +# # plt.plot(xi, np.real(yi), 'r.') # #axis equal -# plt.title('Complex plot of integrand') #,**font_opts) -# plt.show('hold') - ## - # and of just the real part -# intgr = np.real(intg) -# xs, ys, xi, yi, d = intgr.plot_data(1000) - #intgr.plot() -# plt.plot(xs, np.real(ys)) # , **line_opts) -# plt.plot(xi, np.real(yi), 'r.') - #axis equal -# plt.title('Real part of integrand') #,**font_opts) -# plt.show('hold') - - ## - # The Levin method will be accurate for large and small w, and the time - # taken is independent of w. Here we take a reasonably large value of w. - w = 1000; - intg = f*exp(1j*w*g) - val0 = np.sum(intg) - # val1 = sum(intg) - print(val0) - ## - # Start timing - #tic - - ## - # Construct the operator L - L = D + 1j*w*np.diag(g.differentiate()) - - ## - # From asymptotic analysis, we know that there exists a solution to the - # equation which is non-oscillatory, though we do not know what initial - # condition it satisfies. Thus we find a particular solution to this - # equation with no boundary conditions. - - u = L / f - - ## - # Because L is a differential operator with derivative order 1, \ expects - # it to be given a boundary condition, which is why the warning message is - # displayed. However, this doesn't cause any problems: though there are, - # in fact, a family of solutions to the ODE without boundary conditions - # due to the kernel - # - # exp(- 1i * w * g), - # - # it does not actually matter which particular solution is computed. - # Non-uniqueness is also not an issue: \ in matlab is least squares, hence - # does not require uniqueness. The existence of a non-oscillatory solution - # ensures that \ converges to a u with length independent of w. - # - # One could prevent the warning by applying a boundary condition consistent - # with the rest of the system, that is - # L.lbc = {L(1,:),f(0)}; - - ## - # Now we evaluate the antiderivative at the endpoints to obtain the - # integral. - - u(1)*exp(1j*w*g(1)) - u(0)*exp(1j*w*g(0)) - - #toc - - - ## - # Here is a way to compute the integral using Clenshaw--Curtis quadrature. - # As w becomes large, this takes an increasingly long time as the - # oscillations must be resolved. - - #tic - sum( f*exp(1j*w*g) ) - #toc +# # plt.title('Real part of integrand') #,**font_opts) +# # plt.show('hold') +# +# ## +# # The Levin method will be accurate for large and small w, and the time +# # taken is independent of w. Here we take a reasonably large value of w. +# w = 1000; +# intg = f*exp(1j*w*g) +# val0 = np.sum(intg) +# # val1 = sum(intg) +# print(val0) +# ## +# # Start timing +# #tic +# +# ## +# # Construct the operator L +# L = D + 1j*w*np.diag(g.differentiate()) +# +# ## +# # From asymptotic analysis, we know that there exists a solution to the +# # equation which is non-oscillatory, though we do not know what initial +# # condition it satisfies. Thus we find a particular solution to this +# # equation with no boundary conditions. +# +# u = L / f +# +# ## +# # Because L is a differential operator with derivative order 1, \ expects +# # it to be given a boundary condition, which is why the warning message is +# # displayed. However, this doesn't cause any problems: though there are, +# # in fact, a family of solutions to the ODE without boundary conditions +# # due to the kernel +# # +# # exp(- 1i * w * g), +# # +# # it does not actually matter which particular solution is computed. +# # Non-uniqueness is also not an issue: \ in matlab is least squares, hence +# # does not require uniqueness. The existence of a non-oscillatory solution +# # ensures that \ converges to a u with length independent of w. +# # +# # One could prevent the warning by applying a boundary condition consistent +# # with the rest of the system, that is +# # L.lbc = {L(1,:),f(0)}; +# +# ## +# # Now we evaluate the antiderivative at the endpoints to obtain the +# # integral. +# +# u(1)*exp(1j*w*g(1)) - u(0)*exp(1j*w*g(0)) +# +# #toc +# +# +# ## +# # Here is a way to compute the integral using Clenshaw--Curtis quadrature. +# # As w becomes large, this takes an increasingly long time as the +# # oscillations must be resolved. +# +# #tic +# sum( f*exp(1j*w*g) ) +# #toc # aLevinTQ[omega_,a_,b_,f_,g_,nu_,wprec_,prm_,test_,basis_,np_]:= @@ -1624,8 +1624,8 @@ def levin_integrate(): if __name__ == '__main__': - levin_integrate() - # test_docstrings() + # levin_integrate() + test_docstrings() # qdemo(np.exp, 0, 3, plot_error=True) # plt.show('hold') # main() diff --git a/wafo/test/test_padua.py b/wafo/test/test_padua.py index c9372fb..589a93c 100644 --- a/wafo/test/test_padua.py +++ b/wafo/test/test_padua.py @@ -4,33 +4,34 @@ import unittest import numpy as np from numpy import cos, pi from numpy.testing import assert_array_almost_equal -from wafo.padua import (padua_points, testfunct, padua_fit, +from wafo.padua import (padua_points, example_functions, padua_fit, padua_fit2, padua_cubature, padua_val) class PaduaTestCase(unittest.TestCase): + def test_padua_points_degree0(self): pad = padua_points(0) - expected = [[-1],[-1]] + expected = [[-1], [-1]] assert_array_almost_equal(pad, expected, 15) def test_padua_points_degree1(self): pad = padua_points(1) - expected = [cos(np.r_[0,1,1]*pi), - cos(np.r_[1,0,2]*pi/2)] + expected = [cos(np.r_[0, 1, 1] * pi), + cos(np.r_[1, 0, 2] * pi / 2)] assert_array_almost_equal(pad, expected, 15) def test_padua_points_degree2(self): - pad = padua_points(2, domain=[0,1,0,2]) - expected = [(cos(np.r_[0,0,1,1,2,2]*pi/2)+1)/2, - cos(np.r_[1,3,0,2,1,3]*pi/3)+1] + pad = padua_points(2, domain=[0, 1, 0, 2]) + expected = [(cos(np.r_[0, 0, 1, 1, 2, 2] * pi / 2) + 1) / 2, + cos(np.r_[1, 3, 0, 2, 1, 3] * pi / 3) + 1] assert_array_almost_equal(pad, expected, 15) def test_testfunct(self): - vals = [testfunct(0, 0, id) for id in range(12)] + vals = [example_functions(0, 0, id) for id in range(12)] expected = [7.664205912849231e-01, 0.7071067811865476, 0, 1.6487212707001282, 1.9287498479639178e-22, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0] @@ -38,38 +39,38 @@ class PaduaTestCase(unittest.TestCase): def test_padua_fit_even_degree(self): points = padua_points(10) - C0f, abs_error = padua_fit(points, testfunct, 6) + C0f, _abs_error = padua_fit(points, example_functions, 6) expected = np.zeros((11, 11)) - expected[0,0] = 1; + expected[0, 0] = 1 assert_array_almost_equal(C0f, expected, 15) def test_padua_fit_odd_degree(self): points = padua_points(9) - C0f, abs_error = padua_fit(points, testfunct, 6) + C0f, _abs_error = padua_fit(points, example_functions, 6) expected = np.zeros((10, 10)) - expected[0,0] = 1; + expected[0, 0] = 1 assert_array_almost_equal(C0f, expected, 15) def test_padua_fit_odd_degree2(self): points = padua_points(9) - C0f, abs_error = padua_fit2(points, testfunct, 6) + C0f, _abs_error = padua_fit2(points, example_functions, 6) expected = np.zeros((10, 10)) - expected[0,0] = 1; + expected[0, 0] = 1 assert_array_almost_equal(C0f, expected, 15) def test_padua_cubature(self): - domain = [0,1,0,1] + domain = [0, 1, 0, 1] points = padua_points(500, domain) - C0f, abs_error = padua_fit(points, testfunct, 0) + C0f, _abs_error = padua_fit(points, example_functions, 0) val = padua_cubature(C0f, domain) expected = 4.06969589491556e-01 assert_array_almost_equal(val, expected, 15) def test_padua_val_unordered(self): - domain = [0,1,0,1] + domain = [0, 1, 0, 1] points = padua_points(20, domain) - C0f, abs_error = padua_fit(points, testfunct, 0) - X = [0,0.5,1] + C0f, _abs_error = padua_fit(points, example_functions, 0) + X = [0, 0.5, 1] val = padua_val(X, X, C0f, domain) expected = [7.664205912849228e-01, 3.2621734202884815e-01, @@ -77,13 +78,13 @@ class PaduaTestCase(unittest.TestCase): assert_array_almost_equal(val, expected, 14) def test_padua_val_grid(self): - domain = [0,1,0,1] + domain = [0, 1, 0, 1] a, b, c, d = domain points = padua_points(21, domain) - C0f, abs_error = padua_fit(points, testfunct, 0) + C0f, _abs_error = padua_fit(points, example_functions, 0) X1 = np.linspace(a, b, 2) X2 = np.linspace(c, d, 2) - val = padua_val(X1, X2, C0f,domain, use_meshgrid=True); - expected = [[7.664205912849229e-01,1.0757071952145181e-01], - [2.703371615911344e-01,3.5734971024838565e-02]] + val = padua_val(X1, X2, C0f, domain, use_meshgrid=True) + expected = [[7.664205912849229e-01, 1.0757071952145181e-01], + [2.703371615911344e-01, 3.5734971024838565e-02]] assert_array_almost_equal(val, expected, 14)