(PC12xx~14xx) qthsh TanhSinh quadrature

03282021, 02:02 PM
(This post was last modified: 10142021 12:13 AM by robve.)
Post: #1




(PC12xx~14xx) qthsh TanhSinh quadrature
TANHSINH QUADRATURE and COMBINED TANHSINH, EXPSINH, SINHSINH QUADRATURES
Estimate definite integral over open interval with improved and optimized TanhSinh quadrature Credit: Michalski & Mosig TanhSinh rule, improved in this thread and in the accompanying technical paper For details, see the technical paper: https://www.genivia.com/files/qthsh.pdf See also: https://newtonexcelbach.com/2020/10/29/n...urev50/ For SHARP PC12xx to 14xx 381 bytes BASIC image (PC1350) "qthsh" optimized TanhSinh 952 bytes BASIC image (PC1350) "quad" optimized combined TanhSinh, ExpSinh, SinhSinh Precision: 1E9 (adjustable) Levels: N=6 (adjustable) Example: 300 "F1" Y=4/(1+X*X): RETURN RUN f=F1 a=0 b=1 3.141592653 310 "F2" Y=ATN(1/X): RETURN RADIAN RUN f=F2 a=0 b=1 1.131971754 320 "F3" Y=COS X*LN X: RETURN RADIAN RUN f=F3 a=0 b=1 9.460830706E01 330 "F4" Y=X^.8: RETURN RUN f=F4 a=0 b=1 5 Code: ' TANHSINH QUADRATURE Code: ' COMBINED TANHSINH, EXPSINH, SINHSINH QUADRATURE Edit: major update with the latest version. "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03282021, 06:59 PM
(This post was last modified: 03282021 07:00 PM by Albert Chan.)
Post: #2




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
Some issues I noticed for TanhSinh quadrature.
Let x = tanh(pi/2*sinh(t)) ∫(f(x), x=0 .. 1) = pi/2 * ∫(f(tanh(pi/2*sinh(t))) * (cosh(t)/cosh(pi*sinh(t)/2)^2), t = 0 .. inf) Numerically, t does not need to go too high (certainly not inf) lua> function ts(t) return tanh(pi/2*sinh(t)) end lua> ts(2.8), ts(3), ts(3.2) 0.9999999999866902 0.999999999999957 1 At t=3, with double precision, we are very close to evaluating the endpoint, f(1) If f is already utransformed, say, from f0, t is even smaller (to avoid touching f0(1)) ∫(f0(x), x=0..1) = ∫(6*u*(1u) * f0(u*u*(32*u)), u=0..1) // RHS integrand = f(u) lua> function u(u) return u*u*(32*u) end lua> u(ts(2.2)), u(ts(2.4)), u(ts(2.6)) 0.9999999999917426 0.9999999999999855 1 So, what is the big deal ? If f(1) = ∞, utransformation might make the algorithm fail. >>> from mpmath import * >>> f = lambda x: 1/sqrt(ln(x)) >>> quadts(f, [0,1], error=True) (mpf('1.7724538503750329'), mpf('1.0e10')) >>> >>> uf = lambda u: 6*u*(1u) * f(u*u*(32*u)) >>> quadts(uf, [0,1], error=True) Traceback (most recent call last): ... ZeroDivisionError 

