Util

c97efdb4-31cf-48b0-b39e-2e6e902d4d66 9 4 5 false false 0 decimal Util

Euler Mascheroni constant γ 0.577215664901532860606 :

X Cols Ans# 0 : k# 1 k# range Ans# k# el Ans# k# col : for Ans# 31 line :

x# IsZero x# IsDefined x# 0 ≡ (& 11 line :

E# x# DIFF diff(E# x#_DIFF# EcuRep num2str,_DIFF#) concat str2num_DIFF# x# EcuRep 11 line :

E# x# n# DIFF diff(E# x#_DIFF# EcuRep num2str,_DIFF#, n# num2str) concat str2num_DIFF# x# EcuRep 11 line :

XY x cinterp X Y 12 mat XY 1 col XY 2 col 12 mat : U X 1 el UnitsOf Y 1 el UnitsOf 12 mat : U 2 el X U 1 el / Y U 2 el / x U 1 el / cinterp * X min eval x ≤ (* x X max eval ≤ (* 31 line :

Gamma function. Valid for Re(ζ)>0 ζ Γ G 0.99999999999999709182 57.156235665862923517 59.597960355475491248 - 14.136097974741747174 0.49191381609762019978 - 0.33994649984811888699 104 -^* 0.46523628927048575665 104 -^* 0.98374475304879564677 104 -^* - 0.15808870322491248884 103 -^* 0.21026444172410488319 103 -^* - 0.21743961811521264320 103 -^* 0.16431810653676389022 103 -^* - 0.84418223983852743293 104 -^* 0.26190838401581408670 104 -^* - 0.36899182659531622704 105 -^* 151 mat : zg ζ 543128 / + : κ Clear 12 mat 2 π * sqrt zg 0.5 ζ 0.5 - (*^(2^* e zg^/ G 1 el G κ el ζ κ + 2 - / κ 215 sum + (* eval 31 line :

Runge-Kutta 2nd and 3rd order. 2 D to uo N E RK23 n uo length : u uo 1 n range el : U u : t to 1 el : T t 11 mat : v 1 n range : 16 mat ho to 2 el to 1 el - N / eval : h ho : n 1 : 13 mat t u k K h t u D * eval h k * eval try 11 line : t to 2 el < k.1 t u u K : k.2 t h + u k.1 + k.1 K : k.3 t 0.5 h * + u 0.25 k.1 * + 0.25 k.2 * + k.2 K : d k.1 2 k.3 * - k.2 + 3 / normi : (q E u normi ho 12 mat max * : (≤ t t h + : T n n 1 + : el t : 12 mat u u k.1 4 k.3 * + k.2 + 6 / + eval : U v n el u v el : 12 mat 21 line if d 0 ≠ h ho 0.9 h * q d / 3 nthroot * 12 mat min eval : if 51 line while T U transpose augment eval 51 line :

Numerical Inverse Laplace Transform

Based on algorithm of Oscar: https://en.smath.info/forum/yaf_postst1341_Signals-and-Systems.aspx

a ns nd ILT.αβ m ns nd + 1 + : α m 1 matrix : β α : γ nd 1 matrix : 14 mat μ 1 m range : nδ 2 nd range : nσ ns 3 + (m range : 13 mat α μ el a μ 1 - (π * i * + eval : β μ el 1 - (μ 1 +^e a^* eval : 12 mat γ 1 el 2 nd -^: γ nδ el nd 1 + Γ 2 nd^nd 2 + nδ - Γ * nδ Γ * / γ nδ 1 - el + eval : 12 mat β 1 el 0.5 β 1 el * : β nσ el β nσ el γ m nσ - 1 + el * eval : 12 mat α β μ 13 mat 61 line :

Inverse Laplace transform of F(s) 1 F t.i t.e N ILT t N 1 matrix : ft t : Fs αL rows 1 matrix : 13 mat t.o 105 -^t.e t.i - UnitsOf * t.i 12 mat max eval : Δt t.e t.o - N 1 - / : 12 mat n 1 N range t n el t.o n 1 - (Δt * + : Fs μL el βL μL el t n el / αL μL el t n el / F * eval : 12 mat ft n el Fs sum Re eval : 21 line for t ft augment eval 41 line :

α and β values αL βL μL 13 mat 62019 ILT.αβ eval :

equrep

E# a# b# EcuRep equrep(E# num2str, a# num2str, b# num2str) concat str2num 11 line :

E# Rules# EcuRep ans# E# : n# 1 Rules# rows range ans# ans# Rules# n# 1 el Rules# n# 2 el EcuRep : for ans# 31 line :

E# Rules# EcuRepN ans# E# : N# 1 EcuRepN range old# ans# : num2str ans# ans# Rules# EcuRep : num2str ≡ break continue if 11 line for ans# 31 line :

E# equal2EQ ans# E# num2str : ans# \2261\ findstr 1 - ≡ ans# ans# \2261\0 concat : n# ans# \2261\ findstr 1 el : 12 mat n# ans# \2261\ findstr 1 el : if ans# 1 n# 1 - substr str2num ans# n# 1 + ans# strlen n# - substr str2num 12 mat 31 line :

E# x# depends E# x# x#_DEPENDS.X# EcuRep num2str x#_DEPENDS.X# num2str findstr 1 - ≠ 11 line :

EcuRepN 20 :

isol

EQ EQ(-a,b) EQ(a,0-b) minus 13 mat :

EQ EQ EQ(a^b,z) EQ.B(a,b,z,a^b,ln(z)/ln(a),z^(1/b)) power 13 mat stack :

Not usual but better EQ.B(a,b,z,a^b,ln(z)/ln(a),exp(ln(z)/b))

EQ EQ EQ(a/b,z) EQ.B(a,b,z,a/b,z/a,z*b) quotient 13 mat stack :

EQ EQ EQ(a*b,z) EQ.B(a,b,z,a*b,z/a,z/b) times 13 mat stack :

EQ EQ EQ(a+b,z) EQ.B(a,b,z,a+b,z-a,z-b) plus 13 mat stack :

EQA sin cos tan sec csc cot sinh cosh tanh sech csch coth exp_EQ#.10_EQ#.2 erf asin acos atan asec acsc acot asinh acosh atanh asech acsch acoth ln log10 sqrt erfinv 216 mat transpose :

EQA rows 16

k 1 EQA rows range EQ EQ EQ(EQA k 1 el (a),b) concat EQ.A(a,b, EQA k 1 el (a), EQA k 2 el (b) concat EQA k 1 el 13 mat stack : for

k 1 EQA rows range EQ EQ EQ(EQA k 2 el (a),b) concat EQ.A(a,b, EQA k 2 el (a), EQA k 1 el (b) concat EQA k 2 el 13 mat stack : for

E# x# isolate a#_EQ#.10 10 a#^11 line : a#_EQ#.2 a# 2^11 line : 12 mat a# b# f# i# EQ.A a# x# depends a# i# EQ f# b# EQ if 11 line : a# b# z# f# i# j# EQ.B a# x# depends b# x# depends k# f# x# 0 EcuRep : f# k# - z# k# - EQ 21 line a# j# EQ if b# x# depends b# i# EQ z# f# EQ if if 11 line : ans# E# equal2EQ : A# ans# 1 el numden : B# ans# 2 el numden : 13 mat ans# A# 1 el B# 2 el * A# 2 el B# 1 el * - 0 EQ EQ EcuRepN : a# b# EQ a# b# ≡ (: ans# 71 line :

E# x# isol ans# E# x# isolate equal2EQ : ans# 1 el num2str x# num2str ≡ (ans# 2 el x# depends ¬ & ans# 2 el ans# 1 el ans# 2 el - if 21 line :

laplace

LT LT(a) LT.0(a) 013 mat :

LT LT LT(LT(a)) LT(a) LTLT 13 mat stack :

LT LT LT(sin(a)) LT.φω(sin(a),a,(s*sin(lφ)+lω*cos(lφ))/(s^2+lω^2) sin 13 mat stack :

LT LT LT(cos(a)) LT.φω(cos(a),a,(s*cos(lφ)-lω*sin(lφ))/(s^2+lω^2) cos 13 mat stack :

LT LT LT(sinh(a)) LT.φω(sinh(a),a,(s*sinh(lφ)+lω*cosh(lφ))/(s^2-lω^2) sinh 13 mat stack :

LT LT LT(cosh(a)) LT.φω(cosh(a),a,(s*cosh(lφ)+lω*sinh(lφ))/(s^2-lω^2) cosh 13 mat stack :

LT LT LT(exp(a)) LT.φω(exp(a),a,exp(lφ)/(s-lω) exp 13 mat stack :

LT LT LT(ln(a)) LT.φω(ln(a),a,(ln(lω)-ln(s)-'γ)/s ln 13 mat stack :

LT.α sin cos sinh cosh 14 mat :

lk 14 range LT LT LT(c*exp(a)* LT.α lk el (b)) concat LT.α(exp(a), LT.α lk el (b),a,b,c) concat exp* LT.α lk el concat 13 mat stack : for

LT LT LT(-a) 0-LT(a) minus 13 mat stack :

LT LT LT(a^n*b) LT.n(a,n,b) nshift 13 mat stack :

LT LT LT(a^b) LT.pot(a,b) pot 13 mat stack :

LT LT LT(a*b) LT.prod(a,b) product 13 mat stack :

LT LT LT(a*b) LT.prod(b,a) product 13 mat stack :

LT LT LT(a*(b+c)) LT(a*b)+LT(a*c) expand* 13 mat stack :

LT LT LT((a+b)/c) LT(a/c)+LT(b/c) expand/ 13 mat stack :

LT LT LT(a+b) LT(a)+LT(b) plus 13 mat stack :

LT.O t s^(-2) 12 mat :

LT.O LT.O x(t) X(s) 12 mat stack :

LT.O LT.O u(t) 1/s 12 mat stack :

LT.O LT.O δ(t) 112 mat stack :

LT.O LT.O diff(x(t),t) s*X(s)-x.o 12 mat stack :

LT.O LT.O diff(diff(x(t),t),t) s^2*X(s)-x'.o-s*x.o 12 mat stack :

LT.O LT.O diff(diff(diff(x(t),t),t),t) s^3*X(s)-x''.o-s*x'.o-s^2*x.o 12 mat stack :

E# lt t s t x s X lω lφ Clear a# IsT# a# num2str t ≡ 11 line : a# LT.0 a# t depends row# LT.O a# num2str 1 findrows : (0 ≡ a# LT row# 2 el str2num if a# s 1 -^* if 11 line : A# B# a# b# c# LT.α lω b# t DIFF : lα a# t DIFF : IsLT? 13 mat lα t depends lω t depends | (c# t depends | c# A# * B# * LT c# B# lt * s s lα - EcuRep if 21 line : a# n# b# LT.n n# t depends ¬ (a# IsT# & 1 - (n#^b# LT s n# DIFF * a# n#^b# * LT if 11 line : f# a# tr# LT.φω lω a# t DIFF : lφ a# lω t * - : 12 mat a# t depends lω t depends f# LT tr# if a# s 1 -^* if 21 line : a# n# LT.pot n# IsDefined a# IsT# & n#! s n# 1 + (-^* a# n#^LT if 11 line : a# b# LT.prod a# t depends ¬ a# b# LT * a# IsT# b# LT s DIFF - a# b# * LT if if 11 line : ans# E# sqrt(a) a^(1/2) nthroot(a,b) a^(1/b) 22 mat EcuRep : E# num2str \2261\ findstr 1 - ≡ E# LT LT EcuRepN A# B# 12 mat E# equal2EQ : A# LT LT EcuRepN B# LT LT EcuRepN ≡ 21 line if 101 line :

Laplace ODE examples

Solving ODE's by Laplace Transform

Shorthands

plot x x.lt x.rk 1 col x.rk 2 col . 10 green augment 21 sys :

N 25 :

ode a t x t 2 diff * b t x t diff * + c t x * + t φ ≡ :

t.end RK t x D 1 a / a x 2 el * t φ c x 1 el * - b x 2 el * - 21 mat * : t x D 0 t.end 12 mat x.o x'.o 12 mat N 103 -^RK23 21 line :

t φ 1 t - exp 2 / 3 t * cos * - :

a b c 13 mat 21313 mat :

x.o x'.o 12 mat 0112 mat :

t.end 6 :

s X.lt ode lt s X isol : x.lt s X.lt 0 t.end N ILT : t x.lt x.lt t cinterp : x.rk t.end RK : 41 sys

plot

Where

ode 2 t x t diff t diff * 3 t x * + t x t diff + 2 t - exp 3 t * cos * - 2 / ≡

ode lt s X 3 s + (* 21 - s 2^s X * + (* + 291 s + (2^+ (* 1 s + (s * - 291 s + (2^+ (* s * / ≡

ode lt s X isolate s X 291 s + (2^+ (* 12 s * + (* 1 s + (s * - 291 s + (2^+ (* s * 3 s 12 s * + (* + (* / ≡

t φ t 2^t sin + :

a b c 13 mat 10213 mat :

x.o x'.o 12 mat 1 - 0.5 12 mat :

t.end 5 :

s X.lt ode lt s X isol : x.lt s X.lt 0 t.end N ILT : t x.lt x.lt t cinterp : x.rk t.end RK : 41 sys

plot

Where

ode 2 t x * t x t diff t diff + t 2^t sin + ≡

ode lt 1 - 2 s 1 s s X * + (* 2 s X * + (* + 2 / s 3^21 s 2^+ (* + s 3^1 s 2^+ (* / ≡

ode lt s X isolate s X 1 s 2^+ (s 3^* 12 s * - (* 2 s 3^21 s 2^+ (* + (* + 21 s 2^+ (* s 3^* 2 s 2^+ (* / ≡

t φ 0.5 t * - exp 3 t * 2 - sin * :

a b c 13 mat 2 - 13 - 13 mat :

x.o x'.o 12 mat 0012 mat :

t.end 10 :

s X.lt ode lt s X isol : x.lt s X.lt 0 t.end N ILT : t x.lt x.lt t cinterp : x.rk t.end RK : 41 sys

plot

Where

ode 2 t x t diff t diff * 3 t x * + t x t diff - (- t 2 / - exp 2 - 3 t * + sin * ≡

ode lt s X 3 s 1 - 2 s * + (* + (* - 212 s * + (2 sin * - 62 cos * + (* 3612 s * + (2^+ / ≡

ode lt s X isolate s X 212 s * + (2 sin * - 62 cos * + (* 3612 s * + (2^+ (3 s 1 - 2 s * + (* + (* / - ≡

t φ t 2^t 5 / cosh * :

a b c 13 mat 52113 mat :

x.o x'.o 12 mat 1312 mat :

t.end 5 :

s X.lt ode lt s X isol : x.lt s X.lt 0 t.end N ILT : t x.lt x.lt t cinterp : x.rk t.end RK : 41 sys

plot

Where

ode 5 t x t diff t diff * 2 t x t diff * + t x + t 2^t 5 / cosh * ≡

ode lt 21 - s s X * + (* 53 - s 1 - s s X * + (* + (* + s X + 2505 s * 0 cosh * 0 sinh + (1 - 25 s 2^* + (* 101 - 25 s 2^* + (0 cosh * 105 s * 0 cosh * 0 sinh + (* s * - (* s * + (* 1 - 25 s 2^* + (3^/ - ≡

ode lt s X isolate s X 21 - 25 s 2^* + (3^1255 s * 0 cosh * 0 sinh + (1 - 25 s 2^* + (* 101 - 25 s 2^* + (0 cosh * 105 s * 0 cosh * 0 sinh + (* s * - (* s * + (* - (* 53 s + (* 1 - 25 s 2^* + (3^* + 1 - 25 s 2^* + (3^1 s 25 s * + (* + (* / ≡

Alvaro