<content> Utils </content>

3ff9934e-0071-4049-99d1-c9dc1e9d79f9 475 Alvaro Diaz Falconi 4 5 false false 0 0 decimal

Draghilev's method for Solve Blocks

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

e# a# b# Equrep equrep(e# num2str, a# num2str, b# num2str) concat str2num 11 line :

x# a# Find n# 1 x# length range : f# x# n# el n# augment eval a# 1 findrows : (0 ≡ 0 f# 2 col eval if 21 line :

fs# x# feval fs# num2str (x# num2str) concat str2num 11 line :

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

fs# x# y# z# feval fs# num2str (x# num2str, y# num2str, z# num2str) concat str2num 11 line :

<content> Options </content>

O# o# Options Ω# 01 matrix eval : n# 0 : 12 mat ω# o# sys2mat.1 num2str vectorize eval ω# ω# " strrep : n# n# 1 + : 12 mat g# ω# ≡ strsplit : length 1 > f# 1 O# rows range O# augment eval g# 1 el 2 findrows : (0 ≠ (O# f# 1 el 2 el g# 2 el : O# O# g# 1 el g# 2 el augment 11 O# cols range el stack eval : if Ω# Ω# g# 2 el stack eval : 21 line f# O# ω# 1 findrows : (0 ≠ (Ω# Ω# f# 2 el stack eval : Ω# Ω#'NaN stack eval : if if 21 line for Ω# length 1 > o# cols 1 ≡ Ω# Ω# transpose if Ω# 1 el if 31 line :

<content> Example Options </content>

dsol solver Rkadapt steps 100 stepsmax 0 eps 10^{-7} 42 mat :

dsol dummy Options'NaN

dsol solver Rkadapt steps 100 stepsmax 0 eps 10^{-7} 42 mat

dsol steps dummy eps 31 mat Options 100'NaN 10^{-7} 31 mat

dsol steps dummy eps 13 mat Options 100'NaN 10^{-7} 13 mat

dsol steps solver rkfixed ≡ eps 31 mat Options 100 rkfixed 10^{-7} 31 mat

dsol solver rkfixed steps 100 stepsmax 0 eps 10^{-7} 42 mat

dsol Jacobian J# ≡ Options J#

dsol solver rkfixed steps 100 stepsmax 0 eps 10^{-7} Jacobian J# 52 mat

Check inside a procedure

x f dsol dummy x ≡ Options dsol 21 line :

1 f solver rkfixed steps 100 stepsmax 0 eps 10^{-7} Jacobian J# dummy 162 mat

dsol solver rkfixed steps 100 stepsmax 0 eps 10^{-7} Jacobian J# 52 mat

A# tol# RoundZ B# A# abs tol# > vectorize : A# B# * vectorize eval 21 line :

A# RoundZ B# A# abs TOL > vectorize : A# B# * vectorize eval 21 line :

A# tol# UniqueRows B# A# 1 row : n# A# rows : 12 mat r# 2 : r# n# ≤ r# r# 1 + : flag# 1 : k# 1 : k# r# < k# k# 1 + : A# r# row B# k# row - norme tol# < flag# 0 : break 21 line if for flag# 1 ≡ B# B# A# r# row stack eval : if 31 line for B# 31 line :

A# UniqueRows A# TOL UniqueRows 11 line :

<content> pGrid, pMesh and pSurf </content> <content> pUtil </content>

a b c pDiag a b c 13 mat transpose diag 11 line :

Argument handling utility γ pArgγ γ length 2 ≡ γ 1 el γ 2 el pView2 pArgγ η γ 1 col γ 1 col norme / γ 2 col γ 2 col norme / augment : ζ η 32 el η 21 el * η 32 el η 11 el * - η 11 el η 22 el * sign 1 η 32 el 2^- abs sqrt * stack : η ζ augment γ 13 range 12 range el 12 mat eval 31 line if 11 line :

γ 32 el γ 21 el * γ 32 el γ 11 el * - γ 11 el γ 22 el * sign 1 γ 32 el 2^- abs sqrt * 31 mat

Argument handling utility x n pArgN ans x eval : ans sys2mat.1 stack rows 1 ≡ k 1 n range ans k el x : for 11 line if ans 31 line :

x y z pAxis 000 x 000000 y 000000 z 00073 mat :

M# pAxis A# B# C# 13 mat M# 1 col M# 2 col M# 3 col 13 mat : 1.05 A# abs vectorize max B# abs vectorize max C# abs vectorize max pAxis * eval 21 line :

M γ pAxis pShowAxis M cols 3 ≡ A 0 : k 13 range : 12 mat k m 00013 mat 1.05 M k col max * 3 identity k row * stack γ 13 range 12 range el * eval : A k el polyline k m 11 mat red green blue 31 mat k el solid 2 stack : eval 31 line ToDo if 01 matrix if 11 line :

