∫Adaptive(f,a,b,N,ε) ... Carlos

01422d4c-be8a-41ef-9a72-b799cb045de1 15 4 5 decimal

&[DATE] &[TIME] - &[FILENAME]

x f x - exp x sqrt / :

x g 1 x - ln - x / :

ε 1015 -^:

x λ 1 ε x ≤ (x 1 ≤ (& cases :

Singular end point [blue] ... infinite limit [red]

x g x λ * x f x λ * 1 10 1 1 - 2 2 mat 3 1 sys

Two useless Simpson [10000]

x f x ε 1 int 52706

x g x ε 1 ε - int 1.6894

∫Adaptive(f,a,b,N,ε) ... Carlos

f a b N ε ∫Adaptive s 0 : sum 0 : i 1 : tol i el ε : aa i el a : h i el b a - (2 / : 51 line .... fa i el a f : fc i el a h i el + f : fb i el b f : s i el h i el fa i el 4 fc i el * + fb i el + (3 / * : L i el 1 : 51 line 13 mat above => initialise the while loop i 0 > fd aa i el h i el 2 / + f : fe aa i el 3 h i el * 2 / + f : 12 mat s1 h i el fa i el 4 fd * + fc i el + 6 / (* : s2 h i el fc i el 4 fe * + fb i el + 6 / (* : 12 mat v1 aa i el : v2 fa i el : v3 fc i el : v4 fb i el : v5 h i el : v6 ε : v7 s i el : v8 L i el : 18 mat i i 1 - : s1 s2 + v7 - abs v6 < sum sum s1 s2 + (+ : v8 N ≥ Solution Insufficient Number of Intervals : i i 1 + : aa i el v1 v5 + : fa i el v3 : fc i el fe : fb i el v4 : 14 mat h i el v5 2 / : tol i el v6 2 / : s i el s2 : L i el v8 1 + : 14 mat i i 1 + : aa i el v1 : fa i el v2 : fc i el fd : fb i el v3 : h i el h i 1 - el : s i el s1 : 16 mat L i el L i 1 - el : 61 line if if 51 line while sum 51 line eval :

Test function x f 1 x - ln - x / eval :

'N' iterations
precision N 100 : acc 1015 -^: 21 sys

User inputs

lower limit
upper limit a 1015 -^: b 11015 -^- : 21 sys

t0 1 time :

Zeta(2) ... f a b N acc ∫Adaptive 1.644934066848300

Zeta(2) 1.644934066848226

1 time t0 - 32 s *

x - exp x sqrt /

SS 6179 ......... fails miserably
Adaptive MCD .... 1.49364826562485
Wolfram α  ....... 1.49364826562485
sqrt(π)*erf(1)... 1.49364826562485
Simplex Smath ... 1.4936482656238
∫kernel Smath ... 1.49362 else

'N' iterations
precision x f x - exp x sqrt / eval : N 100 : acc 1015 -^: a 1027 -^: b 1 : 51 sys

t0 1 time :

Zeta(2) ... f a b N acc ∫Adaptive 1.49364826562483

1 time t0 - min 1.2

1D simplex integration

∫simplex(f,a,b,s)

k gencm1 N 2 k matrix : i 1 k range N 1 i el i 1 - : N 2 i el k i - : 21 line for N 31 line :

nr nc v genm M nr nc matrix : i 1 nr range j 1 nc range M i j el v : for for M 31 line :

5 gencm1 012344321025 mat

330 genm 00000000033 mat

Carlos added,
enables recursive N 2 k matrix ≡ M nr nc matrix ≡ 21 sys

M idx M transpose 1 col rows :

5 gencm1 idx 5

Carlos  'submat' M 1 n 1 + (range 2 w range el M 1 n 1 + (2 w submatrix ≡

330 genm idx 3

without submatrix(...) n k gencm n 1 ≡ k gencm1 M n 1 + 10 genm : i 0 k range N n 1 - k 1 + (i - gencm : w N idx 1 - : O 1 w 1 + i genm : N O N stack : M M N augment : 51 line for w M idx 0 - : M 1 n 1 + (range 2 w range el 41 line if 11 line :

15 gencm 012344321025 mat

Only true for n=2 and any 'k'

24 gencm 11 mat 000001111222334012340123012010432103210210100315 mat 11 mat

weights(n,s)is universal for any dimension n s weights d 2 s * 1 + : i 0 s range w i 1 + el 1 - (i^22 s * -^* d n + 2 i * - (d^i! (d n + i - (! (* / * : for w 31 line :

Simplex integration algorithm f n s ∫ v 0 : Pt n s pts : w n s 1 - weights : i 1 w rows range P Pt i el : v v w i el P j col f j 1 P idx sum * + : 21 line for v 51 line :

n s pts d 2 s * 1 + : i 0 s range M n s 1 + (i - gencm : M M 2 n 1 + 1 M idx submatrix : Pt i el 2 M * 1 + d n + 2 i * - / : 31 line for Pt 31 line :

Specific for 1D

Recast the simplex domain

Simplex 1D integrator  ... domain [a,b] ... 's' convergence f a b s ∫simplex u b a - (34 / 2 x 1 el * 1 - (* 14 / 2 x 1 el * 1 - (3^* - 12 / + (* a + : x f u f 12 x 1 el * 1 - (2^- (* : simplex 32 / b a - (* f 1 s ∫ * : simplex 41 line :

t0 1 time :

sanity 12 D f 0111 ∫simplex 1.493648265627

π sqrt 1 erf * 1.49364826562485

1 time t0 - 0.3 s *

Explore i 714 range ∫ i 7 - el f 01 i ∫simplex : for ∫ 21 line :

Explore 1.49364827281497 1.49364826611792 1.49364826565769 1.49364826562703 1.49364826563011 1.49364826563789 1.49364826559989 71 mat

Simplex reverses convergence as 's' increases. Once you have determined the direction wrt f(x) use the reported best 's' ... here 's=11'

Search: 's' reverse convergence s k 6 : s k 15 range f 01 s ∫simplex f 01 s 1 - ∫simplex ≥ index k el s : k k 1 + : 21 line k if for s k 1 + : 31 line :

Optimized best simplex 's' s 11

sanity ...  1.49364826562485 ... 12 D f 0111 ∫simplex 1.49364826562703

π sqrt 1 erf * f 0111 ∫simplex - 21012 -^* -

simplex x - exp x sqrt /

adaptive 1 x - ln - x /

Not so sure ... sanity check

∫Adaptive(f,a,b,N,ε)

f a b N ε ∫Adaptive s 0 : sum 0 : i 1 : tol i el ε : aa i el a : h i el b a - (2 / : 51 line .... fa i el a f : fc i el a h i el + f : fb i el b f : s i el h i el fa i el 4 fc i el * + fb i el + (3 / * : L i el 1 : 51 line 13 mat above => initialise the while loop i 0 > fd aa i el h i el 2 / + f : fe aa i el 3 h i el * 2 / + f : 12 mat s1 h i el fa i el 4 fd * + fc i el + 6 / (* : s2 h i el fc i el 4 fe * + fb i el + 6 / (* : 12 mat v1 aa i el : v2 fa i el : v3 fc i el : v4 fb i el : v5 h i el : v6 ε : v7 s i el : v8 L i el : 18 mat i i 1 - : s1 s2 + v7 - abs v6 < sum sum s1 s2 + (+ : v8 N ≥ Solution Insufficient Number of Intervals : i i 1 + : aa i el v1 v5 + : fa i el v3 : fc i el fe : fb i el v4 : 14 mat h i el v5 2 / : tol i el v6 2 / : s i el s2 : L i el v8 1 + : 14 mat i i 1 + : aa i el v1 : fa i el v2 : fc i el fd : fb i el v3 : h i el h i 1 - el : s i el s1 : 16 mat L i el L i 1 - el : 61 line if if 51 line while sum 51 line eval :

β 0.125 :

x f 11 x abs β^+ / :

The mapping function

x f 0.275 0.625 0.275 0.25 2 2 mat 2 1 sys

acc 100 ... NFG x f x 1 - 1 int 1.063963

'N' iterations
precision N 100 : acc 1015 -^: 21 sys

lower limit
upper limit a 1 - : b 1 : 21 sys

t0 1 time :

Zeta(2) ... f a b N acc ∫Adaptive 1.06202606342187

1 time t0 - 36 s *

MCD 11 1.06202606342187

lo 1 - : hi 1 : N 19999 : 13 mat

β 0.125 :

x f 11 x abs β^+ / :

data x lo hi lo hi lo - (N / + range : i 1 x rows range f i el x i el f eval : for x f augment 31 line :

data Σ Finite differences Cumul Integration vx data 1 col : vy data 2 col : Y 1 el 0 : i 2 vx rows range Δx i el vx i el vx i 1 - el - : Y i el Y i 1 - el 12 / vy i 1 - el vy i el + (Δx i el * (* + : 21 line for vx Y augment 61 line :

n N 1 + : 20000

Σ data Σ eval :

NOT deceptive from abundant 'N' Σ n row 1 1.062025 12 mat

MCD 11 1.06202606342187

∫Adaptive(f,a,b,N,ε) Carlos

f a b N ε ∫Adaptive sum 0 : i 1 : tol i el ε : aa i el a : h i el b a - (2 / : 51 line .... fa i el a f : fc i el a h i el + f : fb i el b f : s i el h i el fa i el 4 fc i el * + fb i el + (3 / * : L i el 1 : 51 line 13 mat above => initialise the while loop i 0 > fd aa i el h i el 2 / + f : fe aa i el 3 h i el * 2 / + f : 12 mat s1 h i el fa i el 4 fd * + fc i el + 6 / (* : s2 h i el fc i el 4 fe * + fb i el + 6 / (* : 12 mat v1 aa i el : v2 fa i el : v3 fc i el : v4 fb i el : v5 h i el : v6 ε : v7 s i el : v8 L i el : 18 mat i i 1 - : s1 s2 + v7 - abs v6 < sum sum s1 s2 + (+ : v8 N ≥ Solution Insufficient Number of Intervals : i i 1 + : aa i el v1 v5 + : fa i el v3 : fc i el fe : fb i el v4 : 14 mat h i el v5 2 / : tol i el v6 2 / : s i el s2 : L i el v8 1 + : 14 mat i i 1 + : aa i el v1 : fa i el v2 : fc i el fd : fb i el v3 : h i el h i 1 - el : s i el s1 : 16 mat L i el L i 1 - el : 61 line if if 51 line while sum 41 line eval :

∫kernel(f,a,b,n) ... Simpson/Romberg Kernel [Jean]

Basic Romberg-Simpson ... integration acc from options ... NUMERIC f a b Romberg w s b a - (34 / 2 w * 1 - (* 14 / 2 w * 1 - (3^* - 12 / + (* a + : w G w s f 12 w * 1 - (2^- (* : 32 / b a - (* w G w 01 int * 31 line eval :

f a b n ∫kernel

Simpson Integrator in the Romberg Transform Domain [0,1] <= [a,b] ...  
... integrates  Smath ainterp(X,Y,x), cinterp(X,Y,x), cubic(x) from  l_p_c_spline, Y^X 
Integration accuracy from user command line 'n' ... NUMERIC.
The integration variate is in Romberg [0,1] within 'n' ... 'n' independent of work sheet option. 
Thus, on long range of the integration variate ∫kernel(f,a,b,n) improves the convergence. f a b n ∫kernel w s b a - (34 / 2 w * 1 - (* 14 / 2 w * 1 - (3^* - 12 / + (* a + : w G w s f 12 w * 1 - (2^- (* : ... α 0 : β 1 : h β α - n / : 15 mat Simpson 13 / h 2 / α G β G + (* h α k h * + G k 1 n 1 - sum * + (* 23 / h * α h 2 / - k h * + G k 1 n sum * + eval : 32 / b a - (* Simpson * 41 line :

a x Φ a 2^a x - (a x + (* / :

Part 1: Attempt the symbolic integral curve

a x f a x Φ x int maple : a a x + a x - / ln * 2 /

ε 1012 -^:

x λ 1 0.75 - x ≤ (x 0.75 ≤ (& cases :

a 1 :

a x f x λ * a x f a x Φ 3 1 sys

Part 2: Difficult from insight ........double Singular end points

x integrand a 2^a x - (a x + (* / :

ε 1014 -^:

x f x integrand : iter 1000 : mesh 0.1 : x1 0 : x2 1 ε - : 51 sys

t0 1 time :

Σ f x1 x2 iter mesh ∫Adaptive :

closely convergent Σ 16.278

1 time t0 - 1 s *

x canvas f 0 x iter 0.1 ∫Adaptive :

f x1 x2 1000 ∫kernel Does not converge ≡

t0 1 time :

size canvas < 1 ... otherwise NO plot

x f x canvas 2 1 sys

1 time t0 - min 5.4

x integral a a x + a x - / ln * 2 / :

x2 integral x2 canvas x2 - canvas 31 mat 16.465 16.278 16.278 - 31 mat

Explain Part 2 . . .

1. If you can't find symbolic integral, transit via discrete/spline.

2. 10^-14 is typical safe ln(x) branch margin.

3. X_Y rejects canvas(x).

4. If QuickPlot canvas > 1 . . . takes 15 min to plot

]]>