Int(f,a,b,n) ... Simpson/Romberg Kernel

44946269-ddd5-4483-ae33-949f943354ad 42 4 5 decimal

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

The first tool integrators

4 appVersion 0.98.6179.21440 ≡

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

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

The 3 first tool Integrators

1. Int(f,a,b,n) . . large applications...

2. ∫Adaptive(f,a,b,ε) most difficult

3. FD noisy data sets.

4. ODE solvers are ? ? ? candidates. .

]]>

Test function x f x sqrt exp :

Smath/Maple agree with Mathematica x f x int maple 2 x sqrt exp * 1 - x sqrt + (*

x m 2 x sqrt exp * 1 - x sqrt + (* :

a 0 :

b 2 :

sanity check in  absolute accuracy b m a m abs + 5.40752818465627

Test acc 10000 x f x a b int 5.40752810347196 ≡

x f a b 2000 Int 5.40752818465635

The crucifying region ... peanut !

x f 0 1 o 50 red 1 5 mat 2 1 sys

x Y D x f :

sol 0 a b 200 x Y D dn_GearsBDF :

Could rescue not having better sol sol rows row 2 5.4067 12 mat

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

eval(,) gains 1 min
.........
'a' infinite limit
'b' end canvas plot x f x 0.5 1 -^x - exp * eval : a 1015 -^: b e : 12 mat 21 sys

recast ∫integrator x Canvas f a x 0.01 ∫Adaptive :

t0 1 time :

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

Canvas(x) is scalar wrt pixels runs 126 *∫Adaptive(f,a,x,ε) for canvas plot ... ε 0.01

x f x Canvas x c * 2 1 sys

1 time t0 - min 2.0

t0 1 time :

f a b 100000 Int 1.737493

∫

1 time t0 - min 9.8

acc 10000
 ... freak ! x f x a b int 1434 ≡

productive/slow f a b n Int

t0 1 time :

π sqrt b sqrt erf * 1.73750207472565

for accurate ... ∫ set ε ... 0.000001 f a e 0.00000001 ∫Adaptive 1.73750204

Adaptive Mathcad 1.73750207398205 ≡

1 time t0 - 4.1 s *

x sqrt exp

f a b Int ordinary Romberg depends upon Simpson 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 * 51 line :

Int(f,a,b,n) near exact this demo accuracy Int(f,a,b,n)converges, not limited by plugin integrator ... Int(f,a,b) [10000]

Test function x f x sqrt exp :

t0 1 time :

x f 2 5 2 1 - 2 2 mat 2 1 sys

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

f 02 Int 5.40752818690417

1 time t0 - 0.7 s *

x Φ x sqrt :

acc 10000 back to
diificult region x Φ x 0100 int 666.666637963573 ≡

Declared exact Φ 010010000 Int 666.666666666664 ≡

Conclusion . . .

1. Simpson plugin ranges from exact to pure crap.

2. Sanity check your integrand with these tools.

3. exp(sqrt(x)) looks naive, but difficult.

4. x^{0.5-1}*exp(-x) is one of the most difficult encountered.

5. Both integrators Int(f,a,b,n) ... ∫Adaptive(f,a,x,ε) agree ~

6. Int(f,a,b,n) is a desirable plugin, ∫Adaptive(f,a,x,ε) as well.

]]>

Two exact Simpson integrations ...

