<content> dsol </content>

f9c48a99-5a4a-4395-8b73-738b8f1c8321 49 4 5 false false false 0 3 decimal

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

Campana di Pompelmo

<content> dsol </content> <content> Utils </content>

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

a b n pR a b a - n / 0 n range * + eval 11 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 :

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 Ω# rows 2 > 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 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

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 non-linear BVP y'' = f(x,y,y') in a < x < b with M = [α1 β1 c1; α2 β2 c2] where     α1*y(a) + β1*y'(b) = c1   and  α2*y(a) + β2*y'(b) = c2using central differences. fs# x Y M n ε bvp.2 x# y# y'# f fs# x# y# y'# feval : 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 dsol MaxIters Options str2num 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 dsol MaxIters Options str2num < X Y Y' augment eval MaxIters  rached. error if 91 line :

dsol(DE,X,tf,N) try to solve the Differential Equation system DE for the functions in X. Ev# Yv# TF# N# dsol m kg s K A cd 16 mat 11111116 mat : u# v# t# EQ_F# IG_F# EQ_F# RK_D# Clear ky# 0 : kn# 0 : kv# 0 : 14 mat ys# 0 : U# 0 : V# 0 : K# 0 : NK# 0 : IC# 03 matrix eval : 16 mat M# 0 : ES# Ev# num2str : ES'# ES#' strsplit rows : 13 mat BIC# 0 : bic# 0 : too# : 13 mat s# Yv# sys2mat.1 num2str vectorize eval a# b# 12 mat s# (findstr 1 el s#) findstr 1 el 12 mat : ts# s# a# 1 + b# a# - 1 - substr : Assuming all ≡ 12 mat ys# ky# ky# 1 + : el s# 1 a# 1 - substr : k# 1 ES'# range ES# s# s# ('(strrep : findstr 1 - ≠ M# ky# el k# : if for 41 line for s# Ev# sys2mat.1 num2str vectorize eval g# s# ≡ strsplit : g'# s# ≈ strsplit : 12 mat dsol g# 1 el 1 findrows 0 ≡ NK# NK# 1 + : s# (ts#) concat findstr 1 - ≡ a# b# 12 mat g# 1 el (findstr 1 el g# 1 el) findstr 1 el 12 mat : to# s# a# 1 + b# a# - 1 - substr : too# ≡ (too# to# ≠ (| BIC# bic# bic# 1 + : el to# : too# to# : 12 mat if g'# length 1 ≡ ico# s# 1 a# 1 - substr : S# ico#' strrep (ts#) concat strrep : f# ys# S# Find : (0 ≠ IC# IC# f# 1 el ico#' strsplit rows g# 2 el 13 mat stack : if K# NK# 1 el IC : K# NK# 2 el g# 1 el str2num g# 2 el str2num - : bvp 13 mat 41 line K# NK# 1 el IG : K# NK# 2 el g'# 1 el str2num g'# 2 el str2num - : 12 mat if 41 line g# length 1 ≡ K# NK# 1 el EQ : K# NK# 2 el s# str2num : 12 mat K# NK# 1 el EQ : K# NK# 2 el g# 1 el str2num g# 2 el str2num - : 12 mat if if 21 line dsol s# Options if 21 line for ic# 1000 IC# 1 col * IC# 2 col + IC# 3 col augment : ic# ic# 1 csort 2 col eval : ic# ic# str2num vectorize : 12 mat n# 1 ky# range kn# kn# 1 + : a# ys# n# el (ts#) concat : 12 mat K# K# a# u# kn# el Equrep : m# 1 : m# M# n# el ≤ m# m# 1 + : a# a# ('(strrep : m# M# n# el 1 - ≤ kn# kn# 1 + : kv# kv# 1 + : V# kv# el u# kn# el : 13 mat K# K# a# u# kn# el Equrep : 21 line if 21 line for kv# kv# 1 + : U# n# el v# n# el : V# kv# el v# n# el : 13 mat K# K# a# v# n# el Equrep : 51 line for EQ# K# EQ 1 findrows 2 col : v# EQ_F# EQ# eval : 12 mat bic# 1 > IC# K# IC 1 findrows 2 col : u# FI# IC# : 12 mat M# u1# u2# 21 mat FI# u1# diff vectorize u1# u2# 21 mat FI# u2# diff vectorize 0021 mat FI# - augment eval : dsol dsolver Options lbvp.2 ≡ u# v# FE# EQ# 1 el : dv# u1# u2# 21 mat v1# 11 mat FE# v1# diff : 12 mat L# u1# u2# 21 mat v1# 11 mat FE# u2# diff dv# / u1# u2# 21 mat v1# 11 mat FE# u1# diff dv# / 0021 mat 011 mat FE# dv# / - 13 mat : lf#(ts#): L# num2str concat str2num v1# lf# BIC# str2num vectorize M# N# lbvp.2 41 line dsol dsolver Options bvp.2 ≡ vg# dsol DefaultGuess Options str2num : vg# vg# vg# stack : f# K# IG 1 findrows : num2str 0 ≠ IG# f# 2 col : u# IG_F# IG# eval : 12 mat vg# vg# 01020 -^IG_F# al_nleqsolve 2 el N# 1 + 1 matrix + eval : 21 line vg# 0 : if v# BF# EQ# : ε# dsol eps Options str2num : go# dsol DefaultGuess Options str2num 11 mat : t# y# y'# bf# ts# : t# num2str concat str2num u# y# y'# 12 mat : go# BF# al_nleqsolve 1 el 31 line : bf# BIC# str2num vectorize vg# M# N# ε# bvp.2 81 line Not Implemented. error if 31 line vg# dsol DefaultGuess Options str2num ky# 1 matrix + eval : f# K# IG 1 findrows : num2str 0 ≠ IG# f# 2 col : v# IG_F# IG# eval : vg# vg# 01020 -^IG_F# al_nleqsolve : 13 mat if v# vg# : dsol UpdateGuess Options 1 ≡ t# u# v# RK_D# ts# : t# num2str concat str2num v# v# EQ_F# al_nleqsolve : V# eval 31 line : t# u# RK_D# ts# : t# num2str concat str2num v# vg# EQ_F# al_nleqsolve : V# eval 31 line : if dsol dsolver Options (ic#,str2num(to#),TF#,N#,RK_D#) concat str2num 51 line if 121 line :

DE# Y# N# dsol DE# Y# 0 N# dsol 11 line :

dsol dsolver rkfixed steps 100 nsolver al_nleqsolve stepsmax 0 eps 10^{-7} DefaultGuess 1 UpdateGuess 1 MaxIters 5082 mat :

nsol nsolver al_nleqsolve stepsmax 0 eps 10^{-7} 32 mat :

Jacobi's complete elliptic integral of the first kind

m K π 2 / 12 / 12 / 1 k hyp2f1 * :

or (better approach)

m K m EllipticK evalf rational convert maple 11 line :

// INPUT
tit="Campana C1 (distesa)"
// distanza tra asse di rotazione e baricentro
// complessivo (campana+ ruota +ceppi)u.m. m
d=0.15
// massa totale(campana+ ruota +ceppi)u.m. kg
m=213
// momento d'inerzia rispetto asse rotazione u.m. kg,m
J=23.18
// metà angolo d'oscillazione %pi/2 per distesa
fimax=%pi/2+0.01

d 0.15 m * :

m 213 kg * :

J 23.18 kg * m 2^* :

φ.Max 0.5 π * 0 + :

// FINE INPUT
// lunghezza pendolo equivalente
lr=J/(m*d)
om0=sqrt(g/lr)
//
k=sin(fimax/2)
// periodo
T=4/om0*%k(k)
disp(T,"T=" )
t=linspace(0,T,500)';

lr J m d * / : 0.7255 m *

ω g.e lr / sqrt : 3.6765 Hz *

k 0.5 φ.Max * sin : 0.7071

T 4 ω / k K * : 2.0172 s *

Numeric

t u D φ φ' 12 mat u 1 el u 2 el 12 mat : φ'' eq φ'' ω s * (2^φ sin * + : φ' 0 eq al_nleqsolve stack 31 line : RK φ.Max 021 mat 04 T * 200 D Adams : 21 line

RK t φ 4 T * 200 dsol :

t φ'' ω 2^t φ sin * + 0 ≡

0 φ φ.Max ≡

0 φ' 0 ≡

dsolver Adams ≡

τ RK 1 col s * :

φ' RK 3 col Hz * :

φ RK 2 col :

φ'' ω 2^φ sin vectorize * - :

H m d * φ' 2^φ sin * φ'' φ cos * - (* vectorize :

V m g.e * m d * φ'' φ sin * φ' 2^φ cos * + (* + vectorize :

Symbolic

t φ 2 k ω t 0.25 T * + (* k sn * asin * :

Angular speed and acc

t φ' 2 k * ω * ω t 0.25 T * + (* k cn * ω t 0.25 T * + (* k dn * 1 k 2^ω t 0.25 T * + (* k sn 2^* - sqrt / :

t φ'' ω 2^t φ sin * - :

Forces

t H m d * t φ' 2^t φ sin * t φ'' t φ cos * - (* :

t V m g.e * m d * t φ'' t φ sin * t φ' 2^t φ cos * + (* + :

t s * φ τ φ augment 2 1 sys φ φ' augment 1 1 sys t s * φ' τ φ' augment 2 1 sys t s * H τ H augment 2 1 sys t s * φ'' τ φ'' augment 2 1 sys t s * V τ V augment 2 1 sys

https://en.smath.com/forum/yaf_postsm86613_Problem-in-X-Y-Plot.aspx#post86613

clear
global g d m om0 fi0 lr k T
// accelerazione di gravità u.m. m,s
g=9.806
// INPUT
tit="Campana C1 (distesa)"
// distanza tra asse di rotazione e baricentro complessivo (campana+ ruota +ceppi)u.m. m
d=0.15
// massa totale(campana+ ruota +ceppi)u.m. kg
m=213
// momento d'inerzia rispetto asse rotazione u.m. kg,m
J=23.18
// metà angolo d'oscillazione %pi/2 per distesa
fimax=%pi/2+0.01
// FINE INPUT
// lunghezza pendolo equivalente
lr=J/(m*d)
om0=sqrt(g/lr)
//
k=sin(fimax/2)
// periodo
T=4/om0*%k(k)
disp(T,"T=" )
t=linspace(0,T,500)';

deff('[u]=myplus(x)','u=2*asin(k*ellipj(om0*(x+T/4),k))' )
fi=myplus(t);
//[u]=2*asin(k*ellipj(om0*t,k))
fi1=diag(numderivative(myplus,t));
fi2=-om0^2*sin(fi);
f=1/T
omega=2*%pi/T
// forze all'asse di rotazione
H=-m*d*(fi2.*cos(fi)-fi1.^2 .*sin(fi));
V= m*g+m*d*(fi2.*sin(fi)+fi1.^2 .*cos(fi));
// scomposizione in forzanti sinusoidali per la rappresentazione in serie di Fourier
aH0=1/T*intsplin(t,real(H))
aV0=1/T*intsplin(t,real(V))
n=20
for i=1:n
omegai=i*omega;
integrando1=H .*cos(omegai*t);
integrando2=H .*sin(omegai*t);
aH(i)=2/T*intsplin(t,real(integrando1));
bH(i)=2/T*intsplin(t,real(integrando2));
integrando3=V.*cos(omegai*t);
integrando4=V.*sin(omegai*t);
aV(i)=2/T*intsplin(t,real(integrando3));
bV(i)=2/T*intsplin(t,real(integrando4));
end
aH
bH
aV
bV

Alvaro

4 appVersion 1.73.9126.0