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

590f2a2a-dd9c-4ef5-a5af-2a7884296e14 30 4 5 decimal

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

Simpson-Romberg integration . . .

Exmples [1, 2, 3] were chosen for their known value.

Upper limit [b] may be discretized for full graph. This document

reconciliates the integration operator wrt acc, as command line.

]]>

3 appVersion 0.98.6179 ≡ Timing 3/22 min ≡ 21 mat

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 f a b n Int 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 / : 14 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 :

sanity check the sextuplets Mathcad 11 Adaptiv 5.40752818465632 ≡ Technical accuracy 5.40752818465627 ≡ Simpson Transf'd.. 5.40752818465635 ≡ Smath Simpson 7½ D 5.4075281 ≡ Simplex 10½ D .... 5.40752818462 ≡ 51 line

x f x sqrt exp :

t0 1 time :

transformed domain [0,1] <= Simpson [a,b] f 022000 Int 5.40752818465635

1 time t0 - 9.6 s *

x c 11015 -^x ≤ (x 2 ≤ (& cases :

TOL 1014 -^:

Zeta(2) => 1.644934066848226

1 x - ln - x / x c * x 0.5 1 - ^ x - exp * x c * x sqrt exp x c * 1 0 + 200 black 1 5 mat 4 1 sys

singular end point x f 1 x - ln - x / :

t0 1 time :

Zeta(2) π²/6 ... 1.644934066848226
∫Adaptive(f,TOL,(1-TOL),100000,ε) f TOL 1 TOL - (50000 Int 1.64493406675557

1 time t0 - min 0

infinite limit x f x 0.5 1 -^x - exp * :

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

π sqrt 1 erf * 1.49364826562485

t0 1 time :

'n 100000' => 12 min ...  0.000006 f a b 100000 Int 1.4936425

1 time t0 - min 0

Adaptive integration . . .

Exmples [1, 2, 3] were chosen for their known value.

Upper limit [b] may be discretized for full graph.

]]>

∫Adaptive(f,a,b,N,ε) Program algo style is a bit slower than open style, more convenient. Adapted from Carlos code [20180110]

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

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 :

3 appVersion 0.98.6179 ≡ Timing 3 min ≡ 21 mat

Example_1: Ordinary difficult region

x c 1 0.001 x ≤ (x 2 ≤ (& cases :

Ordinary x f x sqrt exp :

'N' iterations
precision N 1000 : ε 1015 -^: 21 sys

User inputs

lower limit
upper limit a 0 : b 2 : 21 sys

t0 1 time :

Technical acccuracy .....   5.40752818465627 f a b N ε ∫Adaptive 5.40752818465644

1 time t0 - 24 s *

Example_2: Singular end point

x f 1 x - ln - x / :

TOL 1013 -^:

TOL 1013 -^≤

Zeta(2) => 1.644934066848226

x f x c * 1 10 1 1 - 2 2 mat x sqrt exp x c * 3 1 sys

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

t0 1 time :

True value Zeta(2)  [π²/6].......... 1.644934066848226 f a b 200 TOL ∫Adaptive 1.64493406684641

1 time t0 - 22 s *

Example_3: infinite limit

infinite limit x f x 0.5 1 -^x - exp * :

'N' iterations
_ precision _iter 1000 : TOL 1015 -^: 21 sys

x f

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

t0 1 time :

sanity Wolfram ... 1.49364826562485 f a b iter TOL ∫Adaptive 1.49364626562488

1 time t0 - min 1.8

π sqrt 1 erf * 1.49364826562485

an easy test ...

infinite limit x f 1 x / :

'N' iterations
_ precision _iter 1000 : ε 1015 -^: 21 sys

lower limit [any]
_ upper limit _a 104 -^: b 1 : 21 sys

t0 1 time :

sanity 7 D [MCD adaptive] f a b iter ε ∫Adaptive 9.2103404

1 time t0 - min 1.4

Simpson [acc 1000] ! 1 x / x 0.0001 1 int 9.7426574 ≡

Γ(x) [15 D] Alvaro

Γ.κ 0.99999999999999709182 57.156235665862923517 59.597960355475491248 - 14.136097974741747174 0.49191381609762019978 - 0.33994649984811888699 104 -^* 0.46523628927048575665 104 -^* 0.98374475304879564677 104 -^* - 0.15808870322491248884 103 -^* 0.21026444172410488319 103 -^* - 0.21743961811521264320 103 -^* 0.16431810653676389022 103 -^* - 0.84418223983852743293 104 -^* 0.26190838401581408670 104 -^* - 0.36899182659531622704 105 -^* 151 mat :

ζ Γ z ζ 1 - : zh z 0.5 + : zg zh 607128 / + : zp zg 0.5 zh *^: 14 mat 2 π * sqrt Γ.κ 1 el Γ.κ κ 1 + el z κ + / κ 114 sum + (* zp 2^* e zg -^* eval 21 line :

Valid for Re(ζ)>0

infinite limit x f x Γ :

x Γ

'N' iterations
_ precision _iter 1000 : ε 1015 -^: 21 sys

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

t0 1 time :

[1000] sanity  8.9637038 [MCD adaptive] f a b iter ε ∫Adaptive 8.9637042

[10000], 3 min x f x a b int 8.965116 ≡

1 time t0 - min 0

hatching utility

Erzeugt eine Matrix für die Darstellung vertikaler Linien zwischen den Funktionen f(x) und 
g(x) im Bereich von x1 bis x2 mit 2*n Linien. Generates a matrix for 2*N vertical lines (hatching) 
between functions f(x) and g(x) in x-range from x1 to x2 f g x.1 x.2 N hatch dx x.2 x.1 - N / (: h x.1 x.1 f x.1 x.1 g 22 mat : x x.1 x.2 dx - x.1 dx + range h h x dx 2 / + x dx 2 / + g x dx 2 / + x dx 2 / + f x dx + x dx + f x dx + x dx + g 42 mat stack : for h 41 line :

x f 1 x / :

x g 0 :

b 3 :

What is ln(x) ?

The quadratic hyperbola 1/x that generates ln(x)

x f x g x f x g 1 b 10 hatch x f x g 0.1 1 10 hatch 1 0 + 500 black 1 5 mat x ln 6 1 sys

In short math words, ln(x)is the area under 1/x between the reference limits 1 ... ± b. In Smath and any CAS, ln(x) is fetched from Windows at basic 21 Decimals floating point accuracy, of which by convention 15 are displayed.

x f x 1 b int b ln 21 mat 1.098612 1.098612 21 mat

z Ln 1 u / u 12 int :

2 Ln 2 ln 21 mat 0.693147180579478 0.693147180559945 21 mat

sanity Wolfram [10000]
14.0258509299405 x ln x 110 int 14.0258508845545

u b 0 u 2 - ≤ 16 / 2 u + (3^* 2 - u ≤ (u 1 - ≤ (& 16 / 2 u + (3^* 23 / 1 u + (3^* - 1 - u ≤ (u 0 ≤ (& 23 / 1 u - (3^* - 16 / 2 u - (3^* + 0 u ≤ (u 1 ≤ (& 16 / 2 u - (3^* 1 u ≤ (u 2 ≤ (& 0 cases :

The Gaussian b-spline

Gaussian basis 'b-spline'

2 - x ≤ x 2 ≤ & x b Undefined if

exact acc 100 x b x 2 - 2 int 1.000000000000000

Explore/resume a NON-linear integration

1. user model parameter n 2.5 :

ε 1012 -^:

t0 1 time :

2. kernel integrand n t f t n^t exp 1 - / :

3. user non-linear path x p 1 0.5 x + / :

4. non-linear integral x Nlin n t f t ε x p int :

Reconciliate the canvas

Fnct x ε 3 π * (ε 0.025 + (range : i 1 x rows range Fnct i el n x i el f eval : for x Fnct augment 31 line :

NONlin x ε 3 π * (ε 0.025 + (range : i 1 x rows range F i el x i el Nlin eval : for x F augment 31 line :

plot, Thiele

data char size clr plot k 1 data rows range r3 k el char : r4 k el size : r5 k el clr : 13 mat for data r3 r4 r5 augment 21 line :

X Y Δ Thiele ε 1018 -^: n X rows : i 1 n range M i 1 el Y i el : for j 2 n range i j n range dM M i j 1 - el M j 1 - j 1 - el - : dM dM abs ε < ε M i j el X i el X j 1 - el - dM / : if : 21 line for for i 1 n range c i 1 el X i el : c i 2 el M i i el : 21 line for Δ 0 ≡ c M if 61 line :

Create a data set of the NON-linear Integration path

create a data set  =>
for Thiele, to resume
the model as an export XY ε ε Nlin 0.25 0.25 Nlin 0.5 0.5 Nlin 0.75 0.75 Nlin 11 Nlin 1.5 1.5 Nlin 22 Nlin 33 Nlin 44 Nlin 66 Nlin 102 mat :

XY 11012 -^* 1.0475 0.25 0.4963 0.5 0.2755 0.75 0.1703 1 0.1136 1.5 0.0589 2 0.035 3 0.0157 4 0.0086 6 0.0035 102 mat

Setup for Thiele suite, export

XY XY :

vx XY 1 col : vy XY 2 col : 12 mat

evaluate/isolate for the 'K' Thiele vx vy 0 Thiele : 11012 -^* 1.0475 0.25 0.4536 - 0.5 1.2885 - 0.75 5.8454 1 0.2533 1.5 316.1118 - 2 0.0134 - 3 2201.4533 4 0.0012 6 74810.8961 - 102 mat

XY XY . 10 black plot :

Thiele triangular table vx vy 1 Thiele

Thiele 'K' for Cf(x) vx vy 0 Thiele

K Thiele 2 col :

vx vy 0 Thiele

x Cf K 1 el x vx 1 el - K 2 el x vx 2 el - K 3 el x vx 3 el - K 4 el x vx 4 el - K 5 el x vx 5 el - K 6 el x vx 6 el - K 7 el x vx 7 el - K 8 el x vx 8 el - K 9 el x vx 9 el - K 10 el / + / + / + / + / + / + / + / + / + :

L H N xd U 0 : dx H L - N / : i 1 N 1 + (range U i el L dx - (dx i * + : for U 41 line :

L ε : H 8 : N 99 : 13 mat U L H N xd : 21 mat

Thiele continued fraction i 1 U rows range vx i el U i el : vy i el U i el Cf eval : 21 line for vx vy augment 31 line :

Fnct Thiele XY 3 1 sys

Integrand
NL-fitted
support pts Fnct Thiele XY 31 sys

1 time t0 - 4 s *

length Thiele j XY rows : 10

Thiele on length of support points 'j' k X x Cfr N j 1 + el k j 1 + el : c j el N j el : a j el k j el : jj j 1 - (0 range c jj el a jj 1 + el : a jj el k jj el x X jj el - c jj el / + : 21 line for a 1 el 51 line :

symbolic expansion of Thiele continued fraction algo style k X x Cfr k 1 el x X 1 el - k 2 el x X 2 el - k 3 el x X 3 el - k 4 el x X 4 el - k 5 el x X 5 el - k 6 el x X 6 el - k 7 el x X 7 el - k 8 el x X 8 el - k 9 el x X 9 el - k 10 el / + / + / + / + / + / + / + / + / +