x λs 12 - x ≤ (x 2 ≤ (& cases :

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'

x b x λs * 2 - 0 + 40 red 1 5 mat 2 0 + 40 red 1 5 mat 3 1 sys

0 b 0.6667

Exact by design, acc 100
Simpson:  fair all Gaussians u b u 2 - 2 int 1.000000000000000 ≡

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 :

definite integrator algo style # ∫

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 x ln b 0.75 b 0.25 - 2 2 mat 6 1 sys

transit algo wrt scalar 'x'

Smath code ..... z ∫ 1 u / u 1 x int : x Σ z ∫ : 21 line

2 Σ 0.5 Σ 2 ln 0.5 ln 41 mat 0.693147180579478 0.693147180579478 - 0.693147180559945 0.693147180559945 - 41 mat

an exact Simpson
acc  ... 10000 x f x 13 int 3 ln 21 mat 1.09861228866811 1.09861228866811 21 mat ≡

data Σ Finite diff's 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 :

The FD [Finite Differences] is the sole tool to integrate data set. This version is the Riemann rectangle, works well. It is almost immute to noise, exclusive of above examples.

data [XY]

XY pairs raw data data 0 0.1565 0.0982 0.0846 0.1963 0.1657 0.2945 0.5194 0.3927 0.4005 0.4909 0.3374 0.589 0.7377 0.6872 0.8511 0.7854 0.7138 0.8836 0.9814 0.9817 0.6979 1.0799 0.8564 1.1781 0.9703 1.2763 0.964 1.3744 0.8899 1.4726 0.7495 1.5708 1.2045 1.669 1.1278 1.7671 0.9497 1.8653 1.1751 1.9635 1.0413 2.0617 0.7379 2.1598 1.0015 2.258 0.9242 2.3562 0.8504 2.4544 0.6311 2.5525 0.7082 2.6507 0.6188 2.7489 0.263 2.8471 0.3845 2.9452 0.1969 3.0434 0.1199 3.1416 0.2277 3.2398 0.0956 3.3379 0.1321 - 3.4361 0.0562 - 3.5343 0.2352 - 3.6325 0.6592 - 3.7306 0.3314 - 3.8288 0.6555 - 3.927 0.6554 - 4.0252 0.9939 - 4.1233 0.773 - 4.2215 0.7014 - 4.3197 0.6809 - 4.4179 0.8392 - 4.516 0.8484 - 4.6142 1.1527 - 4.7124 1.2204 - 4.8106 1.0773 - 4.9087 1.2067 - 5.0069 0.981 - 5.1051 0.7848 - 5.2033 0.7194 - 5.3014 1.0022 - 5.3996 0.7225 - 5.4978 0.6406 - 5.596 0.565 - 5.6941 0.7197 - 5.7923 0.7004 - 5.8905 0.4193 - 5.9887 0.5055 - 6.0868 0.1086 - 6.185 0.0645 6.2832 0.0435 - 6.3814 0.1926 6.4795 0.1918 6.5777 0.3931 6.6759 0.5767 6.7741 0.3964 6.8722 0.4125 6.9704 0.7767 7.0686 0.8656 7.1668 1.013 7.2649 0.7479 7.3631 0.9702 7.4613 0.9714 7.5595 1.1964 7.6576 1.1391 7.7558 0.7707 7.854 1.1144 7.9522 1.0553 8.0503 1.0986 8.1485 1.1693 8.2467 0.994 8.3449 0.6578 8.443 0.9643 8.5412 0.5882 8.6394 0.8214 8.7376 0.3969 8.8357 0.4731 8.9339 0.2995 9.0321 0.609 9.1303 0.4934 9.2284 0.0394 - 9.3266 0.1149 9.4248 0.0574 - 9.523 0.0479 9.6211 0.2211 - 9.7193 0.3919 - 9.8175 0.5117 - 9.9157 0.6919 - 10.0138 0.6571 - 10.112 0.5298 - 10.2102 0.8961 - 10.3084 0.8912 - 10.4065 0.764 - 10.5047 0.7329 - 10.6029 0.9278 - 10.701 0.966 - 10.7992 1.2277 - 10.8974 0.9181 - 10.9956 0.8389 - 11.0937 1.1849 - 11.1919 1.0741 - 11.2901 0.8036 - 11.3883 1.0069 - 11.4864 0.7579 - 11.5846 0.9891 - 11.6828 0.8423 - 11.781 0.4954 - 11.8791 0.7029 - 11.9773 0.317 - 12.0755 0.2899 - 12.1737 0.6252 - 12.2718 0.5115 - 12.37 0.0291 - 12.4682 0.1073 1282 mat :

data rows 128

Σ data Σ :

data Σ 2 1 sys

Σ Σ rows row 12.4682 0.2811 12 mat

Σ 2 col Mean 1.2

Σ 00 0.0982 0.0118 0.1963 0.0241 0.2945 0.0578 0.3927 0.1029 0.4909 0.1392 0.589 0.1919 0.6872 0.2699 0.7854 0.3467 0.8836 0.43 0.9817 0.5123 1.0799 0.5887 1.1781 0.6783 1.2763 0.7733 142 mat

Exercise the ODE solvers

Revisit ODE solvers as integrators

Test function x f x sqrt exp :

x Y D x f :

a 0 : b 2 : 12 mat

sol 0 a b 1000 x Y D dn_GearsBDF :

Could rescue not having better sol sol rows row 2 5.4065 12 mat

5.40752818465627

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

x Y D x f :

Not much accurate ! sol 0 a b 1000 x Y D dn_GearsBDF :

sol sol rows row 2.7183 1.7383 12 mat

π sqrt b sqrt erf * 1.7375

rkfixed fails miserably ! sol 0 a b 400 x Y D rkfixed :

sol sol rows row 2.7183 214900.6646 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 :

x Y D x b :

rkfixed fails miserably ! sol 02 - 21000 x Y D dn_GearsBDF :

sol sol rows row 2 0.999469673868317 12 mat

Romberg kernel . . .

Export/enable this section ... set acc 100 . . . 10000

]]>

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

a 0 : b 1 : 12 mat

w s b a - (34 / 2 w * 1 - (* 14 / 2 w * 1 - (3^* - 12 / + (* a + :

x f x sqrt :

w G w s f 12 w * 1 - (2^- (* :

x s x c * x G x c * x f x c * 3 1 sys

u f u a b int 0.66663796357003

32 / b a - (* w G w 01 int * 0.666666667142827

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

w s b a - (34 / 2 w * 1 - (* 14 / 2 w * 1 - (3^* - 12 / + (* a + :

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

w G w s f 12 w * 1 - (2^- (* :

x s x λ * x G x λ * x f a x ≤ ( x b ≤ ( & * 3 1 sys

u f u a b int 1.4 105^*

32 / b a - (* w G w 01 int * 1.727983

π sqrt b sqrt erf * 1.737502

acc 10000
 ... freak ! x f x a b int 1434.3831 ≡

Integrator refuses to calculate  ... "x_ not defined"... 
You can rescue integration from finite differences ... or Romberg. x f 1523043253 - x 43256.75 - (/ (25388.3 -^* 43253 - x 43256.75 - (/ (27354.7 -^- exp * :

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

x f

X^Y  not scalar wrt 'x'
bug ►  Smath Simpson. x f x 010 int #

1/1 Mathcad 11 0 f 3.0597

sanity Mathcad  11 adaptive f 010100 Int 24556

Finite differences will accumulate as well

t0 1 time :

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 Int 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

x f x a b int f a b 10 Int :

x f x 05 int 15110

1 time t0 - 22.729 s *

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 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.