<content> equrep </content>

56dc514e-4886-49ce-bef7-36cf44e85c39 397 4 5 false false 0 0 decimal

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

f.feval# x.feval# feval f.feval# num2str (x.feval# num2str) concat str2num 11 line :

ffeval xfeval yfeval feval ffeval num2str (xfeval num2str, yfeval num2str) concat str2num 11 line :

E\032E\ a\032E\ b\032E\ EquRep equrep(E\032E\ num2str, a\032E\ num2str, b\032E\ num2str) concat str2num 11 line :

E\032E\ a\032E\ b\032E\ q\032E\ EquRep q\032E\ 1 EqualQ E\032E\ a\032E\ b\032E\ EquRep sa\0324\ a\032E\ num2str : sq\0324\ q\032E\ num2str : 12 mat str\0324\ equrep(E\032E\ num2str, sa\0324\, concat : sq\0324\ strlen 9 > (sq\0324\ 19 substr Assuming(≡ (& str\0324\ str\0324\ sq\0324\ 1 sq\0324\ strlen 1 - substr concat : str\0324\ str\0324\, b\032E\ num2str,_sa\0324\ concat : 21 line str\0324\ str\0324\ if(sq\0324\, b\032E\ num2str, sa\0324\ concat : if str\0324\)) concat str2num 41 line if 11 line :

E\032E\ Rules\032E\ EquRep ans\0324\ E\032E\ : Rules\032E\ cols 3 < n\0324\ 1 Rules\032E\ rows range ans\0324\ ans\0324\ Rules\032E\ n\0324\ 1 el Rules\032E\ n\0324\ 2 el EquRep : for n\0324\ 1 Rules\032E\ rows range ans\0324\ ans\0324\ Rules\032E\ n\0324\ 1 el Rules\032E\ n\0324\ 2 el Rules\032E\ n\0324\ 3 el EquRep : for if ans\0324\ 31 line :

E\032E\ Rules\032E\ EquRepN ans\0324\ E\032E\ : N\0324\ 1 EquRepN range ans\0324\ num2str ans\0324\ ans\0324\ Rules\032E\ EquRep : num2str ≡ break continue if 11 line for ans\0324\ 31 line :

EquRepN 20 :

a\032E\ b\032E\ EqualQ assume IsDefined a\032E\ num2str b\032E\ num2str ≡ a\032E\ num2str b\032E\ num2str ≡ if 11 line :

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

S\032E\ Sys2EQ S\0324\ 0 : k\0324\ 0 : 12 mat s\0324\ S\032E\ sys2mat.1 s\0324\ s\0324\ equal2EQ : S\0324\ k\0324\ k\0324\ 1 + : el s\0324\ 1 el s\0324\ 2 el - : 21 line for S\0324\ 31 line :

E\032E\ x\032E\ depends E\032E\ x\032E\ x\032E\_DEPENDS.X# EquRep num2str x\032E\_DEPENDS.X# num2str findstr 1 - ≠ 11 line :

a\032E\ b\032E\ x\032E\ depends a\032E\ x\032E\ depends b\032E\ x\032E\ depends & 11 line :

E\032E\ x\032E\ DIFF diff(E\032E\ x\032E\_DIFF# EquRep num2str,_DIFF#) concat str2num_DIFF# x\032E\ EquRep 11 line :

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

Replacing deep EquRepN 20 :

Runge-Kutta 2nd and 3rd order. 2 DS\032E\ to\032E\ xo\032E\ N\032E\ E\032E\ L.ode nx\0324\ xo\032E\ length eval : c\0324\ 1 nx\0324\ range : n\0324\ 1 N\032E\ 1 + (range : 13 mat ut\0324\ to\032E\ 1 el to\032E\ 2 el - UnitsOf num2str : ux\0324\ 0 : ux\0324\ c\0324\ el xo\032E\ c\0324\ el UnitsOf num2str : 12 mat m kg s K A cd 16 mat 11111116 mat : E\032E\ IsString U\0324\ xo\032E\ c\0324\ el to\032E\ 1 el to\032E\ 2 el N\032E\ XX# YY# DS\032E\ Rkadapt : t\0324\ U\0324\ 1 col : X\0324\ U\0324\ n\0324\ 1 c\0324\ + el : 12 mat 21 line t\0324\ to\032E\ min eval : t2\0324\ to\032E\ max eval : h\0324\ t2\0324\ t\0324\ - N\032E\ / : 13 mat k\0324\ 1 : T\0324\ t\0324\ 11 mat : ho\0324\ h\0324\ : 13 mat x\0324\ xo\032E\ c\0324\ el : X\0324\ x\0324\ transpose : 12 mat t\0324\ x\0324\ err\0324\ h\0324\ K\0324\ h\0324\ t\0324\ x\0324\ DS\032E\ * eval err\0324\ try 11 line : E\032E\ 0 ≡ T\0324\ n\0324\ el t\0324\ n\0324\ 1 - (h\0324\ * + : k\0324\ 2 N\032E\ 1 + (range k\0324\.1 T\0324\ k\0324\ el x\0324\ x\0324\ h\0324\ K\0324\ : k\0324\.2 T\0324\ k\0324\ el h\0324\ + x\0324\ k\0324\.1 + k\0324\.1 h\0324\ K\0324\ : k\0324\.3 T\0324\ k\0324\ el 0.5 h\0324\ * + x\0324\ 0.25 k\0324\.1 * + 0.25 k\0324\.2 * + k\0324\.2 h\0324\ K\0324\ : x\0324\ x\0324\ k\0324\.1 4 k\0324\.3 * + k\0324\.2 + 6 / + eval : X\0324\ k\0324\ c\0324\ el x\0324\ c\0324\ el : 12 mat 41 line for 21 line t\0324\ t2\0324\ ≤ k\0324\.1 t\0324\ x\0324\ x\0324\ h\0324\ K\0324\ : k\0324\.2 t\0324\ h\0324\ + x\0324\ k\0324\.1 + k\0324\.1 h\0324\ K\0324\ : k\0324\.3 t\0324\ 0.5 h\0324\ * + x\0324\ 0.25 k\0324\.1 * + 0.25 k\0324\.2 * + k\0324\.2 h\0324\ K\0324\ : d\0324\ k\0324\.1 2 k\0324\.3 * - k\0324\.2 + 3 / normi eval : q\0324\ E\032E\ x\0324\ normi ho\0324\ 12 mat max * : d\0324\ q\0324\ ≤ t\0324\ t\0324\ h\0324\ + : T\0324\ k\0324\ k\0324\ 1 + : el t\0324\ : 12 mat x\0324\ x\0324\ k\0324\.1 4 k\0324\.3 * + k\0324\.2 + 6 / + eval : X\0324\ k\0324\ c\0324\ el x\0324\ c\0324\ el : 12 mat 21 line if d\0324\ 0 ≠ h\0324\ ho\0324\ 0.9 h\0324\ * q\0324\ d\0324\ / 3 nthroot * 12 mat min eval : if 71 line while if 51 line if m kg s K A cd Clear E\032E\ IsString E\032E\ 0 ≡ (| x\0324\ n\0324\ c\0324\ el X\0324\ n\0324\ c\0324\ el ux\0324\ c\0324\ el str2num * : t\0324\ n\0324\ el T\0324\ 1 el n\0324\ 1 - (ho\0324\ * + : x\0324\ n\0324\ c\0324\ el T\0324\ X\0324\ c\0324\ col t\0324\ n\0324\ el cinterp ux\0324\ c\0324\ el * : 21 line if t\0324\ ut\0324\ str2num * x\0324\ augment eval 81 line :

<content> pMesh </content>

az el pView α az : ε el : 12 mat ε ε 180 ° * + 360 ° * mod 360 ° * + 360 ° * mod 180 ° * - : ε 90 ° * > ε 180 ° * ε - : α α 180 ° * + : 12 mat ε 90 ° * - < ε 180 ° * - ε - : α α 180 ° * + : 12 mat 0 if if α α 360 ° * mod 360 ° * + 360 ° * mod : α cos ε sin α sin * - ε cos α sin * α sin ε sin α cos * ε cos α cos * - 0 ε cos ε sin 33 mat eval 51 line :

az el pView2 az el pView 13 range 12 range el 11 line :

P cyl2xyz P 1 el P 2 el cos * P 1 el P 2 el sin * P 3 el 31 mat :

P ρ φ z 31 mat ≡

x# IsMatrix x# length 0 > eval 0 try 11 line 10 if 11 line :

a b 12 mat GetType unknown

a b 12 mat IsMatrix 1

range utility a b n pR a b a - n / 0 n range * + eval 11 line :

PMesh returns a matrix suitable for 3D plots. F# F'# uo vo pTMesh F F# F# F# cols col augment eval : F' F'# eval : 12 mat nu F rows 1 - : nv F cols 1 - : M 03 matrix eval : r' 0 : 14 mat np uo length : p' np 1 + : r.np 1 np range : r.nv 1 nv 1 - (range : c 13 range : 15 mat u 1 nu range v 1 nv range r M rows : k r r.np + : r' r' 1 + : 13 mat M k c el F uo k r - el u + vo k r - el v + el c el : v nv ≡ M r p' + c el M r np + c el : M r p' + c el F' r' el c el : if 31 line for r M rows : k r r.nv + : 12 mat M k c el M u 1 - (p' nv * nv + 1 - (* p' v 1 - (* + k r - 1 - (p' * - 1 - c el : 31 line for M 51 line :

FS U V pTGrid A pr v A 2 el A 1 el - 2 / A + : v v length el A v length el : eval 21 line : FS U pr V pr pGrid 21 line :

pGrid(F,U,V) returns a length(U)*lenght(V) matrix with [u,v,F(u,v)] FS# X# Y# pGrid r# 1 X# length range : c# 1 Y# length range : M# 0 : 13 mat FS# X# 1 el Y# 1 el feval IsMatrix num2str 1 ≡ eval M# r# c# el FS# X# r# el Y# c# el feval : M# r# c# el X# r# el Y# c# el FS# X# r# el Y# c# el feval 31 mat : if eval 21 line :

pGrid(M) converts M = pGrid(F,U,V) into three matrices [X,Y,Z]with the same shape than matlab's meshgrid function M# pGrid r# 1 X# rows range : c# 1 Y# cols range : X# 0 : Y# 0 : Z# 0 : 15 mat X# r# c# el M# r# c# el 1 el : Y# r# c# el M# r# c# el 2 el : Z# r# c# el M# r# c# el 3 el : 13 mat eval 21 line :

PMesh returns a maitrx suitable for 3D plots. F# uo vo pMesh F F# eval : nu\032E\ F rows 1 - : nv F cols 1 - : 13 mat np uo length : r.np 1 np range : r.nv 1 nv 1 - (range : 13 mat M 03 matrix eval : c 13 range : 12 mat u\032E\ 1 nu\032E\ range v 1 nv range r M rows : k r r.np + : 12 mat M k c el F uo k r - el u\032E\ + vo k r - el v + el c el : 21 line for r M rows : j r r.nv + : 12 mat M j c el M u\032E\ 1 - (np nv * nv + 1 - (* np nv 1 - (* + j r - 1 - (np * - c el : 31 line for M 51 line :

S z pZoom z length 2 ≡ S 1 col z 1 el * S 2 col z 2 el * augment eval S 1 col z 1 el * S 2 col z 2 el * S 3 col z 3 el * augment eval if 11 line :

P cl ls w cs pLine cs 0 ≡ polyline P 11 mat cl ls w stack polygon P 11 mat cl ls w cs stack if 11 line :

S# CM γ2 pMeshShow S S# eval : k 2 S rows (range : cm S 3 col CM pClr : 13 mat γ γ2 13 range 12 range el : S2 S γ * eval : 12 mat L k el S2 k row S2 k 1 - row stack cm k el solid 10 pLine : 11 mat eval 31 line :

M# G pClr M M# eval : C 0 : r 1 M length range : 13 mat G 0 ≡ M length 1 matrix G IsString C r el G : a b n 13 mat M min M max G length 13 mat eval : 11 mat C r el G M r el a - (n 1 - b a - / (* 0 round 1 + el eval : 21 line if if 21 line :

n a pCmapJet r 1 n range : m 0.25 n * trunc : M 0 : 13 mat v z v rows 1 matrix : o 255 m 1 matrix + : 12 mat x 2551 m range m / * trunc vectorize : y 255 m / m 2 / m range (* trunc vectorize : 12 mat R y z o z x o y reverse stack : G y z x o x reverse y z stack : B y o x reverse o z y z stack : k R length n - : M r el B r k + el 256 G r k + el * + 2562^R r k + el * + 2563^255 a * trunc * + : M r el M r el 232^2 M r el * 232^1 - ≥ (* - eval : 91 line :

CM 32 0.8 pCmapJet :

Matlab default view γ 37.5 ° * - 30 ° * pView :

∞ 1010^:

CycleColors 5 :

X X'' X# Y# N# sf# L.VF s1 G({X},{X}') s2n num2str : s2 eval(mat({X}',{X''},2,1)) s2n num2str : s1 : s2 concat trace s1 : s2 concat str2num G X# Y# N# sf# L.VF 51 line :

DS# X# Y# N# sf# L.VF B# X# Y# 10 % * 10 % * pBox eval : DS# num2str B# N# sf# pVField X# Y# augment X# 1 el Y# 1 el o augment 31 sys 21 line :

X# Y# px# py# pBox B# X# min X# max Y# min Y# max 22 mat eval : px# B# 1 el abs * - px# B# 2 el abs * py# B# 3 el abs * - py# B# 4 el abs * 22 mat B# + eval 21 line :

xVF xdVF eqVF X# Y# N# sf# L.VF s1 G({xVF},{xdVF}) s2n num2str : s2 eval(mat({xdVF},{eqVF},2,1)) s2n num2str : s1 : s2 concat str2num s1 : s2 concat trace G X# Y# 10 % * 10 % * pBox N# sf# pVField X# Y# augment X# 1 el Y# 1 el o augment 31 sys 51 line :

Fs B N sf L pVField m kg s K A cd 16 mat 11111116 mat : x1 x2 y1 y2 22 mat B eval : nx N 1 el : ny N 2 el : 13 mat M 0 : V 0 : P 0 : X 0 : Y 0 : R 0 : 16 mat Δx Δy 12 mat x2 x1 - nx / y2 y1 - ny / 12 mat eval : r c 12 mat 1 ny range 1 nx range 12 mat : X r c el c 0.5 - (Δx * x1 + : Y r c el r 0.5 - (Δy * y1 + : 12 mat V r c el Fs X r c el Y r c el feval eval X r c el Y r c el 12 mat try : V V Δy V norme vectorize max / * eval : sf IsString a 0.5 0.3 0.15 0.65 0.35 14 mat * eval : a 0.5 0.3 0.15 0.65 0.35 14 mat * sf * eval : if ao L 1 el a 1 el - 000 a 1 el 15 mat * transpose eval : bo L 2 el a 3 el a 4 el - a 3 el a 4 el + a 3 el a 4 el + (- a 3 el a 4 el + a 3 el a 4 el - 15 mat * transpose eval : k 1 nx ny * (range t V k el norme eval : sf IsString x y 12 mat ao bo 12 mat sf t feval * eval : x y 12 mat ao bo 12 mat : if t 0 ≡ s c 12 mat 0112 mat : θ V k el 1 el sign V k el 2 el t / acos * - : s c 12 mat θ sin θ cos 12 mat eval : 12 mat if j 15 range x j el y j el 21 mat X k el Y k el 21 mat c x j el * s y j el * - s x j el * c y j el * + 21 mat + : for M R R 1 + : el x y augment eval : iter 1 CycleColors range M R R 1 + : el ∞ : for 61 line for M mat2sys.1 eval 131 line :

Fs B N sf pVField Fs B N sf 1112 mat pVField 11 line :

Quiver plot of f(x,y)≡[f.1(x,y) f.2(x,y)] FS# X# Y# sf# pQuiver m kg s K A cd 16 mat 11111116 mat : r# 1 X# length range : c# 1 Y# length range : 12 mat Z# 0 : A# 0 : M# 0 : n# 0 : sf2# 1 : 15 mat arr# 01 0.8 1 0.8 15 mat 00 0.08 0 0.08 - 15 mat i * + transpose eval : O5# 151 matrix + eval : FS# IsString Z# r# c# el FS# X# r# el Y# c# el feval : Z# FS# eval : if sf# IsString Z# r# c# el sf# Z# r# c# el norme feval Z# r# c# el * eval : sf2# sf# : if scale# 0.9 Y# 2 el Y# 1 el - (* sf2# * Z# norme vectorize max / eval : A# r# c# el scale# Z# r# c# el 1 el i Z# r# c# el 2 el * + (* arr# * X# r# el Y# c# el i * + (O5# * + eval : k 1 A# length range a# A# k el : M# n# n# 1 + : el a# Re vectorize a# Im vectorize augment eval : iter 1 CycleColors range M# n# n# 1 + : el ∞ : for 31 line for M# mat2sys.1 eval 111 line :

ic to tf N 2 D pRKFam S IC t u Clear k 1 ic length eval range : 12 mat S k el ic k el to tf N t u D rkfixed : IC k el ic k el transpose o augment : mat2sys.1 S k el S k el 2 col S k el 3 col augment : mat2sys.1 21 sys 31 line :

Plots a family of rk by various ic

FamilyPlot(f(α,x),α,x) plots f for the ranges α and x. 2 f a b c d M N pFamily x# y# Clear x a b M LinSpace : y c d N LinSpace : 13 mat x# y# f x y pFamily 21 line :

FamilyPlot(f(α,x),α,x) plots f for the ranges α and x. 2 f α x pFamily M α length : N x length : P 0 : F 0 : 14 mat c 1 M range r 1 N range F r c el α c el x r el f Re eval 0 try : for for c 1 M range P c el x F c col augment eval : for P mat2sys.1 eval 41 line :

1 f 1 g T N pParametric t T 1 el T 2 el N LinSpace : t f vectorize t g vectorize 12 mat 21 line :

1 f T N pParametric t T 1 el T 2 el N LinSpace : X N 1 matrix : Y X : F t f vectorize : 13 mat k 1 N range f F k el : X k el f 1 el : Y k el f 2 el : 13 mat for X Y augment 41 line :

x Round x Re 6 round i x Im 6 round * + :

Derivative utility X\032E\ t\032E\ n\032E\ p\032E\ DIFF R\0324\ 0 : k\0324\ 0 : 12 mat x\0324\ X\032E\ str\0324\ diff(x\0324\ num2str, t\032E\ num2str, n\032E\ num2str) concat : R\0324\ k\0324\ k\0324\ 1 + : el str\0324\ str2num p\032E\ * : 21 line for k\0324\ 1 ≡ R\0324\ 1 el n\0324\ 1 k\0324\ range : A\0324\ X\032E\ : 12 mat A\0324\ n\0324\ el R\0324\ n\0324\ el : 21 line if 31 line :

∂(x,t) returns x'*∂t x\032E\ t\032E\ ∂ U\0324\ t\032E\ Unknowns : R\0324\ 0 : D\0324\ x\032E\ : 13 mat k\0324\ 1 : k\0324\ U\0324\ length ≤ k\0324\ k\0324\ 1 + : str\0324\ U\0324\ k\0324\ el num2str : Unknowns(dSymbol str\0324\) concat str2num num2str 0 ≡ 0 l\0324\ 05 range t\032E\ U\0324\ k\0324\ el l\0324\ 1 + 1 DIFF num2str 0 ≡ break continue if for p\0324\ dSymbol str\0324\ concat dSymbol dSymbol concat dSymbol strrep str2num : n\0324\ str\0324\ dSymbol findstr : (1 - ≡ (R\0324\ k\0324\ 1 el R\0324\ k\0324\ 2 el R\0324\ k\0324\ 3 el 13 mat U\0324\ k\0324\ el l\0324\ p\0324\ 13 mat : m\0324\ str\0324\ 1 n\0324\ rows substr dSymbol findstr : (1 - ≡ R\0324\ k\0324\ 1 el R\0324\ k\0324\ 2 el R\0324\ k\0324\ 3 el 13 mat U\0324\ k\0324\ el l\0324\ n\0324\ rows + 1 - p\0324\ 13 mat : R\0324\ k\0324\ 1 el R\0324\ k\0324\ 2 el R\0324\ k\0324\ 3 el 13 mat U\0324\ k\0324\ el l\0324\ m\0324\ rows + 1 - p\0324\ 13 mat : if if 31 line if 21 line for R\0324\ num2str 0 ≡ 0 k\0324\ 1 R\0324\ rows range D\0324\ D\0324\ R\0324\ k\0324\ 1 el R\0324\ k\0324\ 2 el R\0324\ k\0324\ 3 el (R\0324\ k\0324\ 2 el^DIFF : for D\0324\ 21 line if 31 line :

d(x) returns the total differential of x x\032E\ d U\0324\ x\032E\ Unknowns : s\0324\ 0 : 12 mat k\0324\ 1 : k\0324\ U\0324\ length ≤ k\0324\ k\0324\ 1 + : s\0324\ s\0324\ x\032E\ U\0324\ k\0324\ el ∂ + : for s\0324\ 31 line :

F# q# L.eq dt# dSymbol dTime concat str2num : q#' q# num2str' concat str2num : q#'' q# num2str'' concat str2num : \2202\q# dSymbol q# num2str concat str2num : \2202\q#' dSymbol q# num2str' concat str2num : 1 dt# / F# q#' ∂ \2202\q#' / d * F# q# ∂ \2202\q# / - q#'' isol num2str str2num 61 line :

gc# V# K# ∂Vars dt# dSymbol dTime concat : tt# (dTime) concat : 12 mat var# mat(dTime, concat : m# 1 : 12 mat str# line(: n# 1 : 12 mat k# K# n# n# 2 + : str# str# Clear(k# num2str), dSymbol k# num2str :0, concat : 21 line for k# V# n# n# 4 + : m# m# 2 + : 12 mat v# k# num2str : v'# v#' concat : v''# v'#' concat : 13 mat str# str# Clear(v#, v'#, v''#), concat : str# str# Clear(v# tt#, v'# tt#, v''# tt#), concat : str# str# dSymbol v# : v'# * dt#, concat : str# str# dSymbol v'# : v''# * dt#, concat : var# var# v#, v'#, concat : 71 line for var# gc# num2str : var# 1, m# num2str), concat : str# var# n# num2str,1) concat 71 line :

V# K# ∂Vars dt# dSymbol dTime concat : str# line(: n# 2 : 13 mat tt# (dTime) concat : k# K# n# n# 2 + : str# str# Clear(k# num2str), dSymbol k# num2str :0, concat : 21 line for k# V# n# n# 4 + : v# k# num2str : v'# v#' concat : v''# v'#' concat : 13 mat str# str# Clear(v#, v'#, v''#), concat : str# str# Clear(v# tt#, v'# tt#, v''# tt#), concat : str# str# dSymbol v# : v'# * dt#, concat : str# str# dSymbol v'# : v''# * dt#, concat : 61 line for str# 1, n# num2str,1) concat 51 line :

E# V# τ# L.eq str# line(: n# 1 : t# (τ# num2str) concat : 13 mat k# V# n# n# 3 + : v# k# num2str : v'# v#' concat : v''# v'#' concat : 13 mat Clear(v# t#, v'# t#, v''# t#) concat str2num str# str# v# : v# t#, concat : str# str# v'# : v'# t#, concat : str# str# v''# : v''# t#, concat : 61 line for str# 1, n# num2str,1) concat str2num E# 41 line :

θ t L'' θ δt + t L' θ t L' - δt / 11 line :

δt 109 -^:

T# A# L.at ans# 0 : r# 1 T# length range : c# 1 A# length range : 13 mat r# c# fun# dTime : T# r# el num2str concat str2num A# c# el 21 line : ans# r# c# el r# c# fun# eval : Cols 31 line :

A# c# el t# T# r# el ≡ at

at A# c# el t T# r# el ≡ feval

T# A# L.ls t# dTime str2num : Can't use dTime ans# 0 : r# 1 T# length range : c# 1 A# length range : 15 mat ans# 1 c# el T# A# c# el t lspline : 21 line :

P\032E\ Q\032E\ U\032E\ L.At n\0324\ 1 P\032E\ length range : r\0324\ 1 U\032E\ rows range : c\0324\ 1 U\032E\ cols range : v\0324\ 0 : V\0324\ 0 : C\0324\ 0 : 16 mat V\0324\ r\0324\ n\0324\ el P\032E\ n\0324\ el v\0324\ r\0324\ c\0324\ el Q\032E\ c\0324\ el U\032E\ r\0324\ c\0324\ el ≡ (: at eval : C\0324\ 1 n\0324\ el V\0324\ r\0324\ n\0324\ el : 31 line :

U\032E\ frames\032E\ digits\032E\ L.Frames Fn\0324\ 1 U\032E\ rows frames\032E\ pR : T\0324\ U\032E\ 1 col : Fn\0324\ t = T\0324\ digits\032E\ var2str concat vectorize 12 mat eval 31 line :

Fn\0324\ Fn\0324\ Fn\0324\ reverse stack eval :

dTime t :

dSymbol ∂ :

<content> equrep </content>

\032E\a TrueQ \032E\a num2str 1 ≡ 11 line :

\032E\a FalseQ \032E\a num2str 0 ≡ 11 line :

\032E\a b\032E\ EqualQ \032E\a num2str b\032E\ num2str ≡ 11 line :

E\032E\ \032E\a b\032E\ EquRep equrep(E\032E\ num2str, \032E\a num2str, b\032E\ num2str) concat str2num 11 line :

E\032E\ \032E\a b\032E\ q\032E\ EquRep equrep(E\032E\ num2str, \032E\a num2str,if(TrueQ(q\032E\ num2str), b\032E\ num2str, a# num2str)) concat str2num 11 line :

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

a#

b#

q#

E\032E\ Rules\032E\ EquRep ans\0324\ E\032E\ : n\0324\ 1 Rules\032E\ rows range ans\0324\ ans\0324\ Rules\032E\ n\0324\ 1 el Rules\032E\ n\0324\ 2 el EquRep : for ans\0324\ 31 line :

E\032E\ Rules\032E\ EquRepN ans\0324\ E\032E\ : N\0324\ 1 EquRepN range ans\0324\ ans\0324\ ans\0324\ Rules\032E\ EquRep : EqualQ break continue if 11 line for ans\0324\ 31 line :

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

E\032E\ x\032E\ depends E\032E\ x\032E\ x\032E\_DEPENDS.X# EquRep num2str x\032E\_DEPENDS.X# num2str findstr 1 - ≠ 11 line :

\032E\a b\032E\ x\032E\ depends \032E\a x\032E\ depends b\032E\ x\032E\ depends & 11 line :

E\032E\ x\032E\ DIFF diff(E\032E\ x\032E\_DIFF# EquRep num2str,_DIFF#) concat str2num_DIFF# x\032E\ EquRep 11 line :

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

EquRepN 20 :

EQ EQ(-a,b) EQ(a,0-b) 1 minus 14 mat :

EQ EQ EQ(a^b,z) EQ.B(a,b,z,a^b,ln(z)/ln(a),'sign*z^(1/b)) 1 power 14 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) 1 quotient 14 mat stack :

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

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

EQ EQ EQ(a,z) EQ.L(a,z) 1 lineal 14 mat stack :

EQ EQ EQ(x*(b+a*x),z) EQ.P2(a,b,x,z) 1 Poly2 14 mat stack :

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

EQA rows 15

_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 1 EQA_k# 1 el 14 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 1 EQA_k# 2 el 14 mat stack : for

E\032E\ x\032E\ isolate a\0324\_EQ#.10 10 a\0324\^11 line : a\0324\ b\0324\ f\0324\ a\0324\ EQ.A a\0324\ x\032E\ depends a\0324\ a\0324\ EQ f\0324\ b\0324\ EQ if 11 line : a\0324\ b\0324\ z\0324\ f\0324\ a\0324\ j\0324\ EQ.B a\0324\ x\032E\ depends b\0324\ x\032E\ depends \0324\k f\0324\ x\032E\ 0 EquRep : f\0324\ \0324\k - z\0324\ \0324\k - EQ 21 line a\0324\ j\0324\ EQ if b\0324\ x\032E\ depends b\0324\ a\0324\ EQ z\0324\ f\0324\ EQ if if 11 line : a\0324\ z\0324\ EQ.L a\0324\ x\032E\ 2 DIFF FalseQ x\032E\ z\0324\ a\0324\ x\032E\ 0 EquRep - a\0324\ x\032E\ DIFF / EQ a\0324\ z\0324\ EQ if 11 line : a\0324\ b\0324\ \0324\X z\0324\ EQ.P2 \0324\X x\032E\ EqualQ a\0324\ b\0324\ x\032E\ depends ¬ (& a\0324\ \0324\X b\0324\ 2 a\0324\ * / + (2^* z\0324\ b\0324\ 2^4 a\0324\ * / + EQ \0324\X b\0324\ a\0324\ \0324\X * + (* z\0324\ EQ if 11 line : ans\0324\ E\032E\ equal2EQ : A\0324\ ans\0324\ 1 el numden : B\0324\ ans\0324\ 2 el numden : 13 mat ans\0324\ A\0324\ 1 el B\0324\ 2 el * A\0324\ 2 el B\0324\ 1 el * - 0 EQ EQ EquRepN : a\0324\ b\0324\ EQ a\0324\ b\0324\ ≡ (: ans\0324\ 91 line :

E\032E\ x\032E\ isol ans\0324\ E\032E\ x\032E\ isolate equal2EQ : ans\0324\ 1 el x\032E\ EqualQ ans\0324\ 2 el x\032E\ depends ¬ & ans\0324\ 2 el ans\0324\ 1 el ans\0324\ 2 el - if 21 line :

sign 1 :

<content> Numerical Inverse Laplace Transform </content>

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 :

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

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 :

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 :

1 F t t.i t.e N ILT FT_s# F t.i t.e N ILT : U FT 11 el UnitsOf FT 12 el UnitsOf 12 mat : U 2 el FT 1 col U 1 el / FT 2 col U 2 el / t U 1 el / cinterp * t.i t ≤ (* t t.e ≤ (* 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 :

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

<content> laplace </content>

Euler Mascheroni constant γ 0.577215664901532860606 :

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

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

exp

LT LT e exp(1) e 13 mat stack :

LT LT exp(1)^a exp(a) e 13 mat stack :

LT LT a/exp(b) a*exp(-b) e 13 mat stack :

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

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

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

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

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

LT LT LT(ln(a)) LT.φω(ln(a),a,(ln(lω)-ln('svar)-'γ)/'svar 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) if((¬depends(n,'tvar)&IsT#(a),(-1)^n*DIFF(LT(b),'svar,n),LT(a^n*b) nshift 13 mat stack :

LT LT LT(a^n) if(IsDefined(n)&IsT#(a),n!*'svar^(-n-1),LT(a^n)) pot 13 mat stack :

LT LT LT(a*b) if(¬depends(a,t),a*LT(b),if(IsT#(a),-DIFF(LT(b),'svar),LT(b*a))) product 13 mat stack :

Conmute!

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 :

E# lt lω lφ Clear a# IsT# a# num2str tvar num2str ≡ 11 line : LT.O tvar num2str svar 2^num2str 12 mat : LT.O LT.O x(t) X(s) 12 mat stack : xLT.O LT.O tvar u num2str svar 1 -^num2str 12 mat stack : xLT.O LT.O tvar δ num2str 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 : a# LT.0 a# tvar depends row# LT.O a# num2str 1 findrows : (0 ≡ a# LT row# 2 el str2num if a# svar 1 -^* if 11 line : A# B# a# b# c# LT.α lω b# tvar DIFF : lα a# tvar DIFF : IsLT? 13 mat lα tvar depends lω tvar depends | (c# tvar depends | c# A# * B# * LT c# B# lt * svar svar lα - EquRep if 21 line : f# a# tr# LT.φω lω a# tvar DIFF : lφ a# lω tvar * - : 12 mat a# tvar depends lω tvar depends f# LT tr# if a# svar 1 -^* if 21 line : E# num2str ≡ findstr 1 - ≡ E# LT LT EquRepN A# B# 12 mat E# equal2EQ : A# LT LT EquRepN B# LT LT EquRepN ≡ 21 line if 131 line :

a# b# LT.prod a# t depends ¬ a# b# LT * a# IsT# b# LT s DIFF - a# b# * LT if if 11 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 :

a# n# LT.pot n# IsDefined a# IsT# & n#! (s n# 1 + (-^* a# n#^LT if 11 line :

ans# E# sqrt(a) a^(1/2) nthroot(a,b) a^(1/b) 22 mat EcuRep :

E# t# s# lt tvar svar Clear tvar t# : svar s# : 12 mat E# lt 31 line :

tvar t :

svar s :

xvar t x :

Xvar s X :

M Cols C 0 : k 1 M cols range : 12 mat C 1 k el M k col : 21 line :

A D C B TriDiagSys n B length : v n 1 matrix : y v : 13 mat w A 1 el 1 -^: y 1 el B 1 el w * : 12 mat k 2 n range v k 1 - el C k 1 - el w * eval : w A k el D k el v k 1 - el * - (1 -^eval : y k el B k el D k el y k 1 - el * - (w * eval : 31 line for k n 1 - (1 range y k el y k el v k el y k 1 + el * - eval : for y 51 line :

x y pStem XY 03 matrix : S 02 matrix : 12 mat k 1 x length range XY XY x k el y k el o 13 mat stack : S S x k el 0 x k el y k el x k el 032 mat stack : 21 line for XY S 21 sys eval 31 line :

lbvp.2(f(x),[a,b],M,n) solves the linear BVP y'' + p(x)y' + q(x)y = φ(x) in a < x < b with    f(x) = [p(x), q(x), φ(x)],  M = [α1 β1 c1; α2 β2 c2]  and     α1*y(a) + β1*y'(b) = c1   and  α2*y(a) + β2*y'(b) = c2using central differences. 1 f x M n lbvp.2 h x 2 el x 1 el - n / : X x 1 el 0 n range h * + : ν n 1 + : P 0 : Q 0 : Φ 0 : α1 β1 c1 α2 β2 c2 23 mat M : 17 mat r 1 ν range P r el Q r el Φ r el 13 mat X r el f eval : for A D C B 14 mat 2 - h 2^Q * + (2 h P * - 2 / (2 h P * + 2 / (h 2^Φ * (14 mat eval : A 1 el D 1 el C 1 el B 1 el 14 mat β1 0 ≡ 100 c1 α1 / 14 mat A 1 el 2 h * α1 β1 / * D 1 el * + (02 B 1 el 2 h * c1 β1 / * D 1 el * + (14 mat if eval : A ν el D ν el C ν el B ν el 14 mat β2 0 ≡ 100 c2 α2 / 14 mat A ν el 2 h * α2 β2 / * C ν el * - (20 B ν el 2 h * c2 β2 / * C ν el * - (14 mat if eval : Y A D C B TriDiagSys : Y' Y : j 2 n range : 13 mat Y' j el Y j 1 + el Y j 1 - el - 2 h * / : Y' 1 el β1 0 ≡ Y 2 el Y 1 el - h / c1 β1 / α1 β1 / Y 1 el * - if : Y' ν el β2 0 ≡ Y ν el Y n el - h / c2 β2 / α2 β2 / Y ν el * - if : 13 mat X Y Y' augment eval 81 line :

bvp.2(f(x),[a,b],M,n) solves the linear BVP y'' = f(x,y,y') in a < x < b with    f(x) = [p(x), q(x), φ(x)] , M = [α1 β1 c1; α2 β2 c2]  and     α1*y(a) + β1*y'(b) = c1   and  α2*y(a) + β2*y'(b) = c2using central differences. 3 f x Y M n ε bvp.2 h x 2 el x 1 el - n / : X x 1 el 0 n range h * + : ν n 1 + : P 0 : Q 0 : Φ 0 : α1 β1 c1 α2 β2 c2 23 mat M : 17 mat Y 0 ≡ ψ x 1 el α1 * β1 + α1 x 2 el α2 * β2 + α2 22 mat (1 -^c1 c1 21 mat * : Y ψ 2 el ψ 1 el X * + : 12 mat 0 if Δy Y max Y min - ν / (2^h 2^+ sqrt eval : Δy' Δy h 1 -^* : 12 mat x y y' fy x y Δy + y' f x y Δy - y' f - 2 Δy * / : x y y' fy' x y y' Δy' + f x y y' Δy' - f - 2 Δy' * / : A ν 1 matrix eval : D A : C A : B A : Y' A : j 2 n range : 16 mat Iter 1 NewtonMaxIters range Y' j el Y j 1 + el Y j 1 - el - 2 h * / : Y' 1 el β1 0 ≡ Y 2 el Y 1 el - h / c1 β1 / α1 β1 / Y 1 el * - if : Y' ν el β2 0 ≡ Y ν el Y n el - h / c2 β2 / α2 β2 / Y ν el * - if : 13 mat r 1 ν range P r el Q r el Φ r el 13 mat X r el Y r el Y' r el fy X r el Y r el Y' r el fy' X r el Y r el Y' r el f 13 mat eval : for A D C 13 mat 2 - h 2^P * - (2 h Q * - 2 / (2 h Q * + 2 / (13 mat eval : B j el 2 Y j el * - Y j 1 + el + Y j 1 - el + h 2^Φ j el * - eval : A 1 el D 1 el C 1 el B 1 el 14 mat β1 0 ≡ 100 Y 1 el c1 α1 / - 14 mat A 1 el 2 h * α1 β1 / * B 1 el * + (022 Y 2 el * 21 - h α1 β1 / * + (* Y 1 el * + 2 h * c1 β1 / * - h 2^Φ 1 el * - (14 mat if eval : A ν el D ν el C ν el B ν el 14 mat β2 0 ≡ 100 Y ν el c2 α2 / - 14 mat A ν el 2 h * α2 β2 / * C ν el * - (202 Y n el * 21 h α2 β2 / * + (* Y ν el * - 2 h * c2 β2 / * + h 2^Φ ν el * - (14 mat if eval : Yo A D C B - TriDiagSys : Y Y Yo + eval : 12 mat Yo norme ε Y norme * < break continue if 81 line for Iter NewtonMaxIters < X Y Y' augment eval'NewtonMaxIters  rached. error if 81 line :

NewtonMaxIters 50 :

f.feval# x.feval# feval f.feval# num2str (x.feval# num2str) concat str2num 11 line :

f.feval# x.feval# y.feval# feval f.feval# num2str (x.feval# num2str, y.feval# num2str) concat str2num 11 line :

a\032E\ b\032E\ n\032E\ pR a\032E\ b\032E\ a\032E\ - n\032E\ / 0 n\032E\ range * + eval 11 line :

EQ# Y# X# A# V# Ode2Aux Q# EQ# Sys2EQ : B# 0 : k# 0 : ys# Y# num2str : 14 mat y'# y''# 12 mat ys#''' 12 mat concat str2num vectorize : a# A# at# ys# ('(12 mat a# num2str) concat vectorize B# 1 k# k# 1 + : el at# str2num : Q# Q# at# V# k# el EquRep : 21 line for for k# 0 : M# 0 : 12 mat q# Q# q# y''# depends ODE# q# : M# 1 k# k# 1 + : el q# : if for ODE# M# B# y'# y''# 12 mat 14 mat 61 line :

eq# y# x# a# N# sb_lbvp.2 V# X# Clear ODE# Mq# B# y'# y''# 12 mat 14 mat eq# y# x# a# V# Ode2Aux : p# v1# v2# Z# p# v1# 0 EquRep v2# 0 EquRep - 11 line : M# Mq# 1 el V# 1 el DIFF Mq# 1 el V# 2 el DIFF Mq# 1 el V# 1 el V# 2 el Z# Mq# 2 el V# 3 el DIFF Mq# 2 el V# 4 el DIFF Mq# 2 el V# 3 el V# 4 el Z# 23 mat : Lq# ODE# y''# isol : L# Lq# y'# DIFF Lq# y# DIFF Lq# y# y'# Z# 13 mat x# X# EquRep : X# f# L# - : X# f# a# M# eval N# lbvp.2 81 line :

eq# y# x# a# N# YO# E# sb_bvp.2 V# X# Y# Y'# Clear ODE# Mq# B# y'# y''# 12 mat 14 mat eq# y# x# a# V# Ode2Aux : p# v1# v2# Z# p# v1# 0 EquRep v2# 0 EquRep - 11 line : M# Mq# 1 el V# 1 el DIFF Mq# 1 el V# 2 el DIFF Mq# 1 el V# 1 el V# 2 el Z# Mq# 2 el V# 3 el DIFF Mq# 2 el V# 4 el DIFF Mq# 2 el V# 3 el V# 4 el Z# 23 mat : L# ODE# y''# isol x# X# EquRep y# Y# EquRep y'# Y'# EquRep : X# Y# Y'# f# L# : X# Y# Y'# f# a# YO# M# eval N# E# bvp.2 71 line :

<content> Short hands for plots </content>

Short hands for plots

x xb X Y Y' x y x y' G Ys Y's Clear 1

Err X Y Y' 13 mat sol Cols : Y ψ X feval vectorize - norme 21 line :

x ψ' x ψ x diff 11 line :

str sol Plot X Y Y' 13 mat sol Cols : Xs X 1 el X X rows el 300 pR : str XY ≡ u v us vs 14 mat X Y Xs ψ Xs feval vectorize 14 mat eval : u v us vs 14 mat Y Y' ψ Xs feval vectorize ψ' Xs feval vectorize 14 mat eval : if us vs augment u 1 el u u rows el 12 mat v 1 el v v rows el 12 mat pStem u v . 8 augment 31 sys 41 line :

EQ EQ#(-a,b) EQ#(a,0-b) 1 minus 14 mat :

EQ EQ EQ#(a^b,z) EQ#.B(a,b,z,a^b,ln(z)/ln(a),'sign*z^(1/b)) 1 power 14 mat stack :

EQ EQ EQ#(a/b,z) EQ#.B(a,b,z,a/b,z/a,z*b) 1 quotient 14 mat stack :

EQ EQ EQ#(a*b,z) EQ#.B(a,b,z,a*b,z/a,z/b) 1 times 14 mat stack :

EQ EQ EQ#(a+b,z) EQ#.B(a,b,z,a+b,z-a,z-b) 1 plus 14 mat stack :

EQ EQ EQ#(a,z) EQ#.L(a,z) 1 lineal 14 mat stack :

EQ EQ EQ#(x*(a*x+b),z) EQ#.P2(a,b,x,z) 1 Poly2 14 mat stack :

EQ EQ EQ#(x*(b+a*x),z) EQ#.P2(a,b,x,z) 1 Poly2 14 mat stack :

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

EQA rows 15

_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 1 EQA_k# 1 el 14 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 1 EQA_k# 2 el 14 mat stack : for

E\032E\ x\032E\ RU\032E\ isolate a\0324\_EQ#.10 10 a\0324\^11 line : a\0324\ b\0324\ f\0324\ h\0324\ EQ#.A a\0324\ #ISOL# depends a\0324\ h\0324\ EQ# f\0324\ b\0324\ EQ# if 11 line : a\0324\ b\0324\ z\0324\ f\0324\ h\0324\ j\0324\ EQ#.B a\0324\ #ISOL# depends b\0324\ #ISOL# depends k\0324\ f\0324\ #ISOL# 0 EquRep : f\0324\ k\0324\ - z\0324\ k\0324\ - EQ# 21 line a\0324\ j\0324\ EQ# if b\0324\ #ISOL# depends b\0324\ h\0324\ EQ# z\0324\ f\0324\ EQ# if if 11 line : a\0324\ z\0324\ EQ#.L a\0324\ #ISOL# 2 DIFF 0 EqualQ a\0324\ #ISOL# depends & #ISOL# z\0324\ a\0324\ #ISOL# 0 EquRep - a\0324\ #ISOL# DIFF / EQ# a\0324\ z\0324\ EQ# if 11 line : a\0324\ b\0324\ X\0324\ z\0324\ EQ#.P2 X\0324\ #ISOL# - EqualQ a\0324\ b\0324\ #ISOL# depends ¬ (& a\0324\ - b\0324\ X\0324\ - z\0324\ EQ#.P2 X\0324\ #ISOL# EqualQ a\0324\ b\0324\ #ISOL# depends ¬ (& X\0324\ b\0324\ - sign 4 z\0324\ * a\0324\ * b\0324\ 2^+ sqrt * + 2 a\0324\ * / EQ# X\0324\ b\0324\ a\0324\ X\0324\ * + (* z\0324\ EQ# if 11 line if 11 line : ans\0324\ E\032E\ x\032E\ #ISOL# 12 mat EquRep equal2EQ : A\0324\ ans\0324\ 1 el numden : B\0324\ ans\0324\ 2 el numden : 13 mat ans\0324\ A\0324\ 1 el B\0324\ 2 el * A\0324\ 2 el B\0324\ 1 el * - 0 EQ# RU\032E\ EquRepN : a\0324\ b\0324\ EQ# a\0324\ b\0324\ ≡ (: #ISOL# x\032E\ : 12 mat ans\0324\ 91 line :

E\032E\ x\032E\ isolate E\032E\ x\032E\ EQ isolate 11 line :

E\032E\ x\032E\ isol ans\0324\ E\032E\ x\032E\ isolate equal2EQ : ans\0324\ 1 el x\032E\ EqualQ ans\0324\ 2 el x\032E\ depends ¬ & ans\0324\ 2 el ans\0324\ 1 el ans\0324\ 2 el - if 21 line :

sign 1 :

Converts an expression with binary operators to its functional representation E\032E\ N2E k 1 N2E cols range ans Clear(N2E 2 k el (a,b)) concat : ans str2num 12 mat for T\032E\ b\032E\ #N2E# equrep(T\032E\,a¤b, b\032E\) concat str2num : X\032E\ a\032E\ b\032E\ #N2E# t\0324\ X\032E\ num2str {(strrep}) strrep : t\0324\ t\0324\ (#N2E#(strrep a\032E\ ¤ strrep : N\0324\ 1 EquRepN range t\0324\ num2str t\0324\ t\0324\ num2str b\032E\ #N2E# : num2str ≡ break continue if 11 line for t\0324\ num2str #N2E#((strrep str2num 41 line : ans\0324\ E\032E\ a b ¤ a b #Xor# equrep : ans\0324\ ans\0324\ a b & a b #And# equrep : ans\0324\ ans\0324\ a b | a b #Or# equrep : ans\0324\ ans\0324\ a b > b a < equrep : ans\0324\ ans\0324\ a b ≥ b a ≤ equrep : ans\0324\ ans\0324\ a b ≡ a b #Equals# equrep : ans\0324\ ans\0324\ a b ≠ a b #NotEquals# equrep : ans\0324\ ans\0324\ a b < a b #LessThan# equrep : ans\0324\ ans\0324\ a b ≤ a b #LessOrEquals# equrep : ans\0324\ ans\0324\ num2str - +#neg#* strrep str2num : k\0324\ 13 range ans\0324\ ans\0324\ num2str N2E 1 k\0324\ el # N2E 2 k\0324\ el #(a,b) concat #N2E# : for ans\0324\ ans\0324\ / #Times#(a,#Inv#(b) #N2E# : k\0324\ 1 N2E cols range # N2E 2 k\0324\ el #(a,b): N2E 2 k\0324\ el (a,b) concat str2num for b #Inv# b #neg# Pow : #neg# 1 - : b\032E\ #N2E# b\032E\ : 13 mat ans\0324\ num2str #N2E# strrep num2str str2num 181 line :

This avoid the bug:

a b & (c | (d ≡ num2str str2num a b & (c d ≡ (|

This is, can't use num2str with booleans

a b & (c | (d ≡ num2str a&b|c≡d

Matrix of binary operators N2E + *^& | ≡ ≠ < ≤ ¤ Plus Times Pow And Or Equals NotEquals LessThan LessOrEquals Xor 210 mat :

E\032E\ x\032E\ isolate a\0324\_EQ#.10 10 a\0324\ Pow : a\0324\ b\0324\ f\0324\ h\0324\ EQ#.A a\0324\ #ISOL# depends a\0324\ h\0324\ EQ# f\0324\ b\0324\ EQ# if 11 line : a\0324\ b\0324\ z\0324\ f\0324\ h\0324\ j\0324\ EQ#.B a\0324\ #ISOL# depends b\0324\ #ISOL# depends k\0324\ f\0324\ #ISOL# 0 EquRep : f\0324\ #neg# k\0324\ Times Plus z\0324\ #neg# k\0324\ Times Plus EQ# 21 line a\0324\ j\0324\ EQ# if b\0324\ #ISOL# depends b\0324\ h\0324\ EQ# z\0324\ f\0324\ EQ# if if 11 line : ans\0324\ E\032E\ x\032E\ #ISOL# 12 mat EquRep equal2EQ : A\0324\ ans\0324\ 1 el numden : B\0324\ ans\0324\ 2 el numden : ans\0324\ A\0324\ 1 el B\0324\ 2 el * A\0324\ 2 el B\0324\ 1 el * - : 13 mat ans\0324\ #ISOL# 2 DIFF 0 EqualQ ans\0324\ #ISOL# depends & ans\0324\ #ISOL# ans\0324\ #ISOL# 0 EquRep ans\0324\ #ISOL# DIFF / - EQ# : ans\0324\ ans\0324\ N2E 0 EQ# isol EquRepN : if a\0324\ b\0324\ EQ# a\0324\ b\0324\ ≡ (: #ISOL# x\032E\ : #neg# 1 - : 13 mat k\0324\ 1 N2E cols range N2E 2 k\0324\ el (a,b): a N2E 1 k\0324\ el b concat str2num for ans\0324\ 91 line :

E\032E\ x\032E\ isol ans\0324\ E\032E\ x\032E\ isolate equal2EQ : ans\0324\ 1 el x\032E\ EqualQ ans\0324\ 2 el x\032E\ depends ¬ & ans\0324\ 2 el ans\0324\ 1 el ans\0324\ 2 el - if 21 line :

isol 03 matrix eval :

isol isol EQ#(Pow(a,b),z) EQ#.B(a,b,z,a^b,ln(z)/ln(a),exp(ln(z)/b)) 113 mat stack :

Needs sign

isol isol EQ#(Pow(a,b),z) EQ#.B(a,b,z,Pow(a,b),Times(ln(z),Pow(ln(a),#neg#)),Pow(z,Times(1,Pow(b,#neg#)))) 113 mat stack :

isol isol EQ#(Times(a,b),z) EQ#.B(a,b,z,Times(a,b),Times(z,Pow(a,#neg#)),Times(z,Pow(b,#neg#))) 113 mat stack :

isol isol EQ#(Plus(a,b),z) EQ#.B(a,b,z,Plus(a,b),Plus(z,Times(#neg#,a)),Plus(z,Times(#neg#,b))) 113 mat stack :

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

_k# 1 EQA rows range isol isol EQ#(EQA_k# 1 el (a),b) concat EQ#.A(a,b, EQA_k# 1 el (a), EQA_k# 2 el (b) concat 113 mat stack : for

_k# 1 EQA rows range isol isol EQ#(EQA_k# 2 el (a),b) concat EQ#.A(a,b, EQA_k# 2 el (a), EQA_k# 1 el (b) concat 113 mat stack : for

ison isol eval :

_k# Clear 1

s# s2n lp# rp# 12 mat s# {findstr s#} findstr 12 mat : out# r# 12 mat 112 mat : n# 1 lp# rows range out# out# s# r# lp# n# el r# - substr concat : var# s# lp# n# el 1 + rp# n# el lp# n# el - 1 - substr str2num : out# out# var# num2str concat : r# rp# n# el 1 + : 41 line for out# s# r# s# strlen r# - 1 + substr concat str2num 41 line :

eq# x# msol2 maple(evalf(solve({eq#},{x#}))) s2n 11 line :

X_Concat# Concat str_Concat# : x_Concat# X_Concat# str_Concat# str_Concat# x_Concat# num2str concat : for str_Concat# 31 line :

eq# u# t# U# ode2 f# υ# u2f# concat(" f# ",num2str(υ#),"_",num2str(t#) 17 mat Concat 11 line : uk# u# tmp uk# (t#):str2num(fun_uk# u2f#) 16 mat Concat str2num : tmp uk#'(t#):str2num(fun_' uk# u2f#) 16 mat Concat str2num : 21 line for k# 0 : D# 01 matrix eval : ic# D# : ic'# D# : 14 mat vk# eq# Unknowns vk# num2str fun_findstr 1 - ≡ D# k# k# 1 + : el at vk# eq# feval : vk# num2str fun_' findstr 1 - ≡ ic# ic# at vk# eq# feval stack eval : ic'# ic'# at vk# eq# feval stack eval : if if for D# at D# U# 1 el U# 3 el 21 sys feval : U# 2 el D# 2 el U# 4 el D# 4 el 41 mat ic# 1 el ic'# 1 el ic# 2 el ic'# 2 el 41 mat 12 mat 61 line :

<content> Differentials samples </content>

Differentials

dSymbol d :

Enter an expression, and say that have some constants

z x y 2^* y a x * / + :

da 0 :

dz z d : y 1 x 2^y * a * - (* dx * - 12 x 2^* y * a * + (dy * x * + a x 2^* /

Total differential

Partial derivatives

dz dx ∂ dx / y 1 x 2^y * a * - (* a x 2^* / -

dz dy ∂ dy / 12 x 2^* y * a * + x a * /

z x ∂ dx / y 1 x 2^y * a * - (* a x 2^* / -

z y ∂ dy / 12 x 2^* y * a * + a x * /

Some examples

dSymbol ∂. :

∂.a 0 :

r θ cos * a cos x * r θ sin * a x * 22 mat d θ cos ∂.r * r θ sin * ∂.θ * - a cos ∂.x * θ sin ∂.r * r θ cos * ∂.θ * + a ∂.x * 22 mat

dSymbol ∂ :

∂a 0 :

a e b x *^* d e b x *^a * x ∂b * b ∂x * + (*

a e x sin^* d e x sin^x cos * a * ∂x *

x 2^x y * ln * d x 12 x y * ln * + (∂x * y * x ∂y * + (* y /

∂a Clear 1

d1 a b * d : b ∂a * a ∂b * +

d2 d1 d : ∂a b ∂b + ∂b + (* a ∂b * +

d3 d2 d : ∂a b 3 ∂b * + ∂b + (* a 2 ∂a * + (∂b * +

∂a 0 :

d1 a b * d : a ∂b *

d2 d1 d : a ∂b *

d3 d2 d : a ∂b *

Chain rule

c a x g f * d : a x g f x diff * ∂x *

x f x ln :

x g x cos :

c a x tan * ∂x * -

a x cos ln * d x sin a * ∂x * x cos / -

Total differential

t x x' f d t x x' f t diff ∂t * t x x' f x diff ∂x * + t x x' f x' diff ∂x' * +

Alvaro

∞ 1010^:

CycleColors 5 :

X X' X'' TX N pL s1 G({X},{X'}) s2n num2str : s2 eval(mat({X'},{X''},2,1)) s2n num2str : s1 : s2 concat str2num G TX N pL 41 line :

Ds TX N pL B TX 2 col min TX 2 col max TX 3 col min TX 3 col max 22 mat eval : Ds num2str B N pVField TX 2 col TX 3 col augment TX 12 el TX 13 el o 10 red augment 31 sys 21 line :

¦AuxIPE 0931061551212515610390116 mat :

¦AuxIPΔ 0110 0.5 1100 0.5 25 mat :

¦AuxIPT 0112211331122111042434344323424000000200003000000000030000400000416 mat :

¦AuxIPB 0110011001100110001111000011110000001111111100000101010110101010416 mat :

2 F B N pIPlot x1 x2 y1 y2 22 mat B : nx N 1 el : ny N 2 el : M 0 : R 0 : v 0 : P 51 matrix : ϕ P : 18 mat Δx x2 x1 - nx / eval : Δy y2 y1 - ny / eval : κ 15 range : 13 mat c 0 nx 1 - (range r 0 ny 1 - (range P κ el ¦AuxIPΔ 1 κ el c + (Δx * x1 + ¦AuxIPΔ 2 κ el r + (Δy * y1 + 21 mat eval : ϕ κ el P κ el 1 el P κ el 2 el F eval 0 try : n 1 : k 14 range ϕ k el 0 < n n 2 k 1 -^+ : continue if for 031 line ϕ 5 el 0 < n 11 ≡ n 6 : n 6 ≡ n 11 : 0 if if 0 if 12 mat ¦AuxIPE n el 0 ≠ k 14 range ¦AuxIPB k n el 0 ≠ k1 k 1 + 4 k 4 ≡ (* - eval : v k el P k el ϕ k el ϕ k1 el ϕ k el - / P k1 el P k el - (* - eval : 21 line continue if for k 1 : k 4 < (¦AuxIPT k n el 0 ≠ (& k k 2 + : M R R 1 + : el v ¦AuxIPT k n el el v ¦AuxIPT k 1 + n el el augment : η 1 CycleColors range M R R 1 + : el ∞ ∞ 12 mat : for 21 line for 21 line continue if 41 line for for M transpose vectorize mat2sys.1 eval 41 line :

P ϕ k k1 Interp P k el ϕ k el ϕ k1 el ϕ k el - / P k1 el P k el - (* - eval 11 line :

Ber7 Cycle colors

2 F B N Λ pIPlot x# y# Clear M 0 : R 0 : 13 mat λ Λ x y Φ x y F λ eval - 11 line : M R R 1 + : el x# y# Φ B N pIPlot : M R R 1 + : el 11 sys : 31 line for M mat2sys.1 eval 31 line :

P V L sf pVec t V norme eval : a 0.5 0.3 0.15 0.65 0.35 14 mat * sf t feval * eval : x L 1 el a 1 el - 000 a 1 el 51 mat * eval : y L 2 el a 3 el a 4 el - a 3 el a 4 el + a 3 el a 4 el + (- a 3 el a 4 el + a 3 el a 4 el - 51 mat * eval : 12 mat θ t 0 ≡ 0 V 1 el sign V 2 el t / acos * - if eval : c θ cos eval : s θ sin eval : 13 mat k 15 range x k el y k el 21 mat P c x k el * s y k el * - s x k el * c y k el * + 21 mat + eval : for x y augment eval 61 line :

Fs B N L sf pVField m kg s K A cd 16 mat 11111116 mat : x1 x2 y1 y2 22 mat B eval : nx N 1 el : ny N 2 el : n nx ny * : R 0 : 15 mat M 5 n * 5 + 1 matrix : V n 1 matrix : P V : 13 mat Δx x2 x1 - nx / eval : Δy y2 y1 - ny / eval : 12 mat c 0 nx 1 - (range r 0 ny 1 - (range x y 21 mat c 0.5 + (Δx * x1 + r 0.5 + (Δy * y1 + 21 mat : ν r c ny * + 1 + : P ν el x y 21 mat : V ν el Fs x y feval eval P ν el try : 13 mat 21 line for for V V Δy V norme vectorize max / * eval : k 1 n range M R R 1 + : el P k el V k el L sf pVec : η 1 CycleColors range M R R 1 + : el ∞ : for 21 line for M mat2sys.1 eval 81 line :

Fs B N pVField t sf# t : 11 mat Fs B N 1112 mat sf# pVField 21 line :

ic to tf N 2 D pRKFam S IC t u Clear k 1 ic length eval range : 12 mat S k el ic k el to tf N t u D rkfixed : IC k el ic k el transpose o augment : mat2sys.1 S k el S k el 2 col S k el 3 col augment : mat2sys.1 21 sys 31 line :

Plots a family of rk by various ic

FamilyPlot(f(α,x),α,x) plots f for the ranges α and x. 2 f a b c d M N pFamily x# y# Clear x a b M LinSpace : y c d N LinSpace : 13 mat x# y# f x y pFamily 21 line :

FamilyPlot(f(α,x),α,x) plots f for the ranges α and x. 2 f α x pFamily M α length : N x length : P 0 : F 0 : 14 mat c 1 M range r 1 N range F r c el α c el x r el f Re eval 0 try : for for c 1 M range P c el x F c col augment eval : for P mat2sys.1 eval 41 line :

1 f 1 g T N pParametric t T 1 el T 2 el N LinSpace : t f vectorize t g vectorize 12 mat 21 line :

1 f T N pParametric t T 1 el T 2 el N LinSpace : X N 1 matrix : Y X : F t f vectorize : 13 mat k 1 N range f F k el : X k el f 1 el : Y k el f 2 el : 13 mat for X Y augment 41 line :

x Round x Re 6 round i x Im 6 round * + :

<content> Lagrangians </content>

Lagrangian Dynamics

1 \2202\t / L q' ∂ \2202\q' / d * L q ∂ \2202\q / - q'' isol

x id x :

too 0 time :

frames 99 :

NN 25512 mat * :

D Clear 1

01033 -^:

m kg s K A cd 16 mat 11111116 mat :

dTime t :

dSymbol ∂ :

sign 1 :

2 D U# a# b# N# sf# L.vf m kg s K A cd 16 mat 11111116 mat : α# β# 12 mat a# b# 12 mat 1 + : 12 mat N# 0 ≡ (sf# 0 ≡ (| U# α# col U# β# col augment U# 1 α# el U# 1 β# el o augment 21 sys eval xx# xx'# Clear u# U# 12 U# cols range el : u# a# el xx# : u# b# el xx'# : 14 mat G# U# 11 el u# D a# b# 12 mat el : xx# xx'# G G# : 12 mat t SF# sf# str2num : U# α# col U# β# col augment U# 1 α# el U# 1 β# el o augment G U# α# col U# β# col 10 % * 10 % * pBox N# 1112 mat SF# pVField 31 sys eval 41 line if 21 line :

<content> Mechanics Braintwister Problem </content>

Mechanics Braintwister Problem

https://fr.maplesoft.com/applications/view.aspx?SID=131117&view=html

Gen coords, constants & position

gc r θ 12 mat a m.1 m.2 g 14 mat ∂Vars str2num t r r' θ θ' 15 mat

x1 r 2 a * - :

x2 r θ sin * :

y2 r θ cos * - :

x1' x1 d ∂t / :

x2' x2 d ∂t / :

y2' y2 d ∂t / :

Lagrangian

L 12 / m.1 * x1' 2^* 12 / m.2 * x2' 2^y2' 2^+ (* + m.2 g * y2 * - :

Equations of motion

r''.eq L r L.eq :

θ''.eq L θ L.eq :

L m.1 r' 2^* m.2 θ sin r' * r θ cos * θ' * + (2^θ cos r' * - r θ sin * θ' * + (2^+ (* + 2 m.2 * g * r * θ cos * + 2 /

r''.eq θ' θ cos θ sin r' * r θ cos * θ' * + (* θ sin θ cos r' * - r θ sin * θ' * + (* + (* g θ cos * + (m.2 * m.1 m.2 + /

θ''.eq θ cos θ sin r' * r θ cos * θ' * + (* θ sin θ cos r' * - r θ sin * θ' * + (* + θ' r * + (r' * r g * θ sin * + r 2^/ -

Solutions

N 200 :

te 2 s * :

a m.1 m.2 g 14 mat 50 cm * 50 g * 50 g * g.e 14 mat :

ic a 0 m s / * 90 ° * 0 rpm * 14 mat :

t u D r r' θ θ' 14 mat u 1 el u 2 el u 3 el u 4 el 14 mat : r' r''.eq θ' θ''.eq stack 21 line :

U t u D 0 te 12 mat ic N Rkadapt L.ode :

Find when m1 hit the pulley

X1 11 mat x1 11 mat gc U L.At :

n X1 0 < vectorize eval 11 findrows sum : 52

Solve for a simple pendulum, keeping r constant

t u D2 r r' θ θ' 14 mat u 1 el u 2 el u 3 el u 4 el 14 mat : 00 θ' g θ sin * a / - stack 21 line :

Use U for the start time and IC's

U2 t u D2 U n 1 el te 12 mat U n 2 U cols range el N n - Rkadapt L.ode :

Merge into one solution

U U 1 n range 1 U cols range el U2 stack :

U rows 201

U 1 col 1.0 U 2 col * augment U 1 col 1.0 U 3 col * augment 2 1 sys t u D U 1 2 NN 0 L.vf U 1 col 1.0 U 4 col * augment U 1 col 1.0 U 5 col * augment 2 1 sys t u D U 3 4 2 NN * 0.5*t L.vf

X1 X2 Y2 13 mat x1 x2 y2 13 mat gc U L.At : Fn Tn 12 mat U 1002 L.Frames : t A X2 Y2 augment eval X1 t el 000 X2 t el Y2 t el 32 mat eval 00 + 1000 green 15 mat eval X1 t el 0 X2 t el Y2 t el 0032 mat o augment eval a - 2.2 a * - Tn t el 13 mat eval 51 sys : 31 line

t A

Generalized coordinates

40 A

x1

x2

m.1

gc t r r' θ θ' 15 mat

θ

r

Parameters

y2

m.2

a cm 50

m.1 g 50

m.2 g 50

Initial conditions

ro vo θo ωo 14 mat ic :

1 A

ro cm 50

vo m s / 11033 -^*

θo ° 90

ωo rpm 11033 -^*

<content> Spherical Pendulum </content>

Spherical Pendulum

Gen coords, constants & position

gc θ φ 12 mat r m g 13 mat ∂Vars str2num t θ θ' φ φ' 15 mat

θ φ sph r θ sin φ cos * * r θ sin φ sin * * r 1 θ cos - (* 13 mat :

x y z 13 mat θ φ sph :

x' x d ∂t / :

y' y d ∂t / :

z' z d ∂t / :

Lagrangian

L 12 / m * x' 2^y' 2^+ z' 2^+ (* m g * z * - :

Symplified to

L 12 / m * r 2^* φ' 2^1 θ cos (2^- (* θ' 2^+ (* m r * g * 1 - θ cos + (* + :

Equations of motion

θ''.eq L θ L.eq :

θ''.eq θ sin g - θ cos φ' 2^* r * + (* r /

φ''.eq L φ L.eq :

φ''.eq 2 φ' * θ sin * θ cos * θ' * 1 θ cos 2^- / -

Cyclic coord & conservation

L φ diff 0

p.φ L φ' diff :

p.φ m r 2^* φ' * 1 θ cos 2^- (*

Solutions

N 200 :

te 10 s * :

r m g 13 mat 1 m * 1 kg * g.e 13 mat :

ic 40 ° * 5 rpm * 30 ° * 10 - (rpm * 14 mat :

θo θo' φo φo' 14 mat ic :

t u D θ θ' φ φ' 14 mat u 1 el u 2 el u 3 el u 4 el 14 mat : θ' θ''.eq φ' φ''.eq stack 21 line :

U t u D 0 te 12 mat ic N Rkadapt L.ode :

U 1 col 1.0 U 2 col * augment U 1 col 1.0 U 3 col * augment 2 1 sys t u D U 1 2 NN 0.5*t L.vf U 1 col 0.5 U 4 col * augment U 1 col 1.0 U 5 col * augment 2 1 sys t u D U 3 4 NN 1 L.vf

https://handwiki.org/wiki/Physics:Spherical_pendulum

γ.2 37.5 ° * - 30 ° * pView2 : G sph 0 0.3 π * 20 pR 0 2.1 π * 40 pR pGrid : X Y Z Pφ 14 mat x y z p.φ 14 mat gc U L.At : Fn Tn 12 mat U 1001 L.Frames : A X Y Z augment γ.2 * eval 00 r X t el Y t el Z t el 23 mat γ.2 * eval X t el Y t el Z t el 13 mat γ.2 * o augment eval 00 r 13 mat γ.2 * Tn t el augment eval G 0012 mat 1112 mat pMesh γ.2 * eval 51 sys : 51 line

Φ Pφ 12 mat φ p.φ 12 mat gc U L.At : 11 line

U 1 col Pφ augment <content> Elastic Pendulum </content>

Elastic Pendulum

Gen coords, constants & position

gc x y 12 mat lo m κ g 14 mat ∂Vars str2num t x x' y y' 15 mat

x,y

Lagrangian

L 12 / m * x' 2^y' 2^+ (* m g * y * 12 / κ * x 2^y 2^+ sqrt lo - (2^* + (- :

Equations of motion

x''.eq L x L.eq :

x''.eq x lo - x 2^y 2^+ sqrt + (* κ * m x 2^y 2^+ sqrt * / -

y''.eq L y L.eq :

y''.eq m x 2^y 2^+ sqrt * g * y lo - x 2^y 2^+ sqrt + (* κ * + m x 2^y 2^+ sqrt * / -

Solutions

N 200 :

te 10 s * :

lo m κ g 14 mat 1 m * 1 kg * 20 N m / * g.e 14 mat :

ic 1 m * 0 m s / * 2 m * - 2 m s / * 14 mat :

xo vxo yo vyo 14 mat

t u D x x' y y' 14 mat u 1 el u 2 el u 3 el u 4 el 14 mat : x' x''.eq y' y''.eq stack 21 line :

U t u D 0 te 12 mat ic N Rkadapt L.ode :

U 1 col 1.0 U 2 col * augment U 1 col 1.0 U 3 col * augment 2 1 sys t u D U 1 2 NN 0.1 L.vf U 1 col 1.0 U 4 col * augment U 1 col 1.0 U 5 col * augment 2 1 sys t u D U 3 4 NN 0.5*t L.vf

T X X' Y Y' 15 mat U Cols : Fn Tn 12 mat U 1001 L.Frames : A X Y augment 00 X t el Y t el 22 mat X t el Y t el 12 mat o augment 00 Tn t el augment 41 sys : 31 line

A <content> Double Pendulum </content>

Double Pendulum

Gen coords, constants & position

gc α β 12 mat m.1 m.2 g a b 15 mat ∂Vars str2num t α α' β β' 15 mat

x1 a α sin * :

y1 a α cos * - :

x2 x1 b β sin * + :

y2 y1 b β cos * - :

x1' x1 d ∂t / :

y1' y1 d ∂t / :

x2' x2 d ∂t / :

y2' y2 d ∂t / :

Lagrangian

L 12 / m.1 * x1' 2^y1' 2^+ (* 12 / m.2 * x2' 2^y2' 2^+ (* + m.1 g * y1 * - m.2 g * y2 * - :

Equations of motion

α''.eq L α L.eq :

β''.eq L β L.eq :

L m.1 a * a α' 2^* 1 * 2 g * α cos * + (* m.2 a α cos * α' * b β cos * β' * + (2^a α sin * α' * b β sin * β' * + (2^+ (* + 2 m.2 * g * a α cos * b β cos * + (* + 2 /

α''.eq m.2 b β' 2^β sin α cos * - β cos α sin * + (* β cos α cos * β sin α sin * + (β'' * + (* g α sin * + (* m.1 g * α sin * + a m.1 m.2 + (* / -

β''.eq a α' 2^α sin β cos * - α cos β sin * + (* α cos β cos * α sin β sin * + (α'' * + (* g β sin * + b / -

Push β'' into α''

α''.eq α'' α''.eq β'' β''.eq ≡ at ≡ α'' isol :

Now α''.eq is free of β''. We use that in D

α''.eq m.2 β' 2^β sin α cos * - β cos α sin * + (* b * β cos α cos * β sin α sin * + (a α' 2^* α sin β cos * - α cos β sin * + (* g β sin * + (* - g α sin * + (* m.1 g * α sin * + a m.1 m.2 1 α cos β cos * α sin β sin * + (2^- (* + (* / -

Solutions

N 200 :

te 10 s * :

m.1 m.2 a b g 15 mat 800 g * 1600 g * 120 cm * 250 cm * g.e 15 mat :

ic 30 ° * 5 rpm * 60 ° * 10 rpm * 14 mat :

αo αo' βo βo' 14 mat

t u D α α' β β' 14 mat u 1 el u 2 el u 3 el u 4 el 14 mat : α'' α''.eq : α' α'' β' β''.eq stack eval 31 line :

u x x' 21 mat ≡

Du x' x'' 21 mat ≡

t u D U 34 NN L.vf

U t u D 0 te 12 mat ic N Rkadapt L.ode :

U 1 col 1.0 U 2 col * augment U 1 col 1.0 U 3 col * augment 2 1 sys t u D U 1 2 NN 0.3 L.vf U 1 col 1.0 U 4 col * augment U 1 col 1.0 U 5 col * augment 2 1 sys t u D U 3 4 NN t L.vf

X1 X2 Y1 Y2 14 mat x1 x2 y1 y2 14 mat gc U L.At : Fn Tn 12 mat U 1001 L.Frames : M X1 t el Y1 t el X2 t el Y2 t el 22 mat : A X2 Y2 augment eval 0012 mat M stack eval M o augment eval X1 Y1 augment eval X2 min 0 Tn t el augment eval 51 sys : 41 line

A <content> Particle in a Cone </content>

Particle in a Cone

Gen coords, constants & position

gc r φ 12 mat α m g 13 mat ∂Vars str2num t r r' φ φ' 15 mat

r φ cone r φ cos * r φ sin * r α cot * 13 mat :

x y z 13 mat r φ cone :

x' x d ∂t / :

y' y d ∂t / :

z' z d ∂t / :

Lagrangian

L 12 / m * x' 2^y' 2^+ z' 2^+ (* m g * z * - :

Equations of motion

r''.eq L r L.eq :

φ''.eq L φ L.eq :

L m φ cos r' * r φ sin * φ' * - (2^φ sin r' * r φ cos * φ' * + (2^+ α cot 2^r' 2^* + 2 g * r * α cot * - (* 2 /

r''.eq φ' φ sin φ cos r' * r φ sin * φ' * - (* - φ cos φ sin r' * r φ cos * φ' * + (* + (* - g α cot * + φ cos 2^φ sin 2^+ α cot 2^+ / -

φ''.eq φ sin φ cos r' * r φ sin * φ' * - (* - φ cos φ sin r' * r φ cos * φ' * + (* + r φ' * + (r' * r 2^/ -

Solutions

N 300 :

te 10 s * :

α m g 13 mat 40 ° * 1200 g * g.e 13 mat :

ic 30 ° * 0.3 m s / * 60 ° * 40 rpm * 14 mat :

ro vo φo ωo 14 mat ic :

5 rpm * 0.5236 Hz *

t u D r r' φ φ' 14 mat u 1 el u 2 el u 3 el u 4 el 14 mat : r' r''.eq φ' φ''.eq stack 21 line :

U t u D 0 te 12 mat ic N Rkadapt L.ode :

U 1 col 1.0 U 2 col * augment U 1 col 1.0 U 3 col * augment 2 1 sys t u D U 1 2 NN 0.5*t L.vf U 1 col 0.2 U 4 col * augment U 1 col 1.0 U 5 col * augment 2 1 sys t u D U 3 4 NN 0.1*t L.vf

γ.2 37.5 ° * - 30 ° * pView2 : G cone 0 1.2 U 2 col max * 14 pR 0 2.1 π * 40 pR pGrid : X Y Z 13 mat x y z 13 mat gc U L.At : Fn Tn 12 mat U 1001 L.Frames : A X Y Z augment γ.2 * eval X t el Y t el Z t el 13 mat γ.2 * o augment eval 00013 mat γ.2 * Tn t el augment eval G 0012 mat 1112 mat pMesh γ.2 * eval 51 sys : 51 line

A <content> Simple pendulum </content>

Simple pendulum

Gen coords, constants & position

gc θ a m g θo 14 mat ∂Vars str2num t θ θ' 13 mat

x a θ sin * :

y a θ cos * - :

x' x d ∂t / :

y' y d ∂t / :

Lagrangian

L 12 / m * x' 2^y' 2^+ (* m g * y * - :

L m a * a θ' 2^* 1 * 2 g * θ cos * + (* 2 /

Equations of motion

θ''.eq L θ L.eq : g θ sin * a / -

θ''.eq g θ sin * a / -

Steps

N 200 :

te 10 s * :

Solutions

a m g θo ωo 15 mat 211045 ° * 315 mat :

a m g 13 mat 200 cm * 1200 g * g.e 13 mat :

ic 50 ° * 30 rpm * 12 mat :

θo ωo 12 mat ic :

t u D θ θ' 12 mat u 1 el u 2 el 12 mat : θ' θ''.eq stack 21 line :

U t u D 0 te 12 mat ic N Rkadapt L.ode :

U 1 col 1.0 U 2 col * augment U 1 col 1.0 U 3 col * augment 2 1 sys t u D U 1 2 NN t L.vf

X Y 12 mat x y 12 mat gc U L.At : Fn Tn 12 mat U 1001 L.Frames : A X Y augment eval 00 X t el Y t el 22 mat X t el Y t el o augment X min Y min Tn t el augment 41 sys : 31 line

A <content> Harmonic Oscillator </content>

Harmonic Oscillator

Gen coords, constants & position

gc x κ m 12 mat ∂Vars str2num t x x' 13 mat

x

Lagrangian

L 12 / m * x' 2^* 12 / κ * x 2^* - :

L m x' 2^* κ x 2^* - 2 /

Equations of motion

x''.eq L x L.eq :

x''.eq x κ * m / -

Solutions

N 200 :

te 10 s * :

κ m 12 mat 2 N m / * 1 kg * 12 mat :

ic 10 m * 2 m s / * 12 mat :

xo vo 12 mat

t u D x x' 12 mat u 1 el u 2 el 12 mat : x' x''.eq stack 21 line :

te ic D 2 m s / * 20 m s 2^/ * - 21 mat

U t u D 0 te 12 mat ic N Rkadapt L.ode :

U 1 col 1.0 U 2 col * augment U 1 col 1.0 U 3 col * augment 2 1 sys t u D U 1 2 NN t L.vf

Fn Tn 12 mat U 1001 L.Frames : A 01 X t el 122 mat X t el 1 ■ augment X min 0.5 Tn t el augment 41 sys : 21 line

A <content> Pendulum with Horizontal support </content>

Pendulum with Horizontal support

Gen coords, constants & position

gc x θ 12 mat a m.1 m.2 g 14 mat ∂Vars str2num t x x' θ θ' 15 mat

x2 x a θ sin * + :

y2 a θ cos * - :

x2' x2 d ∂t / :

y2' y2 d ∂t / :

Lagrangian

L 12 / m.1 * x' 2^* 12 / m.2 * x2' 2^y2' 2^+ (* + m.2 g * y2 * - :

L m.1 x' 2^* m.2 x' a θ cos * θ' * + (2^a 2^θ sin 2^* θ' 2^* + (* + 2 m.2 * g * a * θ cos * + 2 /

Equations of motion

x''.eq L x L.eq :

x''.eq a m.2 * θ sin θ' 2^* - θ cos θ'' * + (* m.1 m.2 + / -

θ''.eq L θ L.eq :

θ''.eq θ cos x'' * θ sin g * + a / -

Cyclic coord & conservation

L x diff 0

thus

p.x L x' diff :

p.x x' m.1 * x' a θ cos * θ' * + (m.2 * +

Push x'' into θ''

θ''.eq θ'' θ''.eq ≡ x'' x''.eq ≡ at θ'' isol :

θ''.eq θ sin θ cos a * m.2 * θ' 2^* g m.1 m.2 + (* + (* a m.1 m.2 1 θ cos 2^- (* + (* / -

Solutions

N 200 :

te 10 s * :

a m.1 m.2 g 14 mat 8 m * 1 kg * 3 kg * g.e 14 mat :

ic 0 m * 1 m s / * 70 ° * 0 rpm * 14 mat :

xo vo θo ωo 14 mat

t u D x x' θ θ' 14 mat u 1 el u 2 el u 3 el u 4 el 14 mat : θ'' θ''.eq : x' x''.eq θ' θ''.eq stack 31 line :

U t u D 0 te 12 mat ic N Rkadapt L.ode :

U 1 col 1.0 U 2 col * augment U 1 col 1.0 U 3 col * augment 2 1 sys t u D U 1 2 NN t L.vf U 1 col 1.0 U 4 col * augment U 1 col 1.0 U 5 col * augment 2 1 sys t u D U 3 4 NN t L.vf

X X2 Y2 Px 14 mat x x2 y2 p.x 14 mat gc U L.At : Fn Tn 12 mat U 1001 L.Frames : M X2 t el Y2 t el X t el 022 mat : A X2 Y2 augment eval M M o augment 00 Tn t el augment 41 sys : 41 line

A U 1 col Px augment

Alvaro