Vrkfixed(y0,t0,t2,N,f)

cdfce7c9-cfb1-4ce5-9ab4-601325aac511 52 4 5 decimal

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

y0 t0 t2 N f Vrkfixed t 1 el t0 : Y 1 el y0 : Δ t2 t0 - N / : 13 mat i 2 N range k1 t i 1 - el Y i 1 - el f Δ * eval break try : k2 t i 1 - el Δ 2 / + Y i 1 - el k1 2 / + f Δ * eval break try : k3 t i 1 - el Δ 2 / + Y i 1 - el k2 2 / + f Δ * eval break try : k4 t i 1 - el Δ + Y i 1 - el k3 + f Δ * eval break try : Y i el Y i 1 - el 16 / k1 2 k2 * + 2 k3 * + k4 + (* + : t i el t i 1 - el Δ + : 61 line for t Y augment 31 line :

Was designed by Valery Ochov [Mathsoft Collaboratory] All Mathcad ODE solvers failed.

ODE parametres β 3.5 12 0.25 41 mat :

Smath 1/1 Valery MCD ... ~ mwode113

t Y f t Y 1 el (β 1 el^^t Y 2 el (β 2 el^^+ t Y 1 el (β 3 el^^t Y 2 el (β 4 el^^+ 21 mat :

t0 0 : t2 2 : N 1000 : 31 line

y0 0.5 121 mat :

[t2, N} govern Δ
reflexive wrt plots Δ t2 t0 - N / : 0.002

rkfΔ y0 t0 t2 N f Vrkfixed :

rkfΔ 0 0.5 121 mat 0.002 0.501 121 mat 0.004 0.502 1.001 21 mat 0.006 0.503 1.001 21 mat 0.008 0.505 1.002 21 mat 0.01 0.506 1.003 21 mat 0.012 0.507 1.003 21 mat 0.014 0.509 1.004 21 mat 0.016 0.51 1.005 21 mat 0.018 0.511 1.005 21 mat 0.02 0.513 1.006 21 mat 0.022 0.514 1.007 21 mat 122 mat

data UnestRow v data 1 el : j 2 data cols range v v data j el stack : for v 31 line :

sol rkfΔ 1 col rkfΔ 2 col transpose UnestRow augment :

XY sol sol rows row (: 1.1800 5.2013 2.5372 13 mat

jnct XY 11 el XY 13 el 12 mat :

jnct jnct eval : 1.18 2.5372 12 mat

plot sol 1 col sol 2 col augment sol 1 col sol 3 col augment 21 sys :

sol 0 0.5 1 0.002 0.5009 1.0003 0.004 0.5021 1.0008 0.006 0.5034 1.0013 0.008 0.5046 1.0019 0.01 0.506 1.0025 0.012 0.5073 1.0032 0.014 0.5087 1.0038 0.016 0.51 1.0045 0.018 0.5114 1.0053 103 mat

plot jnct 2 1 sys

N 100 :

ODE113 y0 t0 t2 N t Y f mwode113 :

ODE23 y0 t0 t2 N t Y f mwode23 :

ODE23s y0 t0 t2 N t Y f mwode23s :

ODE45 y0 t0 t2 N t Y f mwode45 :

}NO complain but fail

]]>

ODDE15s y0 t0 t2 N t Y f mwode15s :

These 3 have survived ODE113 ODE113 rows row ODE23 ODE23 rows row ODE23s ODE23s rows row 31 sys 1.16 2.1537 2.4644 13 mat 1.16 2.1524 2.4643 13 mat 1.16 2.1489 2.4638 13 mat 31 sys

They all fail ▼ but not all complain ...

NOT defined ODE15s ODE15s rows row #

't' NOT defined ODE45 ODE45 rows row #

Rkadapt y0 t0 t2 N t Y f al_rkckadapt :

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

Ad y0 t0 t2 N t Y f dn_AdamsMoulton :

Expl y0 t0 t2 N t Y f dn_ExplicitRK45 :

Gears y0 t0 t2 N t Y f dn_GearsBDF :

Impl y0 t0 t2 N t Y f dn_ImplicitRK5 :

No complain but fails Unk y0 t0 t2 N t Y f mk52lfa :

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

Zombie never stops Unk y0 t0 t2 N t Y f mk52lfn :

Intel y0 t0 t2 N t Y f Rkadapt :

Intel y0 t0 t2 N t Y f rkfixed :

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

Intel y0 t0 t2 N t Y f rkm9mka :

}NO complain but fail

]]>

