danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Removed deprecated import

Browse Source 
Updated test/test_padua.py
master
pbrod 9 years ago
parent dd9cfb63ad
commit e6d51a1fab
2 changed files with 154 additions and 153 deletions
Show Stats Download Patch File Download Diff File

258
wafo/integrate.py
Unescape Escape View File

@ -8,7 +8,7 @@ from scipy import special as sp
from wafo.plotbackend import plotbackend as plt from wafo.plotbackend import plotbackend as plt
from scipy.integrate import simps, trapz from scipy.integrate import simps, trapz
from wafo.demos import humps from wafo.demos import humps
from pychebfun import Chebfun #from pychebfun import Chebfun
_EPS = np.finfo(float).eps _EPS = np.finfo(float).eps
_POINTS_AND_WEIGHTS = {} _POINTS_AND_WEIGHTS = {}
@ -1424,133 +1424,133 @@ def test_docstrings():
doctest.testmod() doctest.testmod()
def levin_integrate(): # def levin_integrate():
''' An oscillatory integral # ''' An oscillatory integral
Sheehan Olver, December 2010 # Sheehan Olver, December 2010
#
#
(Chebfun example quad/LevinIntegrate.m) # (Chebfun example quad/LevinIntegrate.m)
#
This example computes the highly oscillatory integral of # This example computes the highly oscillatory integral of
#
f * exp( 1i * w * g ), # f * exp( 1i * w * g ),
#
over (0,1) using the Levin method [1]. This method computes the integral # over (0,1) using the Levin method [1]. This method computes the integral
by rewriting it as an ODE # by rewriting it as an ODE
#
u' + 1i * w * g' u = f, # u' + 1i * w * g' u = f,
#
so that the indefinite integral of f * exp( 1i * w * g ) is # so that the indefinite integral of f * exp( 1i * w * g ) is
#
u * exp( 1i * w * g ). # u * exp( 1i * w * g ).
#
#
#
We use as an example # We use as an example
#
f = 1 / ( x + 2 ); # f = 1 / ( x + 2 );
g = cos( x - 2 ); # g = cos( x - 2 );
w = 100000; # w = 100000;
#
# # #
References: # References:
#
[1] Levin, D., Procedures for computing one and two-dimensional integrals # [1] Levin, D., Procedures for computing one and two-dimensional integrals
of functions with rapid irregular oscillations, Maths Comp., 38 (1982) 531--538 # of functions with rapid irregular oscillations, Maths Comp., 38 (1982) 531--538
''' # '''
exp = np.exp # exp = np.exp
domain=[0, 1] # domain=[0, 1]
x = Chebfun.identity(domain=domain) # x = Chebfun.identity(domain=domain)
f = 1./(x+2) # f = 1./(x+2)
g = np.cos(x-2) # g = np.cos(x-2)
D = np.diff(domain) # D = np.diff(domain)
#
#
# Here is are plots of this integrand, with w = 100, in complex space # # Here is are plots of this integrand, with w = 100, in complex space
w = 100; # w = 100;
line_opts = dict(line_width=1.6) # line_opts = dict(line_width=1.6)
font_opts = dict(font_size= 14) # font_opts = dict(font_size= 14)
# # #
#
intg = f*exp(1j*w*g) # intg = f*exp(1j*w*g)
xs, ys, xi, yi, d = intg.plot_data(1000) # xs, ys, xi, yi, d = intg.plot_data(1000)
#intg.plot(with_interpolation_points=True) # #intg.plot(with_interpolation_points=True)
#xi = np.linspace(0, 1, 1024) # #xi = np.linspace(0, 1, 1024)
# plt.plot(xs, ys) # , **line_opts) # # plt.plot(xs, ys) # , **line_opts)
# plt.plot(xi, yi, 'r.') # # 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 # #axis equal
# plt.title('Complex plot of integrand') #,**font_opts) # # plt.title('Real part of integrand') #,**font_opts)
# plt.show('hold') # # plt.show('hold')
## #
# and of just the real part # ##
# intgr = np.real(intg) # # The Levin method will be accurate for large and small w, and the time
# xs, ys, xi, yi, d = intgr.plot_data(1000) # # taken is independent of w. Here we take a reasonably large value of w.
#intgr.plot() # w = 1000;
# plt.plot(xs, np.real(ys)) # , **line_opts) # intg = f*exp(1j*w*g)
# plt.plot(xi, np.real(yi), 'r.') # val0 = np.sum(intg)
#axis equal # # val1 = sum(intg)
# plt.title('Real part of integrand') #,**font_opts) # print(val0)
# plt.show('hold') # ##
# # Start timing
## # #tic
# 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; # # Construct the operator L
intg = f*exp(1j*w*g) # L = D + 1j*w*np.diag(g.differentiate())
val0 = np.sum(intg) #
# val1 = sum(intg) # ##
print(val0) # # From asymptotic analysis, we know that there exists a solution to the
## # # equation which is non-oscillatory, though we do not know what initial
# Start timing # # condition it satisfies. Thus we find a particular solution to this
#tic # # equation with no boundary conditions.
#
## # u = L / f
# Construct the operator L #
L = D + 1j*w*np.diag(g.differentiate()) # ##
# # 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
# From asymptotic analysis, we know that there exists a solution to the # # displayed. However, this doesn't cause any problems: though there are,
# equation which is non-oscillatory, though we do not know what initial # # in fact, a family of solutions to the ODE without boundary conditions
# condition it satisfies. Thus we find a particular solution to this # # due to the kernel
# equation with no boundary conditions. # #
# # exp(- 1i * w * g),
u = L / f # #
# # it does not actually matter which particular solution is computed.
## # # Non-uniqueness is also not an issue: \ in matlab is least squares, hence
# Because L is a differential operator with derivative order 1, \ expects # # does not require uniqueness. The existence of a non-oscillatory solution
# it to be given a boundary condition, which is why the warning message is # # ensures that \ converges to a u with length independent of w.
# displayed. However, this doesn't cause any problems: though there are, # #
# in fact, a family of solutions to the ODE without boundary conditions # # One could prevent the warning by applying a boundary condition consistent
# due to the kernel # # with the rest of the system, that is
# # # L.lbc = {L(1,:),f(0)};
# exp(- 1i * w * g), #
# # ##
# it does not actually matter which particular solution is computed. # # Now we evaluate the antiderivative at the endpoints to obtain the
# Non-uniqueness is also not an issue: \ in matlab is least squares, hence # # integral.
# does not require uniqueness. The existence of a non-oscillatory solution #
# ensures that \ converges to a u with length independent of w. # u(1)*exp(1j*w*g(1)) - u(0)*exp(1j*w*g(0))
# #
# One could prevent the warning by applying a boundary condition consistent # #toc
# with the rest of the system, that is #
# L.lbc = {L(1,:),f(0)}; #
# ##
## # # Here is a way to compute the integral using Clenshaw--Curtis quadrature.
# Now we evaluate the antiderivative at the endpoints to obtain the # # As w becomes large, this takes an increasingly long time as the
# integral. # # oscillations must be resolved.
#
u(1)*exp(1j*w*g(1)) - u(0)*exp(1j*w*g(0)) # #tic
# sum( f*exp(1j*w*g) )
#toc # #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_]:= # aLevinTQ[omega_,a_,b_,f_,g_,nu_,wprec_,prm_,test_,basis_,np_]:=
@ -1624,8 +1624,8 @@ def levin_integrate():
if __name__ == '__main__': if __name__ == '__main__':
levin_integrate() # levin_integrate()
# test_docstrings() test_docstrings()
# qdemo(np.exp, 0, 3, plot_error=True) # qdemo(np.exp, 0, 3, plot_error=True)
# plt.show('hold') # plt.show('hold')
# main() # main()

