(PC-12xx~14xx) qthsh Tanh-Sinh quadrature
|
04-08-2021, 01:58 PM
Post: #39
|
|||
|
|||
RE: (PC-12xx~14xx) qthsh Tanh-Sinh quadrature
(04-08-2021 11:20 AM)emece67 Wrote: integrals (1 - sqrt(2)*cosh(x)/sqrt(cosh(2*x)), from -inf to +inf) needed one more level. I believe there is nothing to worry about. Integrand suffered from subtraction cancellations. A slight shift to the nodes/weights (even if more precise) does not translate to better result. Even for 4 good digits, consider yourself lucky. mpmath quadts only get 3. >>> f = lambda x: 1-sqrt(2)*cosh(x)/sqrt(cosh(2*x)) >>> quadts(f, [-inf, inf], error=True) (mpf('-0.69306562029329821'), mpf('0.001')) Actually, mpmath "cheated", calcuating in higher precision. See what happen if we pre-calculated sqrt(2), under mp.dps=15 setting. >>> f2 = lambda x, k=sqrt(2): 1-k*cosh(x)/sqrt(cosh(2*x)) >>> quadts(f2, [-inf, inf], error=True) # failed convergence (mpf('-1886.3975082035695'), mpf('1.0')) This is a weakness of tanh-sinh quadrature. It relied on accurate sample points. Gauss-Legendre quadrature, using curve-fitting of samples, handle "fuzzy" more easily. >>> quadgl(f2, [-inf, inf], error=True) (mpf('-0.69314718056249158'), mpf('1.0e-14')) We could rewrite integrand, to perform better. Let c = cosh(x), d = cosh(2x) With half-angle formula, we have c*c = (1+d)/2 1 - √(2)*c/√(d) = 1 - √(1+1/d) = -(1/d) / (1 + sqrt(1+1/d)) >>> def f3(x): y=1/cosh(2*x); return -y/(1 + sqrt(1+y)) ... >>> quadts(f3, [-inf, inf], error=1) (mpf('-0.69314718055994529'), mpf('1.0e-29')) (04-06-2021 11:36 PM)robve Wrote: When reusing previous points is better as default for f on [a,b]: Same issue with all of these. A simple rewrite "sharpen" integrand greatly. First one had the form 1/(1+u). With identity 1/(1+u) + 1/(1+1/u) = 1/(1+u) + u/(1+u) = 1 ∫(1/(1+cos(x)^x), x=-pi/2 .. pi/2) = ∫(1, x=0 .. pi/2) = pi/2 With identity sin(x/2)^2 = (1-cos(x))/2, we can rewrite the others: x^2/(1-cos(x)) = x^2 / (2*sin(x/2)^2) = (x/sin(x/2))^2 / 2 x*sin(x)/(1-cos(x)) = (x*2*sin(x/2)*cos(x/2) / (2*sin(x/2)^2)) = x / tan(x/2) |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)