diff --git a/wafo/integrate.py b/wafo/integrate.py index 02cca20..320ee29 100644 --- a/wafo/integrate.py +++ b/wafo/integrate.py @@ -1003,163 +1003,170 @@ def richardson(Q, k): return R -def quadgr(fun, a, b, abseps=1e-5, max_iter=17): - ''' - Gauss-Legendre quadrature with Richardson extrapolation. - - [Q,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function - FUN from A to B with an absolute error tolerance TOL. FUN is a function - handle and must accept vector arguments. TOL is 1e-6 by default. Q is - the integral approximation and ERR is an estimate of the absolute error. - - QUADGR uses a 12-point Gauss-Legendre quadrature. The error estimate is - based on successive interval bisection. Richardson extrapolation - accelerates the convergence for some integrals, especially integrals - with endpoint singularities. - - Examples - -------- - >>> import numpy as np - >>> Q, err = quadgr(np.log,0,1) - >>> quadgr(np.exp,0,9999*1j*np.pi) - (-2.0000000000122662, 2.1933237448479304e-09) - - >>> quadgr(lambda x: np.sqrt(4-x**2),0,2,1e-12) - (3.1415926535897811, 1.5809575870662229e-13) - - >>> quadgr(lambda x: x**-0.75,0,1) - (4.0000000000000266, 5.6843418860808015e-14) - - >>> quadgr(lambda x: 1./np.sqrt(1-x**2),-1,1) - (3.141596056985029, 6.2146261559092864e-06) - - >>> quadgr(lambda x: np.exp(-x**2),-np.inf,np.inf,1e-9) #% sqrt(pi) - (1.7724538509055152, 1.9722334876348668e-11) - - >>> quadgr(lambda x: np.cos(x)*np.exp(-x),0,np.inf,1e-9) - (0.50000000000000044, 7.3296813063450372e-11) +class _Quadgr(object): + + def __call__(self, fun, a, b, abseps=1e-5, max_iter=17): + ''' + Gauss-Legendre quadrature with Richardson extrapolation. + + [Q,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function + FUN from A to B with an absolute error tolerance TOL. FUN is a function + handle and must accept vector arguments. TOL is 1e-6 by default. Q is + the integral approximation and ERR is an estimate of the absolute + error. + + QUADGR uses a 12-point Gauss-Legendre quadrature. The error estimate is + based on successive interval bisection. Richardson extrapolation + accelerates the convergence for some integrals, especially integrals + with endpoint singularities. + + Examples + -------- + >>> import numpy as np + >>> Q, err = quadgr(np.log,0,1) + >>> quadgr(np.exp,0,9999*1j*np.pi) + (-2.0000000000122662, 2.1933237448479304e-09) + + >>> quadgr(lambda x: np.sqrt(4-x**2),0,2,1e-12) + (3.1415926535897811, 1.5809575870662229e-13) + + >>> quadgr(lambda x: x**-0.75,0,1) + (4.0000000000000266, 5.6843418860808015e-14) + + >>> quadgr(lambda x: 1./np.sqrt(1-x**2),-1,1) + (3.141596056985029, 6.2146261559092864e-06) + + >>> quadgr(lambda x: np.exp(-x**2),-np.inf,np.inf,1e-9) #% sqrt(pi) + (1.7724538509055152, 1.9722334876348668e-11) + + >>> quadgr(lambda x: np.cos(x)*np.exp(-x),0,np.inf,1e-9) + (0.50000000000000044, 7.3296813063450372e-11) + + See also + -------- + QUAD, + QUADGK + ''' + # Author: jonas.lundgren@saabgroup.com, 2009. license BSD + # Order limits (required if infinite limits) + a = np.asarray(a) + b = np.asarray(b) + if a == b: + Q = b - a + err = b - a + return Q, err + elif np.real(a) > np.real(b): + reverse = True + a, b = b, a + else: + reverse = False + + # Infinite limits + if np.isinf(a) | np.isinf(b): + # Check real limits + if ~ np.isreal(a) | ~np.isreal(b) | np.isnan(a) | np.isnan(b): + raise ValueError('Infinite intervals must be real.') + + # Change of variable + if np.isfinite(a) & np.isinf(b): + # a to inf + [Q, err] = quadgr(lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2, + 0, 1, abseps) + elif np.isinf(a) & np.isfinite(b): + # -inf to b + [Q, err] = quadgr(lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2, + -1, 0, abseps) + else: # -inf to inf + [Q1, err1] = quadgr(lambda t: fun(t / (1 - t)) / (1 - t) ** 2, + 0, 1, abseps / 2) + [Q2, err2] = quadgr(lambda t: fun(t / (1 + t)) / (1 + t) ** 2, + -1, 0, abseps / 2) + Q = Q1 + Q2 + err = err1 + err2 + + # Reverse direction + if reverse: + Q = -Q + return Q, err + + # Gauss-Legendre quadrature (12-point) + xq = np.asarray( + [0.12523340851146894, 0.36783149899818018, 0.58731795428661748, + 0.76990267419430469, 0.9041172563704748, 0.98156063424671924]) + wq = np.asarray( + [0.24914704581340288, 0.23349253653835478, 0.20316742672306584, + 0.16007832854334636, 0.10693932599531818, 0.047175336386511842]) + xq = np.hstack((xq, -xq)) + wq = np.hstack((wq, wq)) + nq = len(xq) + dtype = np.result_type(fun(a), fun(b)) + + # Initiate vectors + Q0 = zeros(max_iter, dtype=dtype) # Quadrature + Q1 = zeros(max_iter, dtype=dtype) # First Richardson extrapolation + Q2 = zeros(max_iter, dtype=dtype) # Second Richardson extrapolation + + # One interval + hh = (b - a) / 2 # Half interval length + x = (a + b) / 2 + hh * xq # Nodes + # Quadrature + Q0[0] = hh * np.sum(wq * fun(x), axis=0) + + # Successive bisection of intervals + for k in range(1, max_iter): + + # Interval bisection + hh = hh / 2 + x = np.hstack([x + a, x + b]) / 2 + # Quadrature + Q0[k] = hh * np.sum(wq * np.sum(np.reshape(fun(x), (-1, nq)), + axis=0), + axis=0) + + # Richardson extrapolation + if k >= 5: + Q1[k] = richardson(Q0, k) + Q2[k] = richardson(Q1, k) + elif k >= 3: + Q1[k] = richardson(Q0, k) + + # Estimate absolute error + if k >= 6: + Qv = np.hstack((Q0[k], Q1[k], Q2[k])) + Qw = np.hstack((Q0[k - 1], Q1[k - 1], Q2[k - 1])) + elif k >= 4: + Qv = np.hstack((Q0[k], Q1[k])) + Qw = np.hstack((Q0[k - 1], Q1[k - 1])) + else: + Qv = np.atleast_1d(Q0[k]) + Qw = Q0[k - 1] + + errors = np.atleast_1d(abs(Qv - Qw)) + j = errors.argmin() + err = errors[j] + Q = Qv[j] + if k >= 2: # and not iscomplex: + _val, err1 = dea3(Q0[k - 2], Q0[k - 1], Q0[k]) + + # Convergence + if (err < abseps) | ~np.isfinite(Q): + break + else: + warnings.warn('Max number of iterations reached without ' + + 'convergence.') - See also - -------- - QUAD, - QUADGK - ''' - # Author: jonas.lundgren@saabgroup.com, 2009. license BSD - # Order limits (required if infinite limits) - a = np.asarray(a) - b = np.asarray(b) - if a == b: - Q = b - a - err = b - a - return Q, err - elif np.real(a) > np.real(b): - reverse = True - a, b = b, a - else: - reverse = False - - # Infinite limits - if np.isinf(a) | np.isinf(b): - # Check real limits - if ~ np.isreal(a) | ~np.isreal(b) | np.isnan(a) | np.isnan(b): - raise ValueError('Infinite intervals must be real.') - - # Change of variable - if np.isfinite(a) & np.isinf(b): - # a to inf - [Q, err] = quadgr(lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2, - 0, 1, abseps) - elif np.isinf(a) & np.isfinite(b): - # -inf to b - [Q, err] = quadgr(lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2, - -1, 0, abseps) - else: # -inf to inf - [Q1, err1] = quadgr(lambda t: fun(t / (1 - t)) / (1 - t) ** 2, - 0, 1, abseps / 2) - [Q2, err2] = quadgr(lambda t: fun(t / (1 + t)) / (1 + t) ** 2, - -1, 0, abseps / 2) - Q = Q1 + Q2 - err = err1 + err2 + if ~ np.isfinite(Q): + warnings.warn('Integral approximation is Infinite or NaN.') + # The error estimate should not be zero + err = err + 2 * np.finfo(Q).eps # Reverse direction if reverse: Q = -Q - return Q, err - # Gauss-Legendre quadrature (12-point) - xq = np.asarray( - [0.12523340851146894, 0.36783149899818018, 0.58731795428661748, - 0.76990267419430469, 0.9041172563704748, 0.98156063424671924]) - wq = np.asarray( - [0.24914704581340288, 0.23349253653835478, 0.20316742672306584, - 0.16007832854334636, 0.10693932599531818, 0.047175336386511842]) - xq = np.hstack((xq, -xq)) - wq = np.hstack((wq, wq)) - nq = len(xq) - dtype = np.result_type(fun(a), fun(b)) - - # Initiate vectors - Q0 = zeros(max_iter, dtype=dtype) # Quadrature - Q1 = zeros(max_iter, dtype=dtype) # First Richardson extrapolation - Q2 = zeros(max_iter, dtype=dtype) # Second Richardson extrapolation - - # One interval - hh = (b - a) / 2 # Half interval length - x = (a + b) / 2 + hh * xq # Nodes - # Quadrature - Q0[0] = hh * np.sum(wq * fun(x), axis=0) - - # Successive bisection of intervals - for k in range(1, max_iter): - - # Interval bisection - hh = hh / 2 - x = np.hstack([x + a, x + b]) / 2 - # Quadrature - Q0[k] = hh * np.sum(wq * np.sum(np.reshape(fun(x), (-1, nq)), axis=0), - axis=0) - - # Richardson extrapolation - if k >= 5: - Q1[k] = richardson(Q0, k) - Q2[k] = richardson(Q1, k) - elif k >= 3: - Q1[k] = richardson(Q0, k) - - # Estimate absolute error - if k >= 6: - Qv = np.hstack((Q0[k], Q1[k], Q2[k])) - Qw = np.hstack((Q0[k - 1], Q1[k - 1], Q2[k - 1])) - elif k >= 4: - Qv = np.hstack((Q0[k], Q1[k])) - Qw = np.hstack((Q0[k - 1], Q1[k - 1])) - else: - Qv = np.atleast_1d(Q0[k]) - Qw = Q0[k - 1] - - errors = np.atleast_1d(abs(Qv - Qw)) - j = errors.argmin() - err = errors[j] - Q = Qv[j] - if k >= 2: # and not iscomplex: - _val, err1 = dea3(Q0[k - 2], Q0[k - 1], Q0[k]) - - # Convergence - if (err < abseps) | ~np.isfinite(Q): - break - else: - warnings.warn('Max number of iterations reached without convergence.') - - if ~ np.isfinite(Q): - warnings.warn('Integral approximation is Infinite or NaN.') - - # The error estimate should not be zero - err = err + 2 * np.finfo(Q).eps - # Reverse direction - if reverse: - Q = -Q + return Q, err - return Q, err +quadgr = _Quadgr() def boole(y, x): diff --git a/wafo/source/c_library/c_functions.c b/wafo/source/c_library/c_functions.c index a9bbdd5..d7df721 100644 --- a/wafo/source/c_library/c_functions.c +++ b/wafo/source/c_library/c_functions.c @@ -1,11 +1,11 @@ -#include "math.h" +#include "math.h" /* * Install gfortran and run the following to build the module on windows: * f2py c_library.pyf c_functions.c -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 */ - + /* - * findrfc.c - + * findrfc.c - * * Returns indices to RFC turningpoints of a vector * of turningpoints @@ -17,9 +17,9 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) { double xminus,xplus,Tpl,Tmi,*y,Tstart; int i,j,ix=0,NC,iy; info[0] = 0; - if (*(y1+0)> *(y1+1)){ + if (*(y1+0)> *(y1+1)){ /* if first is a max , ignore the first max*/ - y=&(*(y1+1)); + y=&(*(y1+1)); NC=floor((n-1)/2); Tstart=1; } @@ -28,85 +28,85 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) { NC=floor(n/2); Tstart=0; } - + if (NC<1){ return; /* No RFC cycles*/ } - + if (( *(y+0) > *(y+1)) && ( *(y+1) > *(y+2)) ){ - info[0] = -1; + info[0] = -1; return; /*This is not a sequence of turningpoints, exit */ } if ((*(y+0) < *(y+1)) && (*(y+1)< *(y+2))){ info[0]=-1; return; /*This is not a sequence of turningpoints, exit */ } - - + + for (i=0; i=0) && (*(y+2*j+1)<=*(y+2*i+1))){ - if( (*(y+2*j)=0) && (*(y+2*j+1)<=*(y+2*i+1))){ + if( (*(y+2*j)= xplus){ - if ( (*(y+2*i+1)-xminus) >= hmin){ - *(ind+ix)=Tmi; - ix++; - *(ind+ix)=(Tstart+2*i+1); - ix++; - } /*if*/ - goto L180; - } - + if ( (*(y+2*i+1)-xminus) >= hmin){ + *(ind+ix)=Tmi; + ix++; + *(ind+ix)=(Tstart+2*i+1); + ix++; + } /*if*/ + goto L180; + } + j=i+1; while((j= *(y+2*i+1)) goto L170; - if( (*(y+2*j+2) <= xplus) ){ - xplus=*(y+2*j+2); - Tpl=(Tstart+2*j+2); - }/*if*/ - j++; + if (*(y+2*j+1) >= *(y+2*i+1)) goto L170; + if( (*(y+2*j+2) <= xplus) ){ + xplus=*(y+2*j+2); + Tpl=(Tstart+2*j+2); + }/*if*/ + j++; } /*while*/ - + if ( (*(y+2*i+1)-xminus) >= hmin) { - *(ind+ix)=Tmi; - ix++; - *(ind+ix)=(Tstart+2*i+1); - ix++; - + *(ind+ix)=Tmi; + ix++; + *(ind+ix)=(Tstart+2*i+1); + ix++; + } /*if*/ goto L180; - L170: + L170: if (xplus <= xminus ) { - if ( (*(y+2*i+1)-xminus) >= hmin){ - *(ind+ix)=Tmi; - ix++; - *(ind+ix)=(Tstart+2*i+1); - ix++; - } /*if*/ - /*goto L180;*/ + if ( (*(y+2*i+1)-xminus) >= hmin){ + *(ind+ix)=Tmi; + ix++; + *(ind+ix)=(Tstart+2*i+1); + ix++; + } /*if*/ + /*goto L180;*/ } - else{ - if ( (*(y+2*i+1)-xplus) >= hmin) { - *(ind+ix)=(Tstart+2*i+1); - ix++; - *(ind+ix)=Tpl; - ix++; - } /*if*/ + else{ + if ( (*(y+2*i+1)-xplus) >= hmin) { + *(ind+ix)=(Tstart+2*i+1); + ix++; + *(ind+ix)=Tpl; + ix++; + } /*if*/ } /*elseif*/ L180: iy=i; @@ -118,52 +118,52 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) { /* - * findcross.c - + * findcross.c - * * Returns indices to level v crossings of argument vector * * 1998 by Per Andreas Brodtkorb. last modified 23.06-98 - */ + */ void findcross(double *y, double v, int *ind, int n, int *info) { int i,start, ix=0,dcross=0; - start=0; - if ( y[0]< v){ - dcross=-1; /* first is a up-crossing*/ - } - else if ( y[0]> v){ - dcross=1; /* first is a down-crossing*/ - } - else if ( y[0]== v){ - /* Find out what type of crossing we have next time.. */ - for (i=1; i v){ - ind[ix] = i-1; /* first crossing is a up-crossing*/ - ix++; - dcross=1; /*The next crossing is a down-crossing*/ - goto L120; - } - } - } - L120: - for (i=start; i v) ) || ((dcross==1 ) && (y[i]>=v) && (y[i+1] < v) ) ) { - - ind[ix] = i; - ix++; - dcross=-dcross; - } - } - info[0] = ix; - return; + start=0; + if ( y[0]< v){ + dcross=-1; /* first is a up-crossing*/ + } + else if ( y[0]> v){ + dcross=1; /* first is a down-crossing*/ + } + else if ( y[0]== v){ + /* Find out what type of crossing we have next time.. */ + for (i=1; i v){ + ind[ix] = i-1; /* first crossing is a up-crossing*/ + ix++; + dcross=1; /*The next crossing is a down-crossing*/ + goto L120; + } + } + } + L120: + for (i=start; i v) ) || ((dcross==1 ) && (y[i]>=v) && (y[i+1] < v) ) ) { + + ind[ix] = i; + ix++; + dcross=-dcross; + } + } + info[0] = ix; + return; } @@ -173,19 +173,19 @@ void findcross(double *y, double v, int *ind, int n, int *info) * CALL: disufq(rvec,ivec,rA,iA, w,kw,h,g,nmin,nmax,m,n) * * rvec, ivec = real and imaginary parts of the resultant (size m X n). - * rA, iA = real and imaginary parts of the amplitudes (size m X n). + * rA, iA = real and imaginary parts of the amplitudes (size m X n). * w = vector with angular frequencies (w>=0) * kw = vector with wavenumbers (kw>=0) * h = water depth (h >=0) * g = constant acceleration of gravity - * nmin = minimum index where rA(:,nmin) and iA(:,nmin) is + * nmin = minimum index where rA(:,nmin) and iA(:,nmin) is * greater than zero. - * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is + * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is * greater than zero. * m = size(rA,1),size(iA,1) * n = size(rA,2),size(iA,2), or size(rvec,2),size(ivec,2) * - * DISUFQ returns the summation of difference frequency and sum + * DISUFQ returns the summation of difference frequency and sum * frequency effects in the vector vec = rvec +sqrt(-1)*ivec. * The 2'nd order contribution to the Stokes wave is then calculated by * a simple 1D Fourier transform, real(FFT(vec)). @@ -195,15 +195,15 @@ void findcross(double *y, double v, int *ind, int n, int *info) * * by Per Andreas Brodtkorb 15.08.2001 * revised pab 14.03.2002, 01.05.2002 22.07.2002, oct 2008 - */ + */ void disufq(double *rvec, double *ivec, - double *rA, double *iA, - double *w, double *kw, - double h, double g, - int nmin, int nmax, - int m, int n) -{ + double *rA, double *iA, + double *w, double *kw, + double h, double g, + int nmin, int nmax, + int m, int n) +{ double Epij, Edij; double tmp1, tmp2, tmp3, tmp4, kfact; double w1, w2, kw1, kw2, Cg; @@ -217,7 +217,7 @@ void disufq(double *rvec, double *ivec, } // kfact is set to 2 in order to exploit the symmetry. - // If you set kfact to 1, you must uncomment all statements + // If you set kfact to 1, you must uncomment all statements // including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2]. kfact = 2.0; @@ -227,70 +227,70 @@ void disufq(double *rvec, double *ivec, iz1 = 2*ixi; //iz2 = n*m-ixi; kw1 = kw[ix]; - Epij = kw1; - for (i=0;i=0) * kw = vector with wavenumbers (kw>=0) * h = water depth (h >=0) * g = constant acceleration of gravity - * nmin = minimum index where rA(:,nmin) and iA(:,nmin) is + * nmin = minimum index where rA(:,nmin) and iA(:,nmin) is * greater than zero. - * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is + * nmax = maximum index where rA(:,nmax) and iA(:,nmax) is * greater than zero. * m = size(rA,1),size(iA,1) * n = size(rA,2),size(iA,2), or size(rvec,2),size(ivec,2) * - * DISUFQ2 returns the summation of sum and difference frequency + * DISUFQ2 returns the summation of sum and difference frequency * frequency effects in the vectors svec = rsvec +sqrt(-1)*isvec and * dvec = rdvec +sqrt(-1)*idvec. * The 2'nd order contribution to the Stokes wave is then calculated by @@ -416,17 +416,17 @@ void disufq(double *rvec, double *ivec, * * This is a MEX-file for MATLAB. * by Per Andreas Brodtkorb 15.08.2001 - * revised pab 14.03.2002, 01.05.2002 - */ + * revised pab 14.03.2002, 01.05.2002 + */ void disufq2(double *rsvec, double *isvec, - double *rdvec, double *idvec, - double *rA, double *iA, - double *w, double *kw, - double h, double g, - int nmin, int nmax, - int m, int n) -{ + double *rdvec, double *idvec, + double *rA, double *iA, + double *w, double *kw, + double h, double g, + int nmin, int nmax, + int m, int n) +{ double Epij, Edij; double tmp1, tmp2, tmp3, tmp4, kfact; double w1, w2, kw1, kw2, Cg; @@ -443,7 +443,7 @@ void disufq2(double *rsvec, double *isvec, } // kfact is set to 2 in order to exploit the symmetry. - // If you set kfact to 1, you must uncomment all statements + // If you set kfact to 1, you must uncomment all statements // including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2]. kfact = 2.0; @@ -453,70 +453,70 @@ void disufq2(double *rsvec, double *isvec, iz1 = 2*ixi; //iz2 = n*m-ixi; kw1 = kw[ix]; - Epij = kw1; - for (i=0;i