Intel y0 t0 t2 N t Y f rkm9mkn :

Intel y0 t0 t2 N t Y f rkm9st :

Unk Unk rows row #

Intel Intel rows row #

ODE integrators pitfall . . .

In this adventure, two exact results are given.

]]>

4 appVersion 0.98.6179.21440 ≡

Test function x f x sqrt exp :

exact from Maple 5.40752818465627

x Y D x f :

a 0 : b 2 : 12 mat

x f

sol 0 a b 1000 x Y D al_rkckadapt :

very good ! sol sol rows row 2 5.40752818218196 12 mat

x f x a b int 5.4075281

sol 0 a b 1000 x Y D rkfixed :

2 ½ D sol sol rows row 2 5.4044 12 mat

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

π sqrt b sqrt erf * 1.73750207472565

x f

x Y D x f :

2½ D sol 0 a b 10000 x Y D dn_GearsBDF :

x f x a b int 1434

sol sol rows row 2.7183 1.7382 12 mat

sol 0 a b 1000 x Y D rkfixed :

.... Oh ! sol sol rows row 2.7183 85961.2803 12 mat

sol 0 a b 1000 x Y D al_rkckadapt :

4 ½  D sol sol rows row 2.71828182845905 1.73752522069827 12 mat

sol 0 a b 10000 x Y D mwode45 :

mwode23, mwode45 sol sol rows row 2.718282 1.737518 12 mat

1.73750207472565

Matlab 0 a b 1000 x Y D mwode23 :

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

Matlab Matlab rows row 2.71828182845905 1.73757852024273 12 mat

Matlab 0 a b 1000 x Y D mwode45 :

Matlab Matlab rows row 2.71828182845905 1.73751781761842 12 mat

unk 0 a b 1000 x Y D mk521fn :

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

unk unk rows row #

Intel 0 a b 1000 x Y D rkm9mka :

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

Intel Intel rows row 2.71828182845902 1.73750985745675 12 mat

Intel 0 a b 1000 x Y D rkm9mkn :

Intel Intel rows row 2.71828182845902 1.73750985745675 12 mat

1.73750207472565

Intel 0 a b 1000 x Y D rkm9st :

Intel Intel rows row 2.718282 1.73751 12 mat

Some are valid ~ 4 D

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

Dot 0 a b 1000 x Y D dn_ImplicitRK5 :

Dot Dot rows row 2.718282 1.737499 12 mat

Revisit Vrkfixed(y0,t0,t2,N,f)

Vrkfixed(y0,t0,t2,N,f)