49
wafo/test/test_padua.py
Unescape Escape View File

@ -4,33 +4,34 @@ import unittest
import numpy as np import numpy as np
from numpy import cos, pi from numpy import cos, pi
from numpy.testing import assert_array_almost_equal 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_fit2,
padua_cubature, padua_val) padua_cubature, padua_val)
class PaduaTestCase(unittest.TestCase): class PaduaTestCase(unittest.TestCase):
def test_padua_points_degree0(self): def test_padua_points_degree0(self):
pad = padua_points(0) pad = padua_points(0)
expected = [[-1],[-1]] expected = [[-1], [-1]]
assert_array_almost_equal(pad, expected, 15) assert_array_almost_equal(pad, expected, 15)
def test_padua_points_degree1(self): def test_padua_points_degree1(self):
pad = padua_points(1) pad = padua_points(1)
expected = [cos(np.r_[0,1,1]*pi), expected = [cos(np.r_[0, 1, 1] * pi),
cos(np.r_[1,0,2]*pi/2)] cos(np.r_[1, 0, 2] * pi / 2)]
assert_array_almost_equal(pad, expected, 15) assert_array_almost_equal(pad, expected, 15)
def test_padua_points_degree2(self): def test_padua_points_degree2(self):
pad = padua_points(2, domain=[0,1,0,2]) pad = padua_points(2, domain=[0, 1, 0, 2])
expected = [(cos(np.r_[0,0,1,1,2,2]*pi/2)+1)/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] cos(np.r_[1, 3, 0, 2, 1, 3] * pi / 3) + 1]
assert_array_almost_equal(pad, expected, 15) assert_array_almost_equal(pad, expected, 15)
def test_testfunct(self): 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, expected = [7.664205912849231e-01, 0.7071067811865476, 0,
1.6487212707001282, 1.9287498479639178e-22, 1.0, 1.6487212707001282, 1.9287498479639178e-22, 1.0,
1.0, 1.0, 1.0, 0.0, 1.0, 0.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): def test_padua_fit_even_degree(self):
points = padua_points(10) 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 = np.zeros((11, 11))
expected[0,0] = 1; expected[0, 0] = 1
assert_array_almost_equal(C0f, expected, 15) assert_array_almost_equal(C0f, expected, 15)
def test_padua_fit_odd_degree(self): def test_padua_fit_odd_degree(self):
points = padua_points(9) 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 = np.zeros((10, 10))
expected[0,0] = 1; expected[0, 0] = 1
assert_array_almost_equal(C0f, expected, 15) assert_array_almost_equal(C0f, expected, 15)
def test_padua_fit_odd_degree2(self): def test_padua_fit_odd_degree2(self):
points = padua_points(9) 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 = np.zeros((10, 10))
expected[0,0] = 1; expected[0, 0] = 1
assert_array_almost_equal(C0f, expected, 15) assert_array_almost_equal(C0f, expected, 15)
def test_padua_cubature(self): def test_padua_cubature(self):
domain = [0,1,0,1] domain = [0, 1, 0, 1]
points = padua_points(500, domain) 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) val = padua_cubature(C0f, domain)
expected = 4.06969589491556e-01 expected = 4.06969589491556e-01
assert_array_almost_equal(val, expected, 15) assert_array_almost_equal(val, expected, 15)
def test_padua_val_unordered(self): def test_padua_val_unordered(self):
domain = [0,1,0,1] domain = [0, 1, 0, 1]
points = padua_points(20, domain) points = padua_points(20, domain)
C0f, abs_error = padua_fit(points, testfunct, 0) C0f, _abs_error = padua_fit(points, example_functions, 0)
X = [0,0.5,1] X = [0, 0.5, 1]
val = padua_val(X, X, C0f, domain) val = padua_val(X, X, C0f, domain)
expected = [7.664205912849228e-01, 3.2621734202884815e-01, expected = [7.664205912849228e-01, 3.2621734202884815e-01,
@ -77,13 +78,13 @@ class PaduaTestCase(unittest.TestCase):
assert_array_almost_equal(val, expected, 14) assert_array_almost_equal(val, expected, 14)
def test_padua_val_grid(self): def test_padua_val_grid(self):
domain = [0,1,0,1] domain = [0, 1, 0, 1]
a, b, c, d = domain a, b, c, d = domain
points = padua_points(21, 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) X1 = np.linspace(a, b, 2)
X2 = np.linspace(c, d, 2) X2 = np.linspace(c, d, 2)
val = padua_val(X1, X2, C0f,domain, use_meshgrid=True); val = padua_val(X1, X2, C0f, domain, use_meshgrid=True)
expected = [[7.664205912849229e-01,1.0757071952145181e-01], expected = [[7.664205912849229e-01, 1.0757071952145181e-01],
[2.703371615911344e-01,3.5734971024838565e-02]] [2.703371615911344e-01, 3.5734971024838565e-02]]
assert_array_almost_equal(val, expected, 14) assert_array_almost_equal(val, expected, 14)

Write Preview
Loading…
Cancel
Save