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

b99edbc6-7fd1-4640-8293-f4b72b6e125a 40 4 5 decimal

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

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 :

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' ... NUMERIC 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 / : 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 :

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

0 f 403

x Frechet f 0 x 10 Int :

t0 1 time :

12 Frechet 29807

x f x c * x Frechet x c * 2 1 sys

1 time t0 - 18 s *

Frechet is essentially an X^Y
Does not assign/plot/evaluate.
NOT scalar wrt X^Y power fnct. t f eval t 0 x int

x Y D x f :

init 0 :

sol init 0 end 120 x Y D dn_GearsBDF :

sol sol rows row 122980412 mat

x f x c * sol 2 1 sys

Advanced Integration . . . multi use

]]>

3 appVersion 0.98.6179 ≡

Example: 1

solving wrt 'y' a k u t f u y a y - (* / y int k t int ≡ y solve maple :

a k u t f u a 1 k t * a * u / exp + / ln * k t * a * + u / exp

solving wrt 't' a k u t f u y a y - (* / y int k t int ≡ t solve maple :

a k u t f u y a y - / ln * k a * /

a k u x f1 u a 1 k x * a * u / exp + / ln * k x * a * + u / exp :

a k u x f2 u x a x - / ln * k a * / :

a 1 : k 1 : u 1 : 13 mat

4 a k u 10 x * f1 * 1 a k u 1 x * f2 * 2 1 sys

4 a k u 0 f1 * a k u 0.5 f2 21 mat 2021 mat

a k u t Clear 1

Example: 2

u 1 u 2^+ / u 1 y int u sin - u 0 t int ≡ y solve maple 21 t cos - (* exp 2 - 21 t cos - (* exp + (* - sqrt 21 t cos - (* exp / 21 t cos - (* exp 2 - 21 t cos - (* exp + (* - sqrt 21 t cos - (* exp / - 21 mat

u 1 u 2^+ / u 1 y int u sin - u 0 t int ≡ t solve maple 1 y 2^+ ln 2 / 2 ln 2 / - 1 + acos

t Φ 21 t cos - (* exp 2 - 21 t cos - (* exp + (* - sqrt 21 t cos - (* exp / ± 1 t 2^+ ln 2 / 2 ln 2 / - 1 + acos 21 sys :

x Φ 1 el x Φ 2 el 2 1 sys

Example: 3

u u 1 u - (* / u y0 y int ρ u t0 t int ≡ y solve maple 1 ρ 1000 t * - 13236380403827 s * + (* 1000 / exp 1 - y0 + (* +

u u a u - (* / u y0 y int ρ u 0 t int ≡ y solve maple a 1 ρ t * - exp - (* ρ t * - exp y0 * +

1 x n^/ x C0 C1 int 1 y 0 t int ≡ C1 solve maple C0 ln n * 1 C0 - t C0 ln n * exp * 1 - n + (* + / - ln + 1 - n + / exp

1 x n^/ x C0 C1 int k - (1 y 0 t int * ≡ k solve maple C0 1 n -^1 C1 1 n -^C0 1 - n +^* - (* t 1 - n + (* / -

x 2^x y (2^+ 4 - (x diff x y x diff solve maple x x y * 1 - (x diff x y x diff solve maple 21 sys x x y / - x y x / - 21 sys

5 x y 2^* x y sin + x 2^- (x diff x y x diff solve maple 7 x y 4^* x 3^x y * + x + 4 - (x diff x y x diff solve maple x 2^x y 2^- 1 - (x diff x y x diff solve maple 31 sys 2 x * 10 x y * x y cos + / 13 x 2^* x y * + 28 x y 3^* x 3^+ / - x x y / 31 sys

x y Sol u 1 - 1 u 2^+ / u 0 y int 1 u + (u 0 x int + 1 u + (u 0 u int - maple :

12 / 1 y 2^+ ln * y atan - x + x 2^2 / + 4 -

x y Sol 12 / 1 y 2^+ ln * y atan - x + x 2^2 / + 4 - :

x 2 x + (* - u 2^+ 2 y atan * + 1 y 2^+ ln - 2 u * + 2 / -

Sol

x y Sol x solve maple 1 - 92 y atan * + 1 y 2^+ ln - sqrt + 192 y atan * + 1 y 2^+ ln - sqrt + (- 21 mat

x plot 1 - 92 x atan * + 1 x 2^+ ln - sqrt + 192 x atan * + 1 x 2^+ ln - sqrt + (- 21 mat :

x plot 1 el x plot 2 el 2 1 sys

Example: 4

v Vg * v 2^2 / - (1 -^v V0 Vt int K 1 t 0 T int * ≡ Vt solve maple 2 V0 * Vg * K T * Vg * exp * V0 1 - K T * Vg * exp + (* 2 Vg * + /

v t diff K v Vg * v 2^2 / - (* ≡

K Vg V0 t q 2 V0 * Vg * K t * Vg * exp * V0 1 - K t * Vg * exp + (* 2 Vg * + / :

x c 10 x ≤ cases :

10 5 1 x q x c * 10 - 5 1 x q x c * 2 1 sys

a 10 : b 4 : c 0.25 : 13 mat

Verhulst a b c x q x c * 0.125 a b c x q x diff (* x c * 21 sys :

X a b c x q x 2 diff 0 ≡ x 01 solve : 0.086

Verhulst X 15 X 5 - 2 2 mat 2 1 sys

Example: 5

DE Verhulst [Logistic]

a k t u y Clear 1

a k t f 1 y a y - (* / y int k t int ≡ y solve maple :

a k t f a k t * a * exp * 1 k t * a * exp + /

a k t f a k t * a * exp * 1 k t * a * exp + / :

a k y0 t f1 1 u a u - (* / u y0 y int k u 0 t int ≡ y solve maple :

a k y0 t f1 y0 a * k t * a * exp * a y0 1 - k t * a * exp + (* + /

a k y0 t f1 y0 a * k t * a * exp * a y0 1 - k t * a * exp + (* + / :

k yo uo t f2 1 u 1 u - (* / u yo y int k u uo t int ≡ y solve maple :

k yo to t f2 yo k t to - (* exp * 1 yo 1 - k t to - (* exp + (* + /

a k yo uo t f3 1 u a u - (* / u yo y int k u uo t int ≡ y solve maple :

a k yo to t f3 yo a * k a * t to - (* exp * a yo 1 - k a * t to - (* exp + (* + /

r to yo t Logistic 111 yo / 1 - (r t to - (* - exp * + / :

yo r t to - (* exp * 1 yo - yo r t to - (* exp * + / 111 yo / 1 - (r t to - (* - exp * + / ≡

NOT from Smath

modelling hysterisis r 1 : yo 0.5 : to 10 : ε 0.75 : 14 mat

a k t f a k t * a * exp * 1 k t * a * exp + / :

a k y0 t f1 y0 a * k t * a * exp * a y0 1 - k t * a * exp + (* + / :

r to yo x Logistic r to ε - yo x Logistic r to ε + yo x Logistic to yo + 15 black 1 5 mat 4 1 sys 1.25 0.5 x f 1.25 0.75 x f 1.25 1 x f 3 1 sys 1.25 0.5 0.025 x f1 1.25 0.75 0.025 x f1 1.25 1 0.025 x f1 3 1 sys

Example: 6

Random example...

u C x Clear 1

SYMBOLIC u 1 + (- 1 u 2^+ / u int 1 x / x int C + ≡ maple 1 u 2^+ ln 2 u atan * + 2 / - C x ln + =

x f 1 x 2^+ ln 2 x atan * + 2 / - : x q 1 x 2^+ ln 2 x atan * + 2 / : 21 sys

X x f x q - 0 ≡ x 4 - 3 - solve : 3.503 -

x f x q X 1 X 1 - 2 2 mat 3 1 sys

Example: 7

Theoretical free fall

... SYMBOLIC ...
parabolic free fall m g * m g * k v * - / v int g 1 t int * - Vo + ≡ v solve maple m g * k Vo - g t * + (* m g * / exp - k /

m g k Vo t freefall m g * k Vo - g t * + (* m g * / exp - k / :

m 1 : g 9.81 : k 0.9071 : Vo 1 : 14 mat

u 0 x ≤ :

X m g k Vo x freefall 0 ≡ x 05 solve : 2.619

m g k Vo x freefall u / 0 g + 40 red 1 5 mat X 0 + 40 red 1 5 mat 3 1 sys

m g k Vo 0 freefall 9.81

m adjust gravity g [m/s²] ≡ k drag coefficient ≡ Vo initial velocity [m/s] ≡ 31 line

b Clear 1

x Γ x Gamma :

α 1.354 :

β 190 :

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

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

cumulative on canvas
depends upon WS 'acc'
& Simpson compatibility x Q u f u 0 x int :

x f x α β PDF :

undetected conflict ??? x cumulative f 0 x 50 Int ≡

t0 1 time :

The incremental cumulative companion L Δ N n Σ X 0 : 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 Int : for X Y augment 31 line :

f 0 x 10 Int x c *

1 time t0 - min 0.3

t0 1 time :

Σincremental 01010010 Σ :

1 time t0 - min 0.1

Project NOTES: 1. Gamma PDF is normalized, i.e: on long range, cumulative is asymptotic to 1. 2. Decimate the quick ODE. 3. n[10] Σincremental ... plenty

Σincremental

x Y D x f :

init 0 :

t0 1 time :

sol_100 init 0100100 x Y D dn_GearsBDF : sol_200 init 0200200 x Y D dn_GearsBDF : sol_2000 init 020002000 x Y D dn_GearsBDF : 31 sys

1 time t0 - 1.2 s *

Smath 6179 Cumulative ODE ... export

sol_200 sol_2000 2 1 sys x f x c *

data n Decimate collect every 'n' pairs range 1 data rows n range : k 1 : i range u k el data i 1 el : v k el data i 2 el : k k 1 + : 31 line for u v augment 31 line 31 line :

Smath 6179 Cumulative ODE ... export

sol_100 1 1 sys

# f Σ ε 1012 -^: U ε 5 0.025 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 x :

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

Σ U 2 - 2 1.95 - 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 :

a t u Clear 1

The widow form ... absent of 'u' in integrand

t0 1 time :

a t γ t a 1 -^t - exp * u 0 t int : x c 10 x ≤ (x 10 ≤ (& cases : 21 sys

plot 2 x γ x c * 1.5 x γ x c * 1 x γ x c * 0.5 x γ x c * 11 e / + 20 black 15 mat 51 sys :

plot

Observe these two models:

Integrator accepts [t, u], as well [x, u]

]]>

1 time t0 - 202 s *

evaluate 1rst @ reduced integral accuracy 50 ... timing is acceptable over tha canvas size.
Te cumulative is about exclusive of running integrator again ... set cumulative via ODE. m n t sol_5 m n C - (24 n * m 2^- (sqrt * u sqrt 0.5 m * t u - (* - exp * 0.5 4 n * m 2^- (sqrt * t u - (* sin * u 0 t int * (* :

t0 1 time :

m 0.5 - : n 1 : 12 mat m n C 14 n * - m 2^+ (/ : x c 10 x ≤ (x 15 ≤ (& cases : x f m n x sol_5 x c * : 41 sys

Scalar wrt 't' explicit form of Dirac integral

x f

1 time t0 - 83 s *

t f m n t sol_5 :

t Y D t f :

init 0 :

sol init 015150 t Y D dn_GearsBDF :

m n C 0.2667 -

solves but slow x f 0 ≡ x 015 solve

roots 0 5.986 8.08 11.98 14.919 51 mat ≡

x f sol 2 1 sys

A difficult/uncertain case

t y f 3250341 y * - 9 t * + (* + 3250341 y * 9 t * + (* - / :

init 0 :

y0 x0 x1 n f rkf4 X 1 el x0 : Y 1 el y0 : Δ x1 x0 - n / : 13 mat i 2 n range k1 X i 1 - el Y i 1 - el f eval : k2 X i 1 - el Δ 2 / + Y i 1 - el k1 2 / Δ * + f eval : k3 X i 1 - el Δ 2 / + Y i 1 - el k2 2 / Δ * + f eval : k4 X i 1 - el Δ + Y i 1 - el k3 Δ * + f eval : Y i el Y i 1 - el k1 2 k2 * + 2 k3 * + k4 + 6 / Δ * + : X i el X i 1 - el Δ + : 61 line for X Y augment 31 line :

rkf4 0030500 t y f rkf4 :

Unknown reaction: Yuk [Ru.] 2004

rkf4 15.57 100 15.57 25 - 2 2 mat 12 22.54 20 22.54 2 2 mat 3 1 sys

Whether the integration is right/wrong ? Good question, apparently YES from lab test. This example from Yuk [Ru.] Mathsoft Collaboratoty, Mathcad 11 [2004].

ODE Libary

Matlab solvers

Intel solvers

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

ONE  possible limited solver. # # # # # mwode113 # # # # # mwode15s # # # # # mwode23 # # # # # mwode23s # # # # # mwode45 51 line

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

# # # # # mk52lfa # # # # # mk52lfn 21 line

Under "while", the module outputs the last RK4 value. 
Computes on-line integral value(s) @ user 'u' ... computes for a given vector. init t1 t2 n f RK4 vx init : t t1 : h t2 t1 - n / : 13 mat t t2 < k1 t vx f h * eval : k2 t h 2 / + vx k1 2 / + f h * eval : 12 mat k3 t h 2 / + vx k2 2 / + f h * eval : k4 t h + vx k3 + f h * eval : 12 mat vx vx 16 / k1 2 k2 * + 2 k3 * + k4 + (* + eval : t t h + : 12 mat 31 line while vx 31 line :

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

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' ... NUMERIC  Respect some rules for ± ∞ [a, b] 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 / : 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 :

x f 1 x 1 - (2^- sqrt :

all dn_ODE solvers fail x Y D x f :

sol 0021000 x Y D rkm9st :

π 2 / 1.570796326794900

sol sol rows row 2 1.570793 12 mat

good for plot ...otherwise 5 D

x f sol 2 1 sys

t0 1 time :

f 022000 Int 1.57079632679491

1 time t0 - min 0.1