y0 t0 t2 N f Vrkfixed t 1 el t0 : Y 1 el y0 : Δ t2 t0 - N / : 13 mat i 2 N range k1 t i 1 - el Y i 1 - el f Δ * eval break try : k2 t i 1 - el Δ 2 / + Y i 1 - el k1 2 / + f Δ * eval break try : k3 t i 1 - el Δ 2 / + Y i 1 - el k2 2 / + f Δ * eval break try : k4 t i 1 - el Δ + Y i 1 - el k3 + f Δ * eval break try : Y i el Y i 1 - el 16 / k1 2 k2 * + 2 k3 * + k4 + (* + : t i el t i 1 - el Δ + : 61 line for t Y augment 31 line :

ODE parametres β 3.5 12 0.25 41 mat :

Smath 1/1 Valery MCD ... ~ mwode113

t Y f t Y 1 el (β 1 el^^t Y 2 el (β 2 el^^+ t Y 1 el (β 3 el^^t Y 2 el (β 4 el^^+ 21 mat :

t0 0 : t2 2 : N 700 : 31 line

y0 0.5 121 mat :

[t2, N} govern Δ
reflexive wrt plots Δ t2 t0 - N / : 0.0029

increasing 'N' extends the blue raction valid up to some unknown limit. No mater 'N' within reason, the solver reconciliates itself @ 1.18

rkfΔ y0 t0 t2 N f Vrkfixed :

data UnestRow v data 1 el : j 2 data cols range v v data j el stack : for v 31 line :

ode113 y0 t0 t2 N t Y f mwode113 :

NOT zombie but doesn't 
reconciliate the project. ode113 rows 413

ode113 ode113 rows row 1.1771 2.5149 2.5242 13 mat

sol rkfΔ 1 col rkfΔ 2 col transpose UnestRow augment :

XY sol sol rows row (: 1.180000 3.604095 2.536967 13 mat

jnct XY 11 el XY 13 el 12 mat :

jnct jnct eval : 1.18 2.537 12 mat

plot sol 1 col sol 2 col augment sol 1 col sol 3 col augment jnct 31 sys :

plot

u sol 2 col : v sol 3 col : Δ u v - vectorize : 31 sys

Δf sol 1 col Δ augment :

Δf 1.18 5 1.18 1 - 2 2 mat 2 1 sys

'N' ... 700 1.1743 0.1237 - 1.1771 0.013 - 1.18 1.0671 32 mat

'N' ... 1000 1.176 0.0672 - 1.178 0.0471 1.18 2.6641 32 mat

Adaptive Quadrature ... yard stick.

Test function x f x 0.5 1 -^x - exp * :

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

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

lower limit
upper limit a 1026 -^: b 0.1 : 21 sys

Adaptive Carlos Algorithm

Initialise the while loop

sum 0 : i 1 : tol i el ε : aa i el a : 14 mat h i el b a - (2 / : 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 : 71 line

The adaptive while loop i 0 > fd aa i el h i el 2 / + f : fe aa i el 3 h i el * 2 / + f : 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 / (* : v1 aa i el : v2 fa i el : v3 fc i el : v4 fb i el : 14 mat v5 h i el : v6 ε : v7 s i el : v8 L i el : 14 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 : h i el v5 2 / : 16 mat tol i el v6 2 / : s i el s2 : L i el v8 1 + : i i 1 + : aa i el v1 : fa i el v2 : 15 mat fc i el fd : fb i el v3 : h i el h i 1 - el : s i el s1 : L i el L i 1 - el : 15 mat 41 line if if Solution sum : 91 line while

a b ʃ Solution :

Collected algorithm is 14 sec compared to open style  4 sec a b ʃ sum 0 : i 1 : tol i el ε : aa i el a : 14 mat h i el b a - (2 / : 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 : i 0 > fd aa i el h i el 2 / + f : fe aa i el 3 h i el * 2 / + f : 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 / (* : v1 aa i el : v2 fa i el : v3 fc i el : v4 fb i el : 14 mat v5 h i el : v6 ε : v7 s i el : v8 L i el : 14 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 : h i el v5 2 / : 16 mat tol i el v6 2 / : s i el s2 : L i el v8 1 + : i i 1 + : aa i el v1 : fa i el v2 : 15 mat fc i el fd : fb i el v3 : h i el h i 1 - el : s i el s1 : L i el L i 1 - el : 15 mat 41 line if if 91 line while sum 91 line eval :

π sqrt b sqrt erf * 0.61199136611177

t0 1 time :

x f x c * b 8 b 2 - 2 2 mat 2 1 sys

12½ D Adaptive Quadrature a b ʃ 0.611991366111614

Observe . . .

1. This infinite limit function is the most extreme encountered.

2. Best known integrator Mathcad/Matsoft . . . here exact !

3. Wolfram 0.611991 ... MCD 0.61199136611177

4. Infinite @ 1 is 'singular end point"

]]>

1 time t0 - 70 s *

Profitable Simpson-Romberg kernel ... large X integrand direction.

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

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

Test function x f x sqrt exp :

Adaptive Carlos Algorithm ... reduced eval(,)

Collected algorithm is 14 sec compared to open style  4 sec a b ʃ sum 0 : i 1 : tol i el ε : aa i el a : 14 mat h i el b a - (2 / : 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 : i 0 > fd aa i el h i el 2 / + f : fe aa i el 3 h i el * 2 / + f : 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 / (* : v1 aa i el : v2 fa i el : v3 fc i el : v4 fb i el : 14 mat v5 h i el : v6 ε : v7 s i el : v8 L i el : 14 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 : h i el v5 2 / : 16 mat tol i el v6 2 / : s i el s2 : L i el v8 1 + : i i 1 + : aa i el v1 : fa i el v2 : 15 mat fc i el fd : fb i el v3 : h i el h i 1 - el : s i el s1 : L i el L i 1 - el : 15 mat 41 line if if 91 line while sum 91 line eval :

t0 1 time :

Adaptive Quadrature a b ʃ 5.407528184656440

1 time t0 - 20.4 s *

5.40752818465627

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

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

7 true digits ... insufficient
integration subspaces. x f x 02 int 5.4075281

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

t0 1 time :

f 022000 ∫kernel 5.40752818465635

1 time t0 - 9.6 s *

The ODE Radovan ... Undocumented

t Y k D k 1 el Y 1 el * Y 2 el * k 2 el Y 2 el + / α k 1 el Y 1 el * Y 2 el * k 2 el Y 2 el + / * - 21 mat : IC 0.05 521 mat : k 0.375 105 -^21 mat : 13 mat

magnitude terminal blue [5,1.0101] α 1.0101 :

◄ slide system ► k 1 el

Undocumented AbsTol 107 -^: RelTol 108 -^: 12 mat

red
red
OK
red
? # # # # # dn_AdamsMoulton # # # # # dn_ExplicitRK45 # # # # # dn_GearsBDF # # # # # dn_ImplicitRK5 51 sys

sol IC 018100 t Y k D dn_GearsBDF : X sol 1 col : Y sol 2 col : Z sol 3 col : 13 mat x f1 X Y x cinterp : x f2 X Z x cinterp : 12 mat U x f1 x f2 - x 912 solve : 41 sys

X Y augment X Z augment 2 1 sys

Y Y rows row 5.000005 11 mat

Experiment configurations ... AbsTol 107 -^: RelTol 107 -^: 12 mat AbsTol 108 -^: RelTol 108 -^: 12 mat AbsTol 108 -^: RelTol 107 -^: 12 mat AbsTol 107 -^: RelTol 108 -^: 12 mat 41 sys

U U f1 12 mat U U f2 12 mat 21 mat 10.4454 2.5126 12 mat 10.4454 2.5126 12 mat 21 mat

t Y Clear 1

Difficult unique setup/solution

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

t y f 0.123 y * - 0.027 t * + 3.25 + 0.123 y * - 0.027 t * - 3.25 + / :

converts to standard form t Y D t Y 1 el f :

only '0' ic 0 :

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

Sensitive to step size 30/500 = 0.06 S ic 030500 t Y D Rkadapt :

S 15.5 50 15.5 25 - 2 2 mat 2 1 sys

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

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

Ω ic 030500 t Y D mwode15s :

S S rows row Ω Ω rows row 21 mat 30 42.2606 12 mat 15.48 22.6187 12 mat 21 mat

S rows Ω rows 21 sys 50125921 sys

Ω 15.5 50 15.5 25 - 2 2 mat 2 1 sys

jmG: 20230129

Rkadapt("1:vector", "2:number", "3:number", "4:number", "5:function") - (init, x1, x2, intvls, D) Uses the Runge-Kutta-Fehlberg Method RKF45

mwode15s("1:vector", "2:number", "3:number", "4:number", "5:function") - (init, x1, x2, intvls, D) solve stiff differential equations and DAEs, variable order method.