03292021, 01:10 AM
Post: #3




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(03282021 06:59 PM)Albert Chan Wrote: Some issues I noticed for TanhSinh quadrature. Thanks for sharing your findings! It's great to play with these methods to see their practical strengths and weaknesses. There are limitations to the applicability of the method that are generally known (see e.g. Boost Math TanhSinh for a decent explanation) and issues that are caused in practice using this method. But I am sure I don't need to tell you that The implementation in C and BASIC I wrote are based on several other implementations (see links). I have not optimized the code or tried to fix problems that I already saw in the code early on (more about this further below). (03282021 06:59 PM)Albert Chan Wrote: >>> from mpmath import * The crash with this function does not happen in the C and BASIC versions I've worked on (here output with k, error, and points for debugging): k=3 err=2.29154e08 pts=32 qthsh = 1.77245384129412 However, it happens to the same function after simply reversing the x range over [0,1]: double f(double x) { return 1/sqrt(log(1x)); } This crashes at k=2: k=2 err=nan pts=17 int qthsh = inf When prepping my implementation, I noticed the following conditions in some existing implementations of the method I compared to look for certain details such as this, which I've adopted in the exact same way to avoid hitting the endpoints: if (a+x > a) fp = f(a+x); if (bx < b) fm = f(bx); Do you spot the problem too? If a=0 or b=0 then x can be arbitrarily small. On the other hand, if a=1 or b=1 then this "fix" works the way it was intended. This happens in the "best TanhSinh implementation on the planet". Further, perhaps in the Python and Lua implementations there is no guard at all, since your case failed either way? Fortunately, we can simply modify the code to something like this: if (1+x > 1) { fp = f(a+x); fm = f(bx); } Hopefully the C/C++ compiler does not optimize the condition to x > 0, which is algebraically identical, but numerically destroys the benefits of having this test (because x > 0 is always true). The type of optimizations were a serious problem in the past. Compiler optimizations wrongly assumed that floating point arithmetic is associative and commutative. Modern C/C++ compilers should not mess this up. With this simple change I propose, it works as expected: k=3 err=2.29154e08 pts=29 int qthsh = 1.77245384129412 I can't think of a problem with this change, though the value 1 could be picked smaller such as eps or MachEps for example. It depends on how close you allow x to approach the endpoints. I also noticed in the public implementations of the method with this incomplete fix that fp and fm are simply the old f(a+x) and f(bx) values i.e. at neighboring points closer to the center of the x range. I believe that should be fine, though I doubt this was researched thoroughly. It has an adhoc stink to it. I wonder how many of these public TanhSinh implementations are numerically flawed.  Rob "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03292021, 02:28 AM
(This post was last modified: 03292021 05:00 PM by Albert Chan.)
Post: #4




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(03292021 01:10 AM)robve Wrote:(03282021 06:59 PM)Albert Chan Wrote: >>> from mpmath import * ∫(f(x), x=0..1) did not fail in Python mpmath either. The problem is uf = lambda u: 6*u*(1u) * f(u*u*(32*u)) ∫(uf(u), u=0..1) will touch f(x=1) = ∞, much earlier than expected. >>> utox = lambda u: u*u*(32*u) >>> utox(0.9999999962747) 0.99999999999999989 >>> utox(0.9999999962748) 1.0 Update: we can fix uf(u), to never evaluate endpoints. >>> def uf(u): # u = 0 .. 1 ... x = u*u*(32*u) ... if x==0 or x==1: return 0 ... return 6*u*(1u) * f(x) ... >>> quadts(uf, [0,1]) mpf('1.7724538508925098') Trivia: if u slightly outside valid range, x will still be inside [0, 1] utox(ε) = ε²*(32ε) ≥ 0, even if ε is negative utox(1ε) = (1ε)²*(32*(1ε)) = 1  ε²*(32ε) ≤ 1, even if ε is negative 

