<p> plotG(vx,vy,a,b,c) </p>

0a4a6b7f-432c-4c55-b130-a01b60f9b054 71 15 12 decimal

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

Integrate a NON-stiff two variables function

y0 x0 x1 n f Euler X 1 el x0 : Y 1 el y0 : Δ x1 x0 - n / : 13 mat i 2 n range Y i el Y i 1 - el X i 1 - el Y i 1 - el f Δ * + eval : X i el X i 1 - el Δ + : 21 line for X Y augment 31 line :

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 :

x0 0 : x1 10 : y0 1 : n 400 : 41 line

x y f x sin y - :

euler y0 x0 x1 n f Euler :

rkf y0 x0 x1 n f rkf4 :

Euler & rkf4 2 variables integrators. "Integrators" in the sense of 'Integral Curve'

euler rkf 2 1 sys

Integrating a "widow" system

A widow system is one that does not contain the variable of integration, i.e: a system containing only relationship between derivatives. You must recast to suit Smath built-in.

init tspan n f Euler a tspan 1 el : b tspan 2 el : 12 mat h b a - (n / : t 1 el a : x 11 el init : i 1 n range t i 1 + el t i el h + : k h x 1 i el f * : x 1 i 1 + (el x 1 i el k + : 31 line for Unest x transpose ≡ data x transpose : S data 1 el : j 2 data rows range S S data j el stack : for t S augment 101 line :

Euler DE solver

1. Euler in programmed code is slow. 2. More profitable solvers have been designed since. 3. However, it belongs to the classroom. 4. Euler is "algo style". 5. Collapse in "Utilities"

x f x 2 el 1 x 1 el (2^- (x 2 el * x 1 el - 21 mat :

Euler algo style

init tspan n f Euler

Sol 2021 mat 01021 mat 200 f Euler :

t0 1 time :

sol_1 Sol 1 col Sol 2 col augment : sol_2 Sol 1 col Sol 3 col augment : 21 line

sol_1 sol_2 2 1 sys

t x D x f : recast the system for rkfixed ≡

sol 2021 mat 010100 t x D rkfixed :

s_1 sol 1 col sol 2 col augment : s_2 sol 1 col sol 3 col augment : 21 line

1 time t0 - 0.6 s *

s_1 s_2 2 1 sys

"one' revolution Van Der Pol

sol 2 col sol 3 col augment

Integrate a Jacobian 3 variables "stiff" system

x y z f x y * x 3^z * - x y * sin z * 21 mat :

f x y z JacobiStiff x y z f 1 el (x diff x y z f 1 el (y diff x y z f 1 el (z diff x y z f 2 el (x diff x y z f 2 el (y diff x y z f 2 el (z diff 23 mat :

Expand symbolic f x y z JacobiStiff y - 3 x 2^* z * + (- x x 3^- y x y * cos * z * x x y * cos * z * x y * sin 23 mat

Substitute by local definitions t Y J x t : y Y 1 el : z Y 2 el : f x y z JacobiStiff 41 line :

t Y J t t 2 Y 2 el * t Y 2 el x diff * + (* Y 1 el x diff - t Y 2 el * + (* - Y 1 el + t 1 t 2^Y 2 el y diff * - (* t t 2^- Y 1 el z diff + (* t Y 1 el x diff * Y 1 el + (t Y 1 el * cos * Y 2 el * t Y 1 el * sin Y 2 el x diff * + t t Y 1 el * cos * Y 2 el * t Y 1 el * sin Y 2 el y diff * + t Y 1 el z diff * t Y 1 el * cos * Y 2 el * t Y 1 el * sin + 23 mat

t0 0 : t1 3 : N 200 : 31 line

IC 5021 mat :

S IC t0 t1 N t Y J Rkadapt :

t S 1 col : U S 2 col : V S 3 col : 31 line

SOL_1 t U augment : SOL_2 t V augment : 21 line

SOL_1 SOL_2 2 1 sys

step i 1 U rows range U i el U i el 32 ≤ U i el 0 : cases for t U augment 21 line :

step SOL_2 2 1 sys <p> plotG(vx,vy,a,b,c) </p>

vx vy a b c plotG a char ≡ b size ≡ c color ≡ 13 mat n vx length : plot vx 1 el vy 1 el a b c augment : i 2 n range plot plot vx i el vy i el a b c augment stack : for plot 51 line :

Solving an "n order" DE is solving it as a "system" i.e: wrt its derivatives components in hierachical descending order. Observe this demo construction.

x y x 2 diff 3 x * x y x diff * + 7 x y * - 4 x * ≡

Y 2 el 4 x * 7 Y 1 el * + 3 x * Y 2 el * - 21 mat

x Y D Y 2 el x 43 Y 2 el * - (* 7 Y 1 el * + 21 mat :

x0 0 : start ≡ x1 2 : end ≡ N 200 : mesh ≡ 31 line

IC 1021 mat :

... implying user initial

Runge-Kutta-Fehlberg Method (RKF45). sol IC x0 x1 N var Deriv D Rkadapt :

x sol 1 col : X sol 2 col : Y sol 3 col : 31 line

U x X augment : V x Y augment : else X Y augment : 31 line

U V 2 1 sys else Y X augment 2 1 sys

Your next visit is about the RK5 "Modeller", explained below : why "Modeller". Observe that previously used 'x,t' have been evaporated by renaming 'xx, tt'.

The RK5 "acc" source: Smath examples

t y t 2 diff 2 t y t diff * + 10 t y * + t sin ≡

A 1210 t sin - 14 mat :

The source DE [NON homogeneous under forcing] t y t 2 diff 3 t * t y t diff * + 7 t y * - 4 t * ≡

Coefficients of the solved  NON-homogeneous DE A 13 tt * 7 - 4 tt * - 14 mat :

n expected number of solutions [Equation Power] n A cols 2 - :

i 2 : i n ≤ i i 1 + : f i 1 - el Ω i el : B i el Ω n i - 2 + el : 21 line for

n 2

fails
succeeds 1210 t sin - 14 mat 13 tt * 7 - 4 tt * - 14 mat 21 sys

B 0 Ω 2 el 21 mat

Construct the D(t,x) 
RK solving argument. B n 1 + el Ω 1 el : B n 2 + el 1 : f n el A B A 11 el / * (11 el - : 31 line

The 'f' RK  D(t,x)solving argument. f Ω 2 el 3 tt * Ω 2 el * - 7 Ω 1 el * + 4 tt * + 21 mat

<= conventionalise the RK solving argument tt Ω D f :

A Cauchy problem xx 1021 mat :

tt 1 el 0 : (start point ≡ T π : (segment length ≡ ε 104 -^: (accuracy ≡ h 0.25 : (step size ≡ 41 line

Метод Рунге-Кутта пятого порядка с автоматическим выбором шага

Metode Runge-Kutta orde lima dengan langkah adaptif

Fifth-order Runge–Kutta method with adaptive step accuracy ... adapted from Smath repository

Adaptive j tt 1 el h + : j T tt 1 el + h + ≤ j j h + : tx xx xx cols col : k1 h 3 / j tx D eval * : k2 h 3 / j h 3 / + tx k1 + D eval * : k3 h 3 / j h 3 / + tx k1 k2 + 2 / + D eval * : k4 h 3 / j h 2 / + tx 3 k1 * 9 k3 * + 8 / + D eval * : k5 h 3 / j h + tx 3 k1 * 9 k3 * - 2 / + 6 k4 * + D eval * : δ k1 92 / k3 * - 4 k4 * + 12 / k5 * - : δ min abs 5 ε * ≤ (δ max abs 5 ε * ≤ (& xx xx tx 12 / k1 4 k4 * + k5 + (* + augment : steps 1 xx cols el j : h h 2 * : 31 line j j h - : h h 2 / : 21 line if 81 line for steps transpose xx transpose augment 21 line :

<= ε & h are mutually compatible. Observe the number of values from tabulating "Adaptive".

assign to speed up Adaptive Adaptive :

Adaptive 1 col rows 43

idx Adaptive 1 col rows : 43

sub Adaptive 1 idx 13 submatrix : rk1 sub 1 col sub 2 col augment : rk2 sub 1 col sub 3 col augment : 31 sys

rk1 rk2 2 1 sys

======================= read more ... =================>>>>>

The parametric Rkadapt

t μ α D α 1 el μ 2 el μ 1 el - (* μ 1 el α 3 el μ 3 el - (* μ 1 el - μ 1 el μ 2 el * α 2 el μ 3 el * - 31 mat :

init 0.1 0.1 0.1 31 mat :

Example collaboration Smath Martin, Viacheslav, Radovan

t.0 0 :

t.max 50 :

n 2000 :

res init t.0 t.max n t μ 2 1.5 2031 mat D Rkadapt :

X res 1 col :

Y res 2 col :

Z res 3 col :

Y Z augment

res init t.0 t.max n t μ 21531 mat D Rkadapt :

X res 1 col :

Y res 2 col :

Z res 3 col :

Y Z augment

res 1012 -^1012 -^1012 -^31 mat t.0 t.max n t μ 21531 mat D Rkadapt :

x y axes x y / - :

X res 1 col :

Y res 2 col :

Z res 3 col :

Y Z augment axes 2 1 sys

===========The parametric option unary free fall under resistance

k 0.025 : g 9.807 : init 0 : par 0.125 : 41 line

t0 0 : t1 10 : N 25 : 41 line

u Φ P free Φ 1 el g - k P * Φ 1 el (2^P 2^+ sqrt * Φ 1 el 1 + (* - :

sol init t0 t1 N var Deriv par free Rkadapt :

X sol 1 col : Y sol 2 col : Ω X Y Y - Max / augment : x ∫ legacy from Maple ≡ x ∫ X Y Y - Max / x ainterp : 51 line

data X Ω 2 col . 10 black plotG :

x c 1 t0 x ≤ (x t1 ≤ (& cases :

inverse plot is more inductive

0 x ∫ x c * data 3 1 sys

10 ∫ 1 -

The state space (phase plot)

r 1.25 : (circle of radius 'r' ≡

xx δ D δ 2 el δ 1 el - 21 mat :

xo 0 : xe 2 π * : N 100 : 31 line

init 0 r 21 mat :

C init xo xe N var Deriv D rkfixed :

Plotting column 2 and 3 versus column 1 gives the functions versus x

C 1 col C 2 col augment C 1 col C 3 col augment 2 1 sys

Plotting column 3 versus column 2 gives a projection of the state space (phase plot). The red circle marks the initial state.

C 2 col C 3 col augment init 1 el init 2 el o 15 red 1 5 mat 2 1 sys

t Y Q 20 Y 1 el * 40 Y 2 el * + 2000 + 20 Y 1 el * 40 Y 2 el * - 2000 + 21 mat :

IC 0021 mat :

sol IC 0 0.1 100 var Deriv Q Rkadapt :

U sol 1 col sol 2 col augment : V sol 1 col sol 3 col augment : 21 line

U V 2 1 sys

This ODE solve memo

Y Clear 1

12 / y''' * - 14 / y'' * + y' - y + 0 ≡ (3rd order homogeneous/linear ≡

t Y Vect Y 2 el Y 3 el 214 / Y 3 el * - Y 2 el + Y 1 el - (* - 31 mat :

0 y 0 y' 0 y'' 31 mat 212 - 031 mat ≡

IC 212 - 031 mat :

soldemo IC 061000 t Y Vect dn_GearsBDF :

X soldemo 1 col : discretised variable of integration ≡ soldemo_1 X soldemo 2 col augment : function ≡ soldemo_2 X soldemo 3 col augment : 1rst derivative ≡ soldemo_3 X soldemo 4 col augment : 2nd derivative ≡ 41 line

soldemo_1 soldemo_2 soldemo_3 3 1 sys

k1 10.052 : k2 14.836 : k3 4.982 : 13 mat

Given a Pharmacokinetics DE system t X t diff k1 t X * - ≡ t Y t diff k1 t X * k2 t Y * - k3 t Z * + ≡ t Z t diff k2 t Y * k3 t Z * - ≡ 31 mat

t ω D k1 ω 1 el * - k1 ω 1 el * k2 ω 2 el * - k3 ω 3 el * + k2 ω 2 el * k3 ω 3 el * - 31 mat :

Φ 10031 mat 0130 t ω D rkfixed :

vx vy a b c plotG a char ≡ b size ≡ c color ≡ 13 mat n vx length : plot vx 1 el vy 1 el a b c augment : i 2 n range plot plot vx i el vy i el a b c augment stack : for plot 51 line :

vx Φ 1 col : X Φ 2 col : Y Φ 3 col : Z Φ 4 col : 41 line

Φ1 vx X o 5 black plotG : Φ2 vx Y o 5 black plotG : Φ3 vx Z o 5 black plotG : 31 line

Φ11 vx X augment : Φ22 vx Y augment : Φ33 vx Z augment : 31 line

Φ 0100 0.0333 0.7153 0.224 0.0607 0.0667 0.5117 0.3107 0.1776 0.1 0.366 0.3346 0.2993 0.1333 0.2618 0.3321 0.4061 0.1667 0.1873 0.32 0.4927 0.2 0.134 0.3061 0.5599 74 mat

Φ11 Φ22 Φ33 Φ1 Φ2 Φ3 61 sys

Φ11 Φ22 Φ33 Φ1 Φ2 Φ3 6 1 sys

The parametric DE solver

Lorentz DE system t L α D α 1 el L 2 el L 1 el - (* L 1 el α 3 el L 3 el - (* L 1 el - L 1 el L 2 el * α 2 el L 3 el * - 31 mat :

init 0.1 0.1 0.1 31 mat :

p 2 1.5 2031 mat :

t.0 0 :

t.max 50 :

n 2000 :

res init t.0 t.max n t L p D dn_GearsBDF :

x y axes x y / - :

X res 1 col :

Y res 2 col :

Z res 3 col :

Y Z augment axes 2 1 sys