r g b a rgb hex 0123456789ABCDEF : d 0 : 12 mat c H n 0255 c 0 round 12 mat min 12 mat max : m n 16 mod : 12 mat n m - 16 mod 0 ≡ hex n m - 16 / 1 + 1 substr hex m 1 + 1 substr concat hex m 1 + 1 substr if 21 line : h 255 a * H r H g H b H concat : n h strlen : 12 mat k\0324\ 1 n range d d hex h k\0324\ 1 substr findstr 1 el 1 - (16 n k\0324\ -^* + : for d 232^2 d * 232^1 - ≥ (* - eval 51 line :

C rgb C 1 el C 2 el C 3 el 1 rgb :

C a rgb C 1 el C 2 el C 3 el a rgb :

T n a pCMap C 0 : r 1 n range : c 13 range : 13 mat z r T f 02550 - n 1 - / r 1 - (* + (3 identity z 1 el row * T + a rgb eval 11 line : n 0 ≡ 0 T IsString T Jet ≡ m 0.25 n * 0 round 112 mat max eval : v z v rows 1 matrix : 12 mat x 2551 m range m / * : y 255 m / m 2 / m range (* : o 255 m 1 matrix + : 13 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 : m R length n - : C 1 r el R r m + el G r m + el B r m + el a rgb eval : 71 line T GS ≡ R n a pCMap G n a pCMap + B n a pCMap + rgb RGB T findstr : (1 - ≡ C 1 r el T : C 1 r el rgb r 0 f : if if if zero T 0 ≡ (vectorize c el c augment eval 11 findrows 2 col eval : C 1 r el zero r T f : 21 line if if 31 line :

T x pCMap T x 1 pCMap 11 line :

T pCMap T 16 pCMap 11 line :

z a b CL pClr CL 0 ≡ (CL IsString | CL CL z a - (CL length 1 - b a - / * trunc 1 + el if 11 line :

pCMJet Jet pCMap :

pCMRed R pCMap :

pCMGreen G pCMap :

pCMBlue B pCMap :

<content> pView </content>

EulerRot(v) returns the Euler rotation matrix. At least in some convention. v pEulerRot 1000 v 1 el cos v 1 el sin - 0 v 1 el sin v 1 el cos 33 mat v 2 el cos 0 v 2 el sin 010 v 2 el sin - 0 v 2 el cos 33 mat * v 3 el cos v 3 el sin - 0 v 3 el sin v 3 el cos 000133 mat * eval 11 line :

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 :

View point at az angle from -y axis and el from +x axis. 3x3 Matrix az el pView a az : l el : 12 mat l l 180 ° * + 360 ° * mod 360 ° * + 360 ° * mod 180 ° * - : l 90 ° * > l 180 ° * l - : a a 180 ° * + : 12 mat l 90 ° * - < l 180 ° * - l - : a a 180 ° * + : 12 mat if a a 360 ° * mod 360 ° * + 360 ° * mod : a cos l sin a sin * - l cos a sin * a sin l sin a cos * l cos a cos * - 0 l cos l sin 33 mat eval 51 line :

View point for proyect into 2D plots. 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 ≡

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

γ 32 el γ 21 el * γ 32 el γ 11 el * - γ 11 el γ 22 el * sign 1 γ 32 el 2^- abs sqrt * 31 mat

<content> pGrid </content>

pR(a,b,n) retunrs a vector of equally spaced n+1 points between a and b. a b n pR a b a - n trunc / 0 n trunc range * + eval 11 line :

pR(B,N) returns equally spaced vectors given a Box and a vector with the number of points. B N pR B 1 el B 2 el N 1 el pR B 3 el B 4 el N 2 el pR 12 mat eval 11 line :

Box x.1 x.2 y.1 y.2 22 mat ≡

G pGridBox C G vectorize : C min vectorize C max vectorize augment transpose eval 21 line :

X pBox B 20 matrix : k 1 X cols range B B X k col min X k col max 21 mat augment eval : for B 31 line :

pGrid(F,U,V) returns a length(U)*lenght(V) matrix with [u,v,F(u,v)] FS# X# Y# pGrid X# length 4 ≡ (Y# length 2 ≡ (& U# V# 12 mat X# Y# pR : FS# U# V# pGrid 21 line U# FS# Unknowns : FS# IsDefined ¬ (U# length 1 ≡ (& FS# num2str X# Y# pGrid r# 1 X# length range : c# 1 Y# length range : M# 0 : 13 mat FS# IsDefined FS# IsString num2str 1 ≡ XS# X# num2str vectorize : YS# Y# num2str vectorize : 12 mat str# FS# (#1#,#2#) concat : r# c# F# str# #1# XS# r# el strrep #2# YS# c# el strrep str2num str# #1# XS# r# 1 + el strrep #2# YS# c# 1 + el strrep str2num try 11 line : 31 line r# c# F# FS# r# c# el : if 11 line str# Fxy#(#1#,#2#): FS# num2str concat : str# #1# U# 1 el num2str strrep #2# U# 2 el num2str strrep str2num r c F# X# r# el Y# c# el Fxy# : 31 line if 11 F# GetType matrix ≡ M# r# c# el r# c# F# : M# r# c# el X# r# el Y# c# el r# c# F# 31 mat : if 31 line if 21 line if 11 line :

pTGrid(F,U,V) returns a [G G'] where G is the pGrid and G' is a grid with the middle points of the grid FS X Y pTGrid X length 4 ≡ (Y length 2 ≡ (& U V 12 mat X Y pR : FS U V pTGrid 21 line A pr n 1 A length 1 - (range : A' 0 : 12 mat A' n el 0.5 A n el A n 1 + el + (* : A' A length el A A length el : 31 line : FS X Y pGrid FS X pr Y pr pGrid 12 mat eval 21 line if 11 line :

<content> pMesh </content>

pMesh returns a matrix suitable for 3D plots. G# h pMesh tri G# length 2 ≡ eval : G G# 1 el G# 1 el G# 1 el cols col augment eval : G' G# 2 el eval : 12 mat G G# eval : if nu G rows 1 - : nv G cols 1 - : M 03 matrix eval : r' 0 : 14 mat nh h cols : nh' nh tri + : r.nh 1 nh 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.nh + : r' r' 1 + : 13 mat M k c el G h 1 k r - el u + h 2 k r - el v + el c el : tri v nv ≡ M r nh' + c el M r nh + c el : M r nh' + c el G' r' el c el : if if 31 line for r M rows : k r r.nv + : 12 mat M k c el M u 1 - (nh' nv * nv + 1 - (* nh' v 1 - (* + k r - 1 - (nh' * - tri - c el : 31 line for M 51 line :

Mesh of a Grid G# pMesh G# 001100110025 mat pMesh 11 line :

Triangular Mesh of a Grid G# pTMesh G# 001010110025 mat pMesh 11 line :

Utill: pMesh(FS,X,Y) Mesh from the function or        pMesh(M,CL,γ) Coloring a Mesh for a XYPlot M# γ CL pMesh M# IsString M# IsDefined 1 ≡ num2str (¬ (| M# γ CL pGrid pMesh γo γ2 12 mat γ pArgγ : M M# eval : M2 M γ2 * eval : k 2 M rows (range : L 0 : 14 mat a b 12 mat M 3 col min M 3 col max 1010 -^+ 12 mat eval : M pBox γ pAxis L k el M2 k row M2 k 1 - row stack M k 3 el a b CL pClr solid 10 pLine : stack 11 mat eval 41 line if 11 line :

<content> pSurf </content>

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

X pCenter X 1 col sum X 2 col sum X 3 col sum 13 mat X rows / eval 11 line :

0.5 M 1 c el M ν c el + (* eval

nG G# γ CL CS h wm LS d pSurf c 13 range : n 0 : S 0 : N 0 : MG 03 matrix eval : γo γ2 12 mat γ pArgγ : 16 mat s 1 nG range G# s el cols 2 ≡ MG MG G# s el 1 el G# s el 2 el stack eval pGridBox stack eval : MG MG G# s el eval pGridBox stack eval : if 11 line for MG MG pBox : a b 12 mat MG 1 d el MG 2 d el 1010 -^+ 12 mat eval : 12 mat Vo 100 MG normi max + 13 matrix + (γo * eval : s 1 nG range tri G# s el cols 2 ≡ eval : G G' 12 mat G# s el 1 el G# s el 2 el 12 mat eval : G G# s el eval : if nu nv 12 mat G rows G cols 12 mat eval : nh mh ν w 14 mat h cols h max h cols 312 mat min wm CL 0 ≠ (* 14 mat eval : M 0 : r 1 nh range : nj 1 tri 0 nh 2 - (range * + eval : 13 mat u 1 nu mh - (range v 1 nv mh - (range j nj tri M G' u v el G u h 1 j el + v h 2 j el + el G u h 1 j 1 + el + v h 2 j 1 + el + el augment transpose eval : M r c el G u h 1 r el + v h 2 r el + el c el : if Mo M pCenter : N n n 1 + : el Mo Vo - norme eval : S n el M γ2 * Mo d el a b CL s el pClr LS w Mo d el a b CS s el pClr pLine eval : 41 line for for for 51 line for N 1 n range N augment eval 2 csort 1 col reverse eval : MG γ2 pAxis S N el stack 11 mat eval 71 line :

CL 0 ≡ (CL IsString | zo CL CL : nL CL length eval : zo CL CL zo a - (nL 1 - b a - / * trunc 1 + el 11 line : 21 line if CS 0 ≡ (CS IsString | zo CS CS : nS CS length eval : zo CS CS zo a - (nS 1 - b a - / * trunc 1 + el 11 line : 21 line if 12 mat

0.5 M 1 c el M ν c el + (* eval

Surface plot M γ CL CS pSurf M sys2mat.1 M ≡ G M 021 mat : n 1 : 12 mat G M : n G sys2mat.1 stack rows : 12 mat if n G γ CL n pArgN CS n pArgN 001100110025 mat 1 solid 3 pSurf 21 line :

This shows or hide the axis.

pShowAxis 1 :

Some colormaps and viewpoints

CM.T 01923213 mat 641 pCMap :

γ.1 120 ° * 30 ° * pView2 :

CM.R R 641 pCMap :

γ.2 30 ° * 60 ° * pView2 :

CM.GS GS 641 pCMap :

γ.3 37.5 ° * - 30 ° * pView :

CM.J Jet 641 pCMap :

γ 120 ° * 30 ° * pView :

TOL 107 -^:

CM CM.J :

γ.2 120 ° * 30 ° * pView2 :

E# λ# N# Z# EQ# Do# nDM m kg s K A cd 16 mat 11111116 mat : ξ# τ# Clear K# 0 : nk# 0 : 13 mat s# E# sys2mat.1 num2str vectorize eval g# s# ≡ strsplit : g'# s# ≈ strsplit : 12 mat Draghilev g# 1 el 1 findrows 0 ≡ nk# nk# 1 + : g# length 1 > K# nk# 1 el EQ : K# nk# 2 el g# 1 el str2num g# 2 el str2num - : 12 mat g'# length 1 > K# nk# 1 el EG : K# nk# 2 el g'# 1 el str2num g'# 2 el str2num - : 12 mat K# nk# 1 el EB : K# nk# 2 el s# str2num : 12 mat if 21 line Draghilev s# Options if 21 line for Opt# DefaultGuess dsolver eps Speed 14 mat : vg# rk# ε# vo# 14 mat Draghilev Opt# Options : x# K# Unknowns : n# x# length : 12 mat k# 1 x# length range K# K# x# k# el u# k# el Equrep : for k# u# u#_k# num2str concat str2num 11 line : EQ# K# EQ 1 findrows 2 col : vg# vg# str2num n# 1 matrix + eval : f# K# EG 1 findrows : num2str 0 ≠ EG# f# 2 col : u# EG_F# EG# eval : 12 mat vg# vg# 01030 -^ξ# EG_F# al_nleqsolve : 21 line if rg# Draghilev OptimizGuess Options str2num : (0 > Z# num2str ≡ rg# 1 n# EQ# length 12 mat min range eval : rg# 1 n# EQ# length rg# 13 mat min range eval : if u# F1# EQ# rg# el : vg# vg# 01030 -^ξ# F1# al_nleqsolve : 12 mat 21 line if EQ# length n# < EQ# n# el n# 1 + u# : if u# F# EQ# : X# 0 : Do# 0 : Do# 2 n# * 1 + el 0 : 13 mat k# 12 n# * 1 + (range : X# k# el k# u# : 12 mat X1# X# 1 n# range el : X2# X# n# 2 + (2 n# * 1 + (range el : 12 mat F2# X2# F# : k# 1 n# range : F2# k# el n# 1 + u# F2# k# el * : Y# X1# F# F2# - : 14 mat k# 1 n# 1 + (range Xo# X1# : k# n# 1 + (< Xo# k# el n# 1 + u# : if Do# k# el 1 - (k# n# 1 + ≠^Y# Xo# Jacob det * : 31 line for k# 12 n# * 1 + (range EQ# Do# 12 mat EQ# Do# 12 mat k# u# u# k# el Equrep : for u# F# EQ# : Draghilev Normalize Options 1 ≡ Dn# Do# Do# norme / : t# u# D# Dn# : 12 mat t# u# D# Do# u# try 11 line : if t# u# D# Unknowns trace Draghilev dsolver Options Internal ≡ h# λ# N# 1 + / : v# vg# 1 vg# stack eval : RK# v# transpose : r# 1 v# length range : 14 mat v# κ# DK# h# 0 v# D# * eval h# κ# * eval try 11 line : c# 2 N# 1 + (range c# trace k.1 v# v# DK# : k.2 v# k.1 + k.1 DK# : k.3 v# 0.25 k.1 * + 0.25 k.2 * + k.2 DK# : v# v# k.1 4 k.3 * + k.2 + 6 / + eval : RK# c# r# el v# r# el : 12 mat 51 line for RK# 0 RK# augment eval : 41 line RK# rk# (eval(stack(vg#, vo#,vg#)),0,λ#,N#,D#(τ#,ξ#)) concat str2num : if Z# num2str ≠ V# RK# n# 2 + col : Z# 0 : Roots# 0 : r# 0 : 14 mat k# 1 : k# V# length < k# k# 1 + : V# k# 1 + el V# k# el * 0 ≤ Roots# r# r# 1 + : el k# : k# k# V# k# el 0 ≡ (+ : 12 mat if 11 line for P# ρ# I# P# ρ# el V# ρ# el V# ρ# 1 + el V# ρ# el - / P# ρ# 1 + el P# ρ# el - (* - eval : r# 0 > (Draghilev OptimizRoots Options all ≡ (& k# 1 r# range : c# 1 n# range : Z# k# c# el RK# c# 1 + col Roots# k# el I# : 13 mat Z# k# c# el Z# k# row transpose 0 ε# str2num ξ# F# al_nleqsolve eval c# el : Z# Z# RoundZ eval : 31 line Z# 0 : if 41 line if RK# 1 RK# rows range 2 n# 2 + (range el 261 line :

E# λ# N# nDM E# λ# N# nDM 11 line :

E# λ# N# Z# nDM E# λ# N# Z# nDM 11 line :

E# λ# N# Z# EQ# nDM E# λ# N# Z# EQ# nDM 11 line :

T γ.2 Ro U Eq pDM sz Draghilev PointSz Options str2num : u# E# Eq : x y F# x y 12 mat E# 1 el : x y G# x y 12 mat E# 2 el : 13 mat R# Ro 0 augment γ.2 * : XY2 U γ.2 * : XYo U 1 col U 2 col augment : 13 mat T 2 ≡ x y F# x y G# XYo Ro 1 Ro rows range num2str vectorize sz black augment Ro o augment XYo 1 row . sz red augment XYo 1 row Start sz black augment XYo XYo rows row . sz red augment XYo XYo rows row End sz black augment 91 sys 11 line T 3 ≡ XYo 0 augment γ.2 * XY2 U pAxis γ.2 * R# o sz black augment R# 1 R# rows range num2str vectorize sz black augment XY2 1 row . sz red augment XY2 1 row Start sz black augment XY2 XY2 rows row . sz red augment XY2 XY2 rows row End sz black augment 101 sys 11 line if 41 line :

f# g# Box nG γ.2 Roo U Eq pDM_LM sz Draghilev PointSz Options str2num : δy# f#(x,y): f# num2str (x,y) concat str2num : δy# g#(x,y): g# num2str (x,y) concat str2num : Roo 0 ≡ Ro 00013 mat : Ro Roo : if u# E# Eq : x y F# x y 12 mat E# 1 el : x y G# x y 12 mat E# 2 el : 13 mat X' U 1 col : Y' U 2 col : Z X' Y' f# vectorize : 13 mat U3 U 1 col U 2 col U 4 col augment : 11 mat XYo U 1 col U 2 col augment : 11 mat x y f_g x y g# 0 ≤ x y x y f# 13 mat u g u 1 el u 2 el g# 11 mat : xy x y 21 mat uu g al_nleqsolve : xy 1 el xy 2 el xy 1 el xy 2 el f# 13 mat 31 line if 11 line : r 1 Ro rows range : Ro r 3 el Ro r 1 el Ro r 2 el f# : 12 mat r 1 XYo rows range : XYo r 3 el XYo r 1 el XYo r 2 el f# : 12 mat X Y 12 mat Box 2 nG * pR : r 1 X length range : c 1 Y length range : Gfg 0 : 14 mat Gfg r c el X r el Y c el f_g : f# Box nG pGrid pMesh γ.2 * Gfg pMesh γ.2 * U3 γ.2 * U3 pAxis γ.2 * XYo γ.2 * Ro γ.2 * . sz green augment Ro γ.2 * 1 Ro rows range num2str vectorize sz black augment U3 1 row γ.2 * o sz black augment U3 U3 rows row γ.2 * o sz black augment U3 1 row γ.2 * Start sz black augment U3 U3 rows row γ.2 * End sz black augment 111 sys 141 line :

Options Draghilev dsolver rkfixed steps 100 nsolver al_nleqsolve stepsmax 0 eps 10^{-7} DefaultGuess 1 UpdateGuess 1 MaxIters 50 Speed 1 OptimizGuess 1 PointSz 8 Normalize 1 OptimizRoots all 132 mat :

<content> Draghilev's method examples </content>

Sintaxis

E λ N nDM

Apply the Draghilev's method to the system E for a path length λ for N steps of a numerical ode solver.

E λ N Z nDM

Uses Z for store the roots of the system E

E λ N Z EQ Δ nDM

Store the equations in EQ and the Draghilev system in Δ

Parametrizing plane curves

<content> Parametrizing plane curves </content>

Examples

nDM try to find the starting point near the provided guess values

a 3 :

b 2 :

f 41 x 2^9 / - (* sqrt :

U 15.8 100 nDM : U 1 col U 2 col augment

x 2^a 2^/ y 2^b 2^/ + 1 ≡

f' f x diff eval :

x 0.5 ≈

y 1 ≈

41 f' 2^+ sqrt x 0 a 0.0001 - int * 15.8006

Arc length of the ellipse

L π 3 a b + (* 3 a * b + (a 3 b * + (* sqrt - (* : 15.87

U k row U k 1 - row - norme k 2 U rows sum 15.8

This is: the argument λ is the path lengh

U.1 10 100 nDM : U.2 12 100 nDM :

4 x 2^* 7 y 5^* - 12 x cos * ≡

4 x 4^* 2 y 4^* + 45 x cos * ≡

x 4 - ≈

y 1 ≈

x 4 - ≈

y 1 ≈

U.1 1 col U.1 2 col augment U.2 1 col U.2 2 col augment

By default nDM uses rkfixed

Draghilev dsolver Options rkfixed

U.2 16 100 nDM : U.1 12 - 100 nDM :

x 1 - (x 2^y 2^+ (* 4 x 2^* + 0 ≡

x 3^y 3^+ 3 x * y * ≡

x 1 ≈

y 4 - ≈

x 4 - ≈

y 2 ≈

dsolver dn_AdamsMoulton ≡

dsolver al_rkckadapt ≡

U.1 1 col U.1 2 col augment U.2 1 col U.2 2 col augment

Root finding

<content> Viacheslav example 1 </content>

Example

A similar scheme can be used to solve a system of equations. The green curve shows the path of the variable introduced by the method.

U 40 200 Ro Eq nDM :

Ro rows 12

1 U 3 col - normi 3.4819

x 2^x y * - y 2^+ 1 ≡

U k row U k 1 - row - norme k 2 U rows sum 39.48

5 x 2^* sin 4 y 2^* sin + 0 ≡

x 0 ≈

y 1 ≈

2 γ.2 Ro U Eq pDM 1 1 line 3 γ.2 Ro U Eq pDM 1 1 line

These seed values are more efficient in the search for roots, since the path traveled is shorter.

U 20 100 Ro Eq nDM :

x 2^x y * - y 2^+ 1 ≡

5 x 2^* sin 4 y 2^* sin + 0 ≡

Ro rows 12

1 U 3 col - normi 2.1573

U k row U k 1 - row - norme k 2 U rows sum 19.57

x 1 ≈

y 1 ≈

2 γ.2 Ro U Eq pDM 1 1 line 3 γ.2 Ro U Eq pDM 1 1 line <content> Viacheslav example 2 </content>

Example

The first equation matter: nDM uses it to set the initial point of the method by resolving it with the seed values.

U 20 - 400 Ro Eq nDM :

Ro rows 16

U 1.4843 1.4375 1 1.4879 1.4365 0.9501 1.4913 1.4355 0.9003 1.4945 1.4346 0.8504 43 mat

4 y 1 - (2^* x 1 - (2^+ 1 ≡

y x 25 x * sin * ≡

x 1.5 ≈

y 1.5 ≈

2 γ.2 Ro U Eq pDM 1 1 line 3 γ.2 Ro U Eq pDM 1 1 line U 20 - 400 Ro Eq nDM :

U 1.5595 1.4981 1 1.5599 1.5018 1.0499 1.5602 1.5055 1.0997 1.5605 1.5092 1.1496 43 mat

y x 25 x * sin * ≡

4 y 1 - (2^* x 1 - (2^+ 1 ≡

x 1.5 ≈

y 1.5 ≈

2 γ.2 Ro U Eq pDM 1 1 line 3 γ.2 Ro U Eq pDM 1 1 line <content> Viacheslav example 3 </content>

Example

U 3 100 Ro Eq nDM :

3 x * y z * cos - 0.5 ≡

Ro 0.5 0 0.5236 - 0.4981 0.1996 - 0.5288 - 23 mat

x 2^81 y 0.1 + (2^* - z sin + 1.06 + 0 ≡

20 z * e x y * -^+ 13 / 3 - 10 π * + (* + 0 ≡

T 1 U rows range :

XYZ U 1 col U 2 col U 3 col augment :

x 1 ≈

y 1 ≈

z 1 ≈

T U 1 col augment T U 2 col augment T U 3 col augment T U 4 col augment 4 1 sys XYZ γ.2 * XYZ 1 row γ.2 * . 12 red augment Ro γ.2 * o augment 3 1 sys <content> Viacheslav example 4 </content>

Example

U 19 200 Ro Eq nDM :

U Δo Eq 13 mat Info : u# E# Eq : x y F# x y 12 mat E# 1 el : x y G# x y 12 mat E# 2 el : 14 mat

x 2^y x * - y 2^+ 1 ≡

Ro rows 12

5 x 2^* sin 4 y 2^* sin + 0 ≡

U 1 row 22113 mat

x 1 ≈

y 1 ≈

U 11 el U 12 el F# 0

2 γ.2 Ro U Eq pDM 1 1 line 3 γ.2 Ro U Eq pDM 1 1 line <content> Viacheslav example 5 </content>

Example

U 120 500 Ro Eq nDM :

x 2^y 2^+ 8 ≡

x y * sin x y * exp sin * 0 ≡

Ro rows 26

x 3 ≈

y 0.5 - ≈

2 γ.2 Ro U Eq pDM 1 1 line 3 γ.2 Ro U Eq pDM 1 1 line

Minimizing functions

<content> Rosenbrock </content>

Example

a 3 : b 100 : 21 sys

Extrema of

x y f a x - (2^b y x 2^- (2^* + :

U 10 - 100 Ro Eq nDM :

Xo Yo 12 mat Ro 1 col Ro 2 col 12 mat :

x y f x diff 0 ≡

x y f y diff 0 ≡

Zo Xo Yo f vectorize RoundZ : 011 mat

x 2 ≈

y 2 ≈

OptimizGuess 0 ≡

Ro 3912 mat

Zo 011 mat

Rfo Xo Yo Zo augment : f 0 4 0 10 2 2 mat 3 12 12 2 1 mat * pMesh γ.2 * Rfo γ.2 * . 8 green augment Rfo pAxis γ.2 * 3 1 sys 2 1 line 2 γ.2 Ro U Eq pDM 1 1 line 3 γ.2 Ro U Eq pDM 1 1 line <content> Himmelblau </content>

Example

Extrema of

x y f x 2^y + 11 - (2^x y 2^+ 7 - (2^+ :

U 30 - 100 Ro Eq nDM :

Xo Yo 12 mat Ro 1 col Ro 2 col 12 mat :

x y f x diff 0 ≡

x y f y diff 0 ≡

Zo Xo Yo f vectorize RoundZ :

x 7 - ≈

y 7 - ≈

OptimizGuess 0 ≡

Ro 3.78 - 3.28 - 3.07 - 0.08 - 2.81 - 3.13 0.09 2.88 0.27 - 0.92 - 0.13 - 1.95 - 3.58 1.85 - 3.39 0.07 3292 mat

Zo 0 104.02 0 67.72 181.62 178.34 0 13.31 091 mat

Rfo Xo Yo Zo augment : f 4.5 1 - 1 1 - 1 2 2 mat * 3 12 12 2 1 mat * pMesh γ.2 * Rfo γ.2 * . 8 green augment Rfo pAxis γ.2 * 3 1 sys 2 1 line

f.max Zo max : 181.62

Ro Zo f.max - RoundZ 0 Find 12 range el 0.27 - 0.92 - 12 mat

f.min Zo min : 0

Ro Zo f.min - RoundZ 0 Find 12 range el 3.78 - 3.28 - 2.81 - 3.13 3.58 1.85 - 3242 mat

2 γ.2 Ro U Eq pDM 1 1 line 3 γ.2 Ro U Eq pDM 1 1 line

Lagrange multipliers

<content> Lagrange multipliers 1 </content>

Example

Extrema of f restricted to g=0

λo xo yo Clear 1

x y f x y + :

x y g x 2^y 2^+ 1 - :

L x y f λ x y g * + :

U 5 200 Ro Eq nDM : f g 1.5 1 - 1 1 - 1 2 2 mat * 12 12 2 1 mat γ.2 Ro U Eq pDM_LM 1 1 line

L λ diff 0 ≡

L x diff 0 ≡

L y diff 0 ≡

λ 1 ≈

x 1 - ≈

y 1 - ≈

OptimizGuess 0 ≡

Ro 0.71 - 0.71 - 0.71 0.71 0.71 0.71 - 23 mat

<content> Lagrange multipliers 2 </content>

Example

Extrema of f restricted to g=0

λo xo yo Clear 1

x y f x y + (2^:

x y g x 2^y 2^+ 1 - :

L x y f λ x y g * + :

L λ diff 0 ≡

L x diff 0 ≡

L y diff 0 ≡

xo yo Ro Eq U E 10100 Ro Eq nDM 11 line :

U 34 R1 0 U :

U 34 - R2 Eq U :

λ 1 ≈

x xo ≈

y yo ≈

Ro R1 R2 stack UniqueRows :

OptimizGuess 0 ≡

Xo Ro 1 col :

Yo Ro 2 col :

Zo Xo Yo f vectorize RoundZ :

Ro 0.71 - 0.71 0 0.71 - 0.71 - 2 - 0.71 0.71 - 0 0.71 0.71 2 - 43 mat

Zo 020241 mat

f g 1.5 1 - 1 1 - 1 2 2 mat * 12 12 2 1 mat γ.2 Ro U Eq pDM_LM 1 1 line

f.max Zo max : 2

Ro Zo f.max - RoundZ 0 Find 12 range el 0.71 - 0.71 - 0.71 0.71 22 mat

f.min Zo min : 0

Ro Zo f.max - RoundZ 0 Find 12 range el 0.71 - 0.71 - 0.71 0.71 22 mat

<content> Lagrange multipliers 3 </content>

Example

Extrema of f restricted to g=0

x y f x 2^y * :

x y g x 2^y 2^+ 3 - :

L x y f λ x y g * + :

U 10 100 Ro Eq nDM : f g 2 1 - 1 1 - 1 2 2 mat * 12 12 2 1 mat γ.2 Ro U Eq pDM_LM 1 1 line

L λ diff 0 ≡

L x diff 0 ≡

L y diff 0 ≡

λ 1 ≈

x 1 ≈

y 1 ≈

OptimizGuess 0 ≡

Ro 0 1.7321 0 1.4142 - 11 - 1.4142 - 1 - 133 mat

<content> Lagrange multipliers 4 </content>

Minimize f subject to g.1, g.2 ≡ 0

a b c 13 mat 25 sqrt 513 mat :

g.1 x 2^a 2^/ y 2^b 2^/ + z 2^c 2^/ + 1 - : g.2 x y + z - : 21 sys

Example

f x 2^y 2^+ z 2^+ :

L f λ.1 g.1 * + λ.2 g.2 * + :

Eo E 15 range el : Uo Eo Unknowns : 12 mat Eo Uo Solve 21 line

λo s R U E s 20 R nDM 11 line :

L x diff 0 ≡

L y diff 0 ≡

L z diff 0 ≡

U 2 - 1 R1 U :

U 11 - R2 U :

25 sqrt * 19 sqrt / 35 sqrt * 19 sqrt / 52 nthroot (3^19 sqrt / 522 nthroot (5^* 323 sqrt / 35646 sqrt / 5646 sqrt / 23 mat 1.03 1.54 2.56 1.57 1.38 0.2 23 mat

L λ.1 diff 0 ≡

L λ.2 diff 0 ≡

Ro R1 R2 stack : 1.03 1.54 2.56 10 - 3.08 1.57 - 1.38 0.2 - 4.41 - 0.32 - 25 mat

x 0.5 ≈

y 1 ≈

z 1 ≈

f Ro 113 range el at 10

f Ro 213 range el at 4.4118

λ.1 λo ≈

λ.2 λo ≈

The option OptimizGuess ≡ 5 forces to nDM uses the five first equations for correct the guess values.

OptimizGuess 5 ≡

<content> Lagrange multipliers 5 </content>

Example

Extrema of f restricted to g=0 and h<0

f x y + z + :

h.1 y 1 - (2^z 2^+ 1 - : h.2 x 2^y 1 - (2^+ z 2^+ 3 - : 21 sys

L f μ.1 h.1 * + μ.2 h.2 * + :

xo s R U E s 20 R nDM 11 line :

L x diff 0 ≡

L y diff 0 ≡

L z diff 0 ≡

U 21 R1 U :

U 2 - 1 - R2 U :

L μ.1 diff 0 ≡

L μ.2 diff 0 ≡

Ro R1 R2 stack : 1.41 1.71 0.71 0.35 - 0.35 - 1.41 - 0.29 0.71 - 0.35 0.35 25 mat

2 sqrt 112 sqrt / + 12 sqrt / 122 sqrt * / - 122 sqrt * / - 2 sqrt - 112 sqrt / - 12 sqrt / - 122 sqrt * / 122 sqrt * / 25 mat 1.4142 1.7071 0.7071 0.3536 - 0.3536 - 1.4142 - 0.2929 0.7071 - 0.3536 0.3536 25 mat

x xo ≈

y 1 ≈

z 2 ≈

μ.1 1 ≈

μ.2 1 ≈

f Ro 113 range el at 3.8284

f Ro 213 range el at 1.8284 -

OptimizGuess 6 ≡

The option OptimizGuess ≡ 6 forces to nDM uses the six first equations for correct the guess values.

Parametrizing space curves

<content> Parametrizing space curves </content>

Examples

U 35 - 200 - eq nDM :

X 2.01 - 9300 pR :

Δ 916 X 3^+ (* X 4^- vectorize :

X Δ augment

x 2^9 / z 2^4 / + y ≡

x 2^x y * - 4 z * ≡

Z.1 212 - Δ sqrt + (* 3 X * / vectorize :

Z.2 212 Δ sqrt + (* 3 X * / vectorize - :

x 7 ≈

y 7 ≈

z 4 - ≈

Y.1 X 2^9 / Z.1 2^4 / + vectorize :

Y.2 X 2^9 / Z.2 2^4 / + vectorize :

dsolver dn_AdamsMoulton ≡

Δf y x 2^9 / - : x y f.11 2 Δf sqrt * Δf 0 ≥ (* - : x y f.12 2 Δf sqrt * Δf 0 ≥ (* : x y f.2 x 2^x y * - 4 / : Box 2 - 92 - 922 mat : nG 2121221 mat * : 61 line

U 1 col U 2 col U 3 col augment γ.2 * X Y.1 Z.1 augment γ.2 * 0.1 + X Y.2 Z.2 augment γ.2 * 0.1 - f.11 Box nG pMesh f.12 Box nG pMesh stack γ.2 * f.2 Box nG pMesh γ.2 * X Y.1 Z.1 augment pAxis γ.2 * U 1 1 3 range el γ.2 * . 20 red augment U U rows 1 3 range el γ.2 * . 20 blue augment 8 1 sys <content> Parametrizing space curves </content>

Example

Disabling normalization of the Draghilev derivatives. Example from Ber7.

https://en.smath.com/forum/yaf_postsm11927_Method-A-B--Draghilev-and-animation-spatial-mechanisms.aspx#post11927

U 0.00864 275 nDM :

x.1 0.5 e 1 -^* y.1 * 7 y.1 * sin * - (2^y.1 0.5 9 z.1 * sin * - (2^+ z.1 0.5 11 x.1 * sin * - (2^+ 12 ≡

x.1 2^0.25 y.1 * + 0.2 z.1 * - 0 ≡

U 1 U rows range 1 3 5 1 3 mat el γ.2 * U 1 U rows range 2 4 6 1 3 mat el γ.2 * 2 1 sys

x.2 4 + (2^y.2 2^+ 9 - 0 ≡

y.2 z.2 4^- 0 ≡

x.1 x.2 - (2^y.1 y.2 - (2^+ z.1 z.2 - (2^+ 20 - 0 ≡

OptimizGuess 0 ≡

Normalize 0 ≡

x.1 0.683 - ≈

y.1 0.476 ≈

z.1 2.926 ≈

x.2 4.005 - ≈

y.2 3.0004 ≈

z.2 1.31607 ≈

Intercepting surfaces

<content> Intercepting surfaces </content>

Example

R 2 : r 1 : 21 sys

x y φ x 2^x y * - x - y + 41 x cos 2^+ (* / :

u v F R r u cos * + (v cos * R r u cos * + (v sin * r u sin * 31 mat :

u v Clear 1

U 20 300 nDM :

Plot code κ 3 : ν 12 : n# 0 : j# 13 range : 14 mat W U : r# 1 W rows ν / trunc range c# 1 ν range n# n# 1 + : W n# j# el U n# 1 el U n# 2 el F j# el : GU r# c# el W n# row transpose : 31 line for for Gφ φ 3 - 33 - 322 mat κ 12821 mat * pGrid : GF F 02 π * 0.5 π * - 0.5 π * 22 mat κ 12821 mat * pGrid : P.1 Gφ pMesh γ.2 * GF pMesh γ.2 * W γ.2 * 2 W * pAxis γ.2 * 41 sys : P.2 2 W * pAxis γ.2 * Gφ GF GU 31 sys γ.2 black black white 31 sys CM.J CM.T white 31 sys pSurf 21 sys : 71 line

P.1 1 1 line

u v F 1 el u v F 2 el φ u v F 3 el ≡

u 1.5 ≈

v 1.5 ≈

dsolver dn_AdamsMoulton ≡

P.2 1 1 line

Alvaro

4 appVersion 1.0.8348.30405