03292021, 07:42 PM
Post: #5




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
Albert, thanks for sharing additional examples and details.
Let me first say that I'm still "playing around" with this method using several implementations to compare (but I do not have much time beyond weekends.) TanhSinh also crashes in Boost Math library for \( \int_0^1 1/\sqrt{\log(1x)} \,dx \) with the Abort message: libc++abi.dylib: terminating with uncaught exception of type boost::wrapexcept<boost::math::evaluation_error>: Error in function tanh_sinh<double>::integrate: The tanh_sinh quadrature evaluated your function at a singular point and got inf. Please narrow the bounds of integration or check your function for singularities. Abort but does not crash for \( \int_0^1 1/\sqrt{\log x} \) which is the symmetric equivalent(!), i.e. x approaching 0 is more problematic than approaching 1. But the reason is probably different, since there is no explicit guard against evaluating at a+x or bx when x is too small (IMHO a practice that can have problems of its own, i.e. to ignore function values may produce inaccurate error estimates and less precise integrals). (03292021 02:28 AM)Albert Chan Wrote: The problem is uf = lambda u: 6*u*(1u) * f(u*u*(32*u)) This function works in Boost Math: Tanhsinh=0.5 err=3.66374e14 points=74 C implementation in this post: qthsh=0.5 k=4 err=1.55431e15 points=57 The Boost Math TanhSinh is a rather straightforward implementation. By comparison, the WP 43S code is more clever, but both are not always robust to functions that have problems in close proximity to the endpoints. The TanhSinh transformation does produce more points closer to the endpoints and performs particularly well for Ushaped functions or Lshaped on the domain of integration. Functions that are sensitive to errors in close proximity to the endpoints may not be usable with this method, depending on how many levels are (internally) used and if the implementation has guards that attempt to prevent this (if there are any at all). For the latter case, the "fix" in the code I looked at is likely a response to this problem (but still fails). (03292021 02:28 AM)Albert Chan Wrote: ∫(f(x), x=0..1) did not fail in Python mpmath either. A quick look at mpmath's quadrature internals reveals some interesting details. The implementation produces a list of (w,x) with weights w and nodes x to sum, see 260 class TanhSinh(QuadratureRule): The mpf precision is set extra high: 311 def calc_nodes(self, degree, prec, verbose=False): ... 329 extra = 20 330 ctx.prec += extra After producing the TanhSinh points in calc_nodes, the function to integrate is evaluated at all x in the node list, regardless how close x is to the endpoints: 308 S += self.ctx.fdot((w,f(x)) for (x,w) in nodes) where fdot ultiplies the weight w with f(x) represented in mpf, e.g. when w and x are mpf reals: 892 def fdot(ctx, A, B=None, conjugate=False): 944 real.append(mpf_mul(a._mpf_, b._mpf_)) ... 972 s = mpf_sum(real, prec, rnd) ... 976 s = ctx.make_mpf(s) 977 return s It looks like there is no error in this case due to higher precision. More luck than wisdom perhaps? Let me add that Boost Math TanhSinh does appear to offer a fix by passing the distance to the endpoint to the integrator: sqrt(x / (1  x * x)) This an example of an integral that has all its area close to a nonzero endpoint, the problem here is that the function being integrated returns "garbage" values for x very close to 1. We can easily fix this issue by passing a 2 argument functor to the integrator: the second argument gives the distance to the nearest endpoint, and we can use that information to return accurate values, and thus fix the integral calculation. As stated in the table with example cases spanning simple to problematic. And, of course: * As always, there are caveats: For instance, if the function you want to integrate is not holomorphic on the unit disk, * then the rapid convergence will be spoiled. In this case, a more appropriate quadrature is (say) Romberg, which does not * require the function to be holomorphic, only differentiable up to some order. "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03302021, 02:00 AM
Post: #6




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
Albert, in the TanhSinh programs that I based my C and BASIC implementation on, I found several issues that can be corrected and optimized. Love to get your feedback on those. I also "played around" a bit with the Python mpmath TanhSinh implementation, which is terrible (slower than GaussLegendre for several test cases when it should be faster and inaccurate results, perhaps more about that later...)
a. the TanhSinh code I looked at attempts to prevent evaluating the function at points a and b. However, as I've explained in my previous reply, these checks do not reliably work. I can think of at least two possible options to fix this, though it might also be fine to leave them out entirely, because it is known and understood that TanhSinh may fail for certain functions with singularities at the endpoints. Option 1: changing the inner loop termination condition from this: if (r == 0  r == 1) break; x = d*r; if (a+x > a) fp = f(a+x); if (bx < b) fm = f(bx); to the following lines of code that check if r<MachEps instead of r==0 and by removing the unnecessary adhoc guards that wrongly used the old fp and fm values (i.e. reused previous points), which I did not find any justification for and I seriously doubt is okay since it may affect the result and the error estimation: if (1+r == 1  r == 1) break; x = d*r; fp = f(a+x); fm = f(bx); This prevents x from getting too close to the endpoints a and b. This change takes the relative distance of the proximity of x to a and b into account, i.e. scaled with the size of the integration domain d=(ba)/2, note that we have x=d*r. Option 2: the first option can still fail if d<<1 and both a and b are large, i.e when ab<<1 and a,b>>1 (or when c>>1 center). So: if (d+r == d  r == 1) break; x = d*r; if (c+x == c) break; fp = f(a+x); fm = f(bx); These are two possible options to reduce the chances of hitting a and b too close. b. The code I looked at lacks some optimizations to further speed up the code, especially when exp() is expensive to compute, e.g. in software on BASIC pocket computers and calculators. We can apply two "loop strength reductions", which do not negatively affect precision. The first strength reduction I already applied to the BASIC program. It is a change from: double h; int k = 0; ... do { double p = 0, q, v; int j = 1; h = pow(2, k); // h=1,1/2,1/4,1/8,... do { double t = exp(h*j); ... j += 1+(k>0); } while (fabs(q) > fabs(p*eps)); to the strength reduced form to remove pow(2, k): double h = 2; int k = 0; ... do { double p = 0, q, v; h /= 2; do { double t = exp(h*j); ... j += 1+(k>0); } while (fabs(q) > fabs(p*eps)); ... } while (e >= tol && k <= n); ... return d*s*h; // we need the final h here Then we apply a second loop strength reduction by observing that exp(h*j)=exp(h)^j then strengthreduce the power ^j away to obtain: double h = 2; int k = 0; ... do { double p = 0, q, v; double exph, exphj; h /= 2; exph = exp(h); exphj = exph; if (k > 0) exph *= exph do { double t = exphj; ... exphj *= exph; } while (fabs(q) > fabs(p*eps)); ... } while (e >= tol && k <= n); ... return d*s*h; // we need the final h here Note that exp(h) in the outer loop has up to k==6 or k==7 values, namely \( \{e, e^{1/2}, e^{1/4}, e^{1/8}, e^{1/16}, e^{1/32}, e^{1/64} \} \) and can be precomputed to speed this up even more. I experimented with these changes to verify them and everything looks good. More testing is necessary and additional improvements might be possible. More to come...  Rob "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

