Σalgo(lo,mesh,hi,∫)

ac43f866-f6c6-4bbc-b328-bfbb585da7c7 29 4 5 decimal

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

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 :

Gaussian basis 'b-spline'

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

0 b 0.6667

x f x b :

Algo style Memo formula of incremental analysis of this demo example h 0.0001 : (increment ≡ x h f deriv 0.5 x h - f x f - x h - (x - / x f x h - f - x x h - (- / + (* : x h f D1 x 2 h * - f x 2 h * + f - 4 h * / - : (1rst Derivative ≡ x h f D2 x 2 h * + f - 16 x h + f * + 30 x f * - 16 x h - f * + x 2 h * - f - 12 h 2^* / : (2nd Derivative ≡ x h f h_int h x h + f x f - (* 2 x f * + (x * : (Integral Curve [antiderivative] ≡ 51 line

constrain the canvas x c 12 - x ≤ (x 2 ≤ (& cases :

x cumul u b u 2 - x int x c * :

x Maths x f x c * x h f D1 x c * x h f D2 x c * x h f h_int x c * x cumul x c * 61 sys :

Σ U 2 - 2 1.99 - range : i 1 U rows range Cumul i el x b x U 1 el U i el int eval : for U Cumul augment 31 line :

x Maths

isolate/accelerate Σ Σ eval :

U 2 - 2 1.99 - range : data U U f vectorize augment : 21 line

NEXT FD [Finite Differences]is ... 1. very fast and accurate. 2. almost immute to noisy data. 3. offered algo style

data Σ FD Cumulative Integration vx data 1 col : vy data 2 col : 12 mat Y 0 : Y 1 el 0 : 12 mat 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 51 line :

data Σ Σalgo(lo,mesh,hi,∫)

plug integral curve

x Clear 1

Cumulative integration algo style ... lo mesh hi ∫ Σalgo cumullative Integrat'n algo style Y 0 : data 0 : 12 mat x lo hi mesh range : data x ∫ vectorize augment : ux data 1 col : uy data 2 col : Y 1 el 0 : i 2 ux rows range Δx i el ux i el ux i 1 - el - : Y i el Y i 1 - el 12 / uy i 1 - el uy i el + (Δx i el * (* + : 21 line for ux Y augment 91 line :

μ 0 : σ 2 : 12 mat ∫ x μ σ pdf.Normal eval : 21 sys

CDF 10 - 9.9 - 10 ∫ Σalgo :

∫ CDF 2 1 sys

Multiple ordinary Simpson ... unpredictable accuracy

Depends upon option Simpson 'acc' # f Σ ε 1012 -^: U ε 5 0.05 range : i 1 U rows range Cumul i el x f x U 1 el U i el int eval : for U Cumul augment 41 line :

Assign user as function

x q 1 x 2^1 + / :

x f 0.25 x * :

x f Σ x ln Σ x q Σ 3 1 sys

will integrate scalar splines fnct(s)... ainterp(X Y x) cinterp(X Y x)

Σ U 2 - 2 1.9 - range : i 1 U rows range Cumul i el x erf vectorize x U 1 el U i el int eval : for U Cumul augment 31 line :

Σ 2 1 sys

∫ x ln eval : lo 1012 -^: mesh 0.01 : hi 5 : 13 mat Σalgo_1 lo mesh hi ∫ Σalgo : 31 sys

∫ 0.25 x * eval : lo 0 : mesh 0.01 : hi 5 : 13 mat Σalgo_2 lo mesh hi ∫ Σalgo : 31 sys

∫ 1 x 2^1 + / eval : lo 0 : mesh 0.01 : hi 5 : 13 mat Σalgo_3 lo mesh hi ∫ Σalgo : 31 sys

collect Σalgo_1 Σalgo_2 Σalgo_3 31 sys :

collect

Numerical pitfall ... Simpson rejects X^Y in f(x)

α 1523043253 - 27354.7 - 43256.75 25388.3 - 40061 mat :

Fréchet: extreme value-II x f α 1 el α 2 el x α 4 el - (/ (α 5 el^* α 2 el x α 4 el - (/ (α 3 el^- exp * α 6 el + :

end 12 :

x c 10 x ≤ (x end ≤ (& cases :

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 :

data U 0 end 0.1 range : i 1 U rows range vy i el U i el f eval : for U vy augment 31 line :

Finite differences Cumulative Integral. @ this point cinterp for in-situ eval.

data data Σ 2 1 sys x f x c * x f 2 1 sys

x Y D x f : sol 005100 x Y D dn_GearsBDF : 21 sys

u sol 101 row 1 el : v sol 101 row 2 el : data Σ 51 row 31 mat 51710951710912 mat 31 mat

sol data Σ u v + 20 black 1 5 mat 3 1 sys

Integrator refuses to calculate  ... "x_ not defined"... 
You can rescue integration from finite differences ... or ʃkernel(f,a,b,n) x f 1523043253 - x 43256.75 - (/ (25388.3 -^* 43253 - x 43256.75 - (/ (27354.7 -^- exp * :

ʃkernel(f,a,b,n) ... Simpson/Romberg 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 
Integration accuracy from user command line 'n' 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 / : 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 * + : 32 / b a - (* Simpson * 41 line :

x f

X^Y  not scalar wrt 'x'
bug ►  Smath Simpson.
rescue ... ʃkernel(f,a,b,n) x f x 010 int #

1/1 Mathcad 11 0 f 3.0597

sanity Mathcad  11 adaptive f 010100 ʃkernel 24556

Finite differences will accumulate as well

plot scale for convenience scale 0.25 :

Σ U 016 0.25 range : j 1 U rows range Cumul j el f 0 U j el 10 ʃkernel eval : for U scale Cumul * augment 31 line :

x f Σ 2 1 sys

vx Σ 1 col : vy Σ 2 col scale 1 -^* : 12 mat x Σ vx vy x ainterp : 21 mat

true Σ @ user sanity above ▲ 10 Σ 24556

The converter offers 3 bonuses ... 1. The kernel reduces long ranges of the variate to [0..1] Thus increasing accuracy vs compatible ordinary Simpson. 2. It enables ordinary Simpson to accept X^Y otherwise rejected. 3. Adds variable 'n' meshing local in-situ use, independent of each application in same document. This Romberg/Simpson companion has been tested exhaustively. Found nearly equivalent to Carlos Adaptive and Mathcad Adaptive. It supports singular end point and infinite limit. As demonstrated, it simply recasts the native Smath simpson.

ODE integrators pitfall ...

4 appVersion 0.98.6179.21440 ≡

Test function x f x sqrt exp :

x Y D x f :

a 0 : b 2 : 12 mat

sol 0 a b 10000 x Y D al_rkckadapt :

exact 5.40752818465627

very good ! sol sol rows row 2 5.40752818457808 12 mat

sol 0 a b 10000 x Y D rkfixed :

3 ½ D sol sol rows row 2 5.40722 12 mat

x f x 0.5 1 -^x - exp * eval : a 1015 -^: b e : 12 mat 21 sys

x Y D x f :

close sol 0 a b 10000 x Y D dn_GearsBDF :

π sqrt b sqrt erf * 1.73750207472565

sol sol rows row 2.7183 1.7382 12 mat

sol 0 a b 10000 x Y D rkfixed :

f....n freak ! sol sol rows row 2.7183 8597.6753 12 mat

sol 0 a b 10000 x Y D al_rkckadapt :

4 ½  D sol sol rows row 2.71828182845905 1.7375246742603 12 mat

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 :

known exact
acc 100 u b u 2 - 2 int 1.000000000000000

x Y D x b :

4 ½ D sol 02 - 210000 x Y D dn_GearsBDF :

sol sol rows row 1.99999999999969 0.999920601698221 12 mat

rkfixed  exact! sol 02 - 2100 x Y D rkfixed :

sol sol rows row 2.000000000000000 1.000000000000000 12 mat

sol 02 - 210000 x Y D al_rkckadapt :

slow but converges exact sol sol rows row 2 0.999999999999998 12 mat

stiff ... NON stiff soup ! # # # # # mwode113 # # # # # mwode15s # # # # # mwode23 # # # # # mwode23s # # # # # mwode45 51 line

Useless solvers. # # # # # al_rkckadapt # # # # # dn_AdamsMoulton # # # # # dn_ExplicitRK45 # # # # # dn_GearsBDF # # # # # dn_ImplicitRK5 51 line

# # # # # Rkadapt # # # # # rkfixed # # # # # rkm9mka # # # # # rkm9mkn # # # # # rkm9st 51 line

1.73750207472565

stiff Jacobian # # # # # mk52lfa # # # # # mk52lfn 21 line

x f x 0.5 1 -^x - exp * eval : a 1018 -^: b e : 12 mat 21 sys

N 1000 :

Fail to ▲ solve

x Y D x f :

sol 0 a b N x Y D al_rkckadapt :

sol sol rows row 2.718282 1.737526 12 mat

sol 0 a b N x Y D dn_AdamsMoulton :

sol sol rows row 2.718282 1.737568 12 mat

sol 0 a b N x Y D dn_ExplicitRK45 :

sol sol rows row 2.718282 1.737519 12 mat

sol 0 a b N x Y D dn_GearsBDF :

sol sol rows row 2.718282 1.738344 12 mat

sol 0 a b N x Y D dn_ImplicitRK5 :

sol sol rows row 2.718282 1.737499 12 mat

x f x 0.5 1 -^x - exp * eval : a 1018 -^: b e : 12 mat 21 sys

# # # # # Rkadapt # # # # # rkfixed # # # # # rkm9mka # # # # # rkm9mkn # # # # # rkm9st 51 line

x Y D x f :

sol 0 a b N x Y D Rkadapt :

freak ! sol sol rows row 2.718282 469.379272 12 mat

sol 0 a b N x Y D rkfixed :

freak ! sol sol rows row 2.718282 2.718283 106^* 12 mat

sol 0 a b N x Y D rkm9mka :

close ! sol sol rows row 2.718282 1.737671 12 mat

sol 0 a b N x Y D rkm9mkn :

close ! sol sol rows row 2.718282 1.737671 12 mat

sol 0 a b N x Y D rkm9st :

close ! sol sol rows row 2.718282 1.737671 12 mat

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

f a b ε ∫Adaptive N 100 : 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 21 line eval :

x f x 0.5 1 -^x - exp * eval : a 1019 -^: b e : 12 mat 21 sys

Adaptive Mathcad 1.73750207398205 ≡

Exact 1.73750207472565

t0 1 time :

Damned good !!! f a e 0.0000000001 ∫Adaptive 1.737502074849300

1 time t0 - 8.2 s *

1.73750207472565 1.7375020748493 - 1.24 1010 -^* -

Like most other solvers, ∫Adaptive(f,a,b,ε) converges first ... then diverges wrt ε. 'a' [zero] has some influence. Mathcad/Mathsoft zero was 'guessed' 10^{-18}

EMERITUS Romberg-Simpson extra long range ...

ʃkernel(f,a,b,n) ... Simpson/Romberg 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 
Integration accuracy from user command line 'n' 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 / : 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 * + : 32 / b a - (* Simpson * 41 line :

Gamma(x) Lanzcos Davide

z Lanczos Davde Carpi 20200514 ε 11015 -^* : p 57.156235665862923517 59.597960355475491248 - 14.136097974741747174 0.49191381609762019978 - .33994649984811888699 104 -^* .46523628927048575665 104 -^* .98374475304879564677 - (104 -^* .15808870322491248884 103 -^* .21026444172410488319 - (103 -^* .21743961811521264320 103 -^* .16431810653676389022 - (103 -^* .84418223983852743293 104 -^* .26190838401581408670 - (104 -^* .36899182659531622704 105 -^* 141 mat : g 4.7421875 : z Re 0.5 < y π π z * sin 1 z - Lanczos * / : z z 1 - : x 0.99999999999999709182 : n 1 p length range x x p n el z n + / + : for t z g + 0.5 + : y 2 π * sqrt t z 0.5 +^* t - exp * x * : 51 line if y Im abs ε ≤ y Re y if 61 line :

g=4.7421875, n=15

x Lanczos

2.5 Lanczos 1.32934038817912

2.5 Gamma 1.32934227535731

x Γ x Lanczos :

x α β PDF 1 α Γ β α^* / x α 1 - (^* x β / - exp * :

x λ 10 x ≤ (x 1000 ≤ (& cases :

1. Gamma PDF is normalized, i.e: on long range, cumulative is asymptotic to 1. 2. [n] wrt range of integration. 3. Simpson NOT sanity low range ?

α 1.35 :

β 190 :

x α β PDF x λ *

x f x α β PDF :

t0 1 time :

f 0 x 100 ʃkernel x λ * 250 x α β PDF x λ * * 2 1 sys

1 time t0 - min 3.6

The incremental cumulative companion L Δ N n Σ X 0 : Y 0 : 12 mat i 1 N range X i 1 + el L Δ i * + : X 1 el L : 21 line for X 31 line i 1 X rows range Y i el f L X i el n ʃkernel : for X Y augment 31 line :

Σincremental 01010010 Σ :

Σincremental

x Y D x α β PDF :

init 0 :

t0 1 time :

sol_300 init 0300300 x Y D dn_GearsBDF : sol_2000 init 020002000 x Y D dn_GearsBDF : 21 sys

1 time t0 - 0.7 s *

Smath 6179 Cumulative ODE ... export

sol_300 sol_2000 2 1 sys

min(,) undocumented x f 2110 x - (* 260 x - 2 / 5128 x - (* 4 / 31 mat min :

x f x bound 412004164 - 22 mat 802008064 - 22 mat 41 sys

x bound 1080 x - (* 7 / :

x f x bound 41 200 41 64 - 2 2 mat 80 200 80 64 - 2 2 mat 4 1 sys

accepted
acc 50 x f x bound - (x 16.923 - 185 int 6911.5

x f x bound - 0 ≡ x 32 - 0 solve x f x bound - 0 ≡ x 160200 solve 21 mat 16.923 - 18521 mat

L H N Φ i 1 N range ux i el L i N / H L - (* + : uy i el ux i el f ux i el bound - (eval : 21 line for ux uy augment 21 line :

data 16.923 - 185400 Φ :

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 :

data Σ

n data Σ rows :

data Σ n row 185 6911.4 12 mat

Adaptive on canvas ... user point value

∫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 :

x 0.5 1 - ^ x - exp * π sqrt x sqrt erf * 2 1 sys

eval(,) gains 1 min
.....
iter 100 plenty
'ε' canvas accuracy
.....
'a' infinite limit
'b' user end canvas plot x f x 0.5 1 -^x - exp * eval : iter 100 : ε 0.01 : a 1015 -^: b e : 51 sys

x Canvas f a x iter ε ∫Adaptive :

x c 1 a x ≤ (x e ≤ (& cases :

x λ 1 a x ≤ (x 5 ≤ (& cases :

x f π sqrt x sqrt erf vectorize * - 0 ≡ x 0.25 0.5 solve 0.384

t0 1 time :

Adapt Discrete cumulative on that range ux a 5 0.1 range : j 1 ux rows range uy j el f a ux j el iter 0.0001 ∫Adaptive eval : for ux uy augment 41 line :

isolate < timing Adapt Adapt :

vx Adapt 1 col : vy Adapt 2 col : 12 mat x spline vx vy x cinterp : 21 sys

1 time t0 - min 0.8

REMEMBER ... The Rieman rectangles FD ranges from very accurate to exact. It is almost insensible to noisy data sets. ODE could be the first companion to experiment. Watch ± ∞ limit/singular generally a good candidate, very fast. Sanity check ODE vs Adaptive.

Scalar Canvas(x)*c(x) plots on each canvas pixel: higher timing

x f x Canvas x c * x spline x λ * e e spline + 40 black 1 5 mat 4 1 sys

1 time t0 - min 3.2

ε 1015 -^:

x Y D x f :

init 0 :

t0 1 time :

sol init ε b 200 x Y D dn_GearsBDF :

1 time t0 - 0.5 s *

Smath 6179 Cumulative ODE ... export

sol

n sol rows :

sol n row 2.7183 1.7382 12 mat

b spline 1.7376

Rombrg kernel test function ...

f a b Int x f x a b int 32 / b a - (* w G w 01 int * ≡ memo ≡ 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 * 41 line :

Test function x f x sqrt exp :

t0 1 time :

acc 10000 7 true digits ... 
insufficient integration subspaces.
acc 100 3 decimals x f x 02 int 5.40744702348009

1 time t0 - 0.1 s *

acc 100 8 decimals f 02 Int 5.40752818690417

x f

5.40752818465638

sanity check the quadruplet Mathcad 11 Romberg 5.40752818465632 ≡ Technical accuracy 5.40752818465627 ≡ Smath Romberg 5.40752818465638 ≡ Smath Integrator 5.4075281 ≡ 41 line

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 / : 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 * + : 32 / b a - (* Simpson * 41 line :

exact  5.40752818465627 f 022000 ʃkernel 5.40752818465638

x f x 0.5 1 - ^ x - exp * : ε 10 15 - ^ : x f x ε e int 143268 = π sqrt e sqrt erf * 1.73750207472565 = SimpsonGloriousFreak 5 1 line

only ... ∫Adaptive(f,a,x ...)

x f