04012021, 04:27 PM
(This post was last modified: 04022021 12:30 AM by robve.)
Post: #7




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
If anyone is interested, I uploaded a draft article here that observes and addresses the following issues in existing (said to be "the best") implementation of TanhSinh:
 fix division by zero in e=fabs(v/s1), which is not safe to execute, e.g. integrating f(x)=x.5 causes floating point exception;  fix underflow in the inner loop condition fabs(q)>fabs(p*eps), e.g. integrating f(x)=sin(pi*x*40) evaluates too many (315) points (versus 15);  by comparison, Python mpmath TanhSinh implementation's error estimation is not as reliable as can be;  performance improvements to speed up qthsh TanhSinh in C and BASIC, runs at least 10% faster on a pocket computer;  permit ±inf singularities close to one or both of the endpoints. This is work in progress and things may change, e.g. the number of functions tested so far is limited to 21 but those reveal some interesting pro/cons successes and failures of the method. I hope this is useful to anyone. Any suggestions and comments are much appreciated.  Rob "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

04022021, 11:05 AM
(This post was last modified: 04022021 12:16 PM by emece67.)
Post: #8




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
Hi,
A question, the integration interval in your calculations, is always [0, 1]? Regards. EDIT: I answer myself: yes. 

04022021, 02:46 PM
(This post was last modified: 04022021 07:35 PM by robve.)
Post: #9




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(04022021 11:05 AM)emece67 Wrote: A question, the integration interval in your calculations, is always [0, 1]? Yes. These are all definite TanhSinh integrals. I have not looked at unbounded domains, for which SinhSinh and ExpSinh are applicable. Functions that are typically integrated over non[0,1] finite domains are normalized to the [0,1] domain in experiments. It is simpler that way to conduct the tests and for comparisons. By implication, this work also applies to non[0,1] finite domains. I've updated the draft article today with the latest improvements and corrections to TanhSinh. Once this is done I will update the original post with the list of improvements and the final code in a couple of days.  Rob Edit: needed clarification: not updating the original post and code until the work is done in a couple of days "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

04022021, 03:16 PM
Post: #10




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
https://www.genivia.com/files/qthsh.pdf, page 10 Wrote:mpmath: mp.dps=15, points=53, est.err=1e14 (incorrect, actual is about 1e9) I was unable to produce above result, for ∫(x**3  2x*x + x, x = 0 .. 1) >>> def f(x): f.n += 1; return x*x*x2*x*x+x ... >>> f.n = 0 >>> quadts(f, [0,1], error=True), f.n ((mpf('0.083333333333333329'), mpf('1.0e22')), 53) Integral look too simple to trip TanhSinh quadrature. With GaussLegendre quadrature, it took even less points. >>> f.n = 0 >>> quadgl(f, [0,1], error=True), f.n ((mpf('0.083333333333333329'), mpf('9.26442286059391e23')), 9) Both algorithms return correct quadrature, 1/4  2/3 + 1/2 = 1/12 

04022021, 03:40 PM
Post: #11




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(04022021 03:16 PM)Albert Chan Wrote:https://www.genivia.com/files/qthsh.pdf, page 10 Wrote:mpmath: mp.dps=15, points=53, est.err=1e14 (incorrect, actual is about 1e9) See the number in the table in the article: (53,1E14#) which is about the same what you found. I neglected to update the title with the correct number. The relative error is about 1E17~1E16 to be fair. The mpmath error estimation is a bit of a mystery to me, given the numbers and the way it is estimated in the code. It seems too high or too low in several cases.  Rob "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

04022021, 04:20 PM
Post: #12




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(04022021 02:46 PM)robve Wrote: I've updated the draft article today with the latest improvements and corrections to TanhSinh. Once this is done I will update the original post with the list of improvements and the final code. I'll also post later my original Python code (from where I later derived the wp34s code). Then I used a battery of 140 integrals (now 161, after adding your 21 functions) for tests. Maybe this set is useful to you. Regards. 

04022021, 05:04 PM
Post: #13




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(04022021 04:20 PM)emece67 Wrote: Then I used a battery of 140 integrals (now 161, after adding your 21 functions) for tests. Maybe this set is useful to you. Thanks for contributing. That would be nice to look at. Eventually I may run the 1200 integrals in the VB spreadsheet. However, it is not about the number of functions per se (though the more the better), the variety matters too, i.e. include nonholomorphic functions, periodic, step functions, etc. Testing these can help to determine if the error estimate is usable, not only accurate. TanhSinh is also not "symmetric", meaning that integrating 1/sqrt(log(x)) over [0,1] should be the same as integrating 1/sqrt(log(1x)). Except for mpmath, that is not the case with the VB article code and Boost Math due to the singularities occurring earlier that apparent to be fatal. I've used the default mpmath TanhSinh parameters with mp.dps=15 in order to compare to C/C++ double fp. Levels is 6 (default) and eps is 1e9. Some observations:  mpmath TanhSinh appears to run up to 427 evaluations, never more with these settings and often too many than necessary. Something might be amiss with the convergence test?  mpmath TanhSinh sometimes produces error estimates that are a mystery when compared to the other TanhSinh implementations. I could be mistaken, but the error estimate is performed in the same way as GausssLegendre in the code. Both methods produce node sets which are used for the quadrature and error estimation.  Boost Math TanhSinh explodes the number of points close to the zero endpoint (that's also one reason to include the zero endpoint in the tests to compare behaviors, not only the integral and error estimates.) For example, f(x)=4/(1+x*x) produces 147 points with a large number close to the 0 and 1 endpoints. See my draft article's graph on this function with Boost Math, and a lot that are very close to the 0 endpoint: 1.69367010239342e108 7.55408323550392e123 4.13900004783521e139 1.55080149097904e157 2.04586734643622e178 4.47908407247427e202 6.93948438714987e229 2.89628710646752e259 Of course, we expect points to be distributed towards the endpoints with TanhSinh, but this is far more than necessary. The weights there are minuscule and negligible. PS. Just to clarify: I am not really biased toward any integration method. They all have strengths and weaknesses that can favor one over the other depending on what the goal is (e.g. CPU time versus accuracy). The best integration results are obtained by applying a transformation when necessary, change variable when helpful, consider integration domain splitting, and then select the method to apply that works best for the integrand. A usable error estimate is just as important (if not more) than the integration result, unless we seek a rough estimate (e.g. 1e3) that is good enough.  Rob "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

04022021, 07:12 PM
Post: #14




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(03292021 07:42 PM)robve Wrote: sqrt(x / (1  x * x)) Can you explain what that meant ? Perhaps an example ? To avoid garbage of x close to 1, I would have rewrite the integrand, with infinity at x=0. (y=1x², dy = 2x dx, would solved it exactly. But, no cheating ...) ∫(x/√(1x²), x=0..1) = ∫((1x)/√(x*(2x)), x=0..1) With infinity moved, at x=0, we can add many, many more points. Since e^x grow faster than x^n, we have e^x shrink faster then x^n, for finite n. ∫(f(x), x=0 .. 1) = ∫f(e^y)*(e^y), y=0 .. inf) Code: def e_transform(f, exp=exp): # interval [0,1] > [0,inf] >>> f = lambda x: (1x)/sqrt(x*(2x)) >>> ef = e_transform(f) >>> quadts(f, [0, 1]) 0.999999999624821 >>> quadts(ef, [0, inf]) 1.0 Similar results, with GaussLegendre quadrature: >>> quadgl(f, [0, 1]) 0.993619917982059 >>> quadgl(ef, [0, inf]) 1.0 

04022021, 08:14 PM
Post: #15




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(04022021 07:12 PM)Albert Chan Wrote:(03292021 07:42 PM)robve Wrote: sqrt(x / (1  x * x)) What is meant is that when x is close to 1 socalled "catastrophic cancellation" occurs with loss of precision. TanhSinh accumulates inaccurate values at this point. I picked two Boost Math TanhSinh points very close to 1 for this function to illustrate: x=0.999999999999991 x^2=0.999999999999982 1x^2=1.77635683940025e14 x=0.999999999999999 x^2=0.999999999999999 1x^2=1.33226762955019e15 The first point for example has about 3~4 digits remaining in IEEE 754 double precision after subtraction 1x^{2} with x^{2}= 0.999999999999982, so roughly 1.77e14 < 1x^2 < 1.78e14. This means that x/(1x*x)) is accurate to 3~4 digits while being fairly large 5.6e13. This happens to a lot of points. The "garbage" is not my wording, see Boost TanhSinh. There could be a problem if the number of points integrated are mostly in the neighborhood of 1 for this function. Thinking that there are better examples to demonstrate this.  Rob "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

04022021, 09:07 PM
Post: #16




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
I also noticed other issues in the VB code in the Article, which I have not mentioned in my writeup, such as computing the relative error with:
errval = Abs(2 * sslast / ss  1) However, the distributive law does not apply to floating point arithmetic. This should be: errval = Abs((2 * sslast  ss) / ss) When 2*sslast ~ ss this is important. However, this is a bit pointless to say, because both pieces of code may cause division by zero when the updated quadrature sum ss is zero. A function such as f(x)=x.5 over [0,1] fails to integrate. It can get worse, because a division by zero produces errval=NaN. NaN is then compared to a threshold value. Comparisons to NaN are always false, even the test x==x is false, which is a neat way to detect that x==NaN by the way. NaN comparisons in loop conditions can lead to early termination (probably good) or nontermination (probably bad). In the VB code there is a WHILE instead of an UNTIL so NaN is "good" instead of "bad": Loop While errval >= tol And k <= maxlevel These kind of issues happen more often than we realize (so I'm hopefully not embarrassing anyone to point this out, because it happens a lot). Because FP is generally hard to use, a while ago I wrote a couple of lecture notes on the subject for graduate students. It's the usual fare on FP and SIMD and there is much more that I would like to add, in particular the reading material "what every computer scientist should know about floating point arithmetic". If you're in to this kind of low level detail and there is not much to do because of the pandemic, then perhaps this is a text to get excited about!  Rob "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

04032021, 01:40 AM
(This post was last modified: 04032021 01:14 PM by Albert Chan.)
Post: #17




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(04022021 02:46 PM)robve Wrote: Functions that are typically integrated over non[0,1] finite domains are normalized to the [0,1] domain in experiments. Here is a link to convert integral limits to [0, 1]: http://fmnt.info/blog/20180818_infiniteintegrals.html Example, my previous post etransform, without changing integral limits. Code: def e_transform2(f, exp=exp): # interval [0,1] > [0,1] >>> f = lambda x: x**0.99 >>> ef2 = e_transform2(f) >>> quad(f, [0,1]) mpf('35.65470993432163') >>> quad(ef2, [0,1]) mpf('99.999999999999915') Comment: although convenient to have unit interval, integral may be harder to compute. Conversion of [0,∞] > [0,1], compress the curve, possibly making it spiky. Cas> f(x) := x^0.99 Cas> ef1(t) := f(e^t)*e^t; // [0,1] > [0,inf] Cas> ef2(t) := ef1(1/t1)/(t*t); // [0,inf] > [0,1] ef1(t) = exp(0.01*t), no spike, thus easy to integrate. ef2(t) has a spike at t = 0.005 Cas> ef2(0.007) → 4939.99121224 Cas> ef2(0.005) → 5467.81701782 Cas> ef2(0.003) → 4003.61366011 Cas> ef2(0.001) → undef Cas> limit(e2(t), t=exact(0.001)) → 1000000*exp(98901/100)/exp(999) Numerically, e2(0.001) will be hit with overflow, ∞/∞ → undef (*) (*) mpmath does not overflow, nor underflow >>> 1000000*exp(98901/100)/exp(999) mpf('45.856206642206899') Thus, mpmath can easily plot the shape of spike (Cas cannot) >>> plot(ef2, [0, 0.01]) 

04032021, 03:17 AM
(This post was last modified: 04032021 01:13 PM by emece67.)
Post: #18




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
Hi,
These are the results for my Python script: Code: TNFE err_r err_t where err_r is the reported error and err_t is the true error. TNFE is the total number of function evaluations. The only adjustable parameter is mp.dps (set to 15, seems to work, at least, with values 8~128). It worked as (I) expected, except for functions 10~11... they hurt! The script is suffering now from some refactoring, so I expect to upload it to GitHub in a few days. Regards. EDIT: results modified as a result of code refactoring. 

04032021, 01:38 PM
Post: #19




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
My original Python code for the doubleexponential quadrature is now on GitHub.
It has suffered some refactoring from its original shape (that I used to derive the wp34s code). Regards. 

04042021, 01:50 AM
Post: #20




RE: (PC12xx~14xx) qthsh TanhSinh quadrature
(04032021 01:40 AM)Albert Chan Wrote:(04022021 02:46 PM)robve Wrote: Functions that are typically integrated over non[0,1] finite domains are normalized to the [0,1] domain in experiments. Very nice to know. Just for clarification: we're only considering definite integrals at this time with TanhSinh. A simple linear transformation suffices to rescale the bounds. The same area is computed with the same function points, but via a transformation. No tricks, just easier to analyze the point distributions in a normalized [0,1] domain. For example, NR2ed p.134 is transformed: $$ \int_0^2 x^4 \log(x+\sqrt{x^2+1})\,dx = \int_0^1 32u^4 \log(2u+\sqrt{4u^2+1})\,du $$ I've tested other integrals including the one above that are not included in the article, but several of these are not very interesting for TanhSinh, e.g. trig and other transcendental functions. We already know that TanhSinh converges with double exponential rate for holomorphic complexdifferentiable functions bounded on the interior of the unit disk within the Hardy space. The point distribution makes this method more resilient to singularities at the endpoints and is especially good when most of the area is located to the endpoint(s), e.g. L and Ushaped functions. As some have said, this method might never reach the level of a true generalpurpose integrator, but who knows... It's interesting to analyze the implementation differences, apply code optimizations and fixes when applicable to public code. Perhaps we can find a way to make the method stronger when exceptions occur due to singularities that kill the method (thanks Albert for pointing this out!). The approach that might be best to do this is to interpolate/extrapolate, see the last part of the (draft) article: "The integral of function 15 is now the same as the integral of function 14 as it should be, up to the last digit 1.77245382409889. Both are the same for the same specified eps." Though anecdotal at this early stage, this corroborates the interpolation approach. Inverting the integration bounds should give the same integration result, which for TanhSinh is NOT necessarily the case for domains such as [0,1]. I did not have time today to look at emece67 Python code and results yet. As always, suggestions and comments are much appreciated.  Rob "I can count on my friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)