Eu criei um TikZ de um sistema massa-mola-amortecedor que é descrito pelas seguintes EDOs
dx1/dt = x2
dx2/dt = -w0^2*x1 -2*zeta*w0*x2 + u
y = K*w0*x1 + 2*zeta*w0*x2
onde ué a entrada e yé a saída. Vamos supor w0=K=1ee zeta=0.1as condições iniciais x1(0)=x2(0)=0.
O TikZ é
\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{arrows, shapes, patterns, calc, decorations.pathmorphing, decorations.markings, positioning, animations}
\begin{document}
\begin{tikzpicture}
\tikzstyle{wheel} = [draw, circle, minimum size = 0.75cm]
\tikzstyle{mass} = [draw, rectangle, minimum height = 1cm, minimum width = 2cm]
\tikzstyle{spring} = [decorate, decoration = {zigzag, pre length = 0.3cm, post length = 0.3cm, segment length = 6}]
\tikzstyle{damper} = [decoration = {markings, mark connection node = dmp, mark = at position 0.5 with
{
\node (dmp) [inner sep = 0pt, transform shape, rotate = -90, minimum width = 5pt, minimum height = 10pt, draw=none] {};
\draw ($(dmp.north east)+(1.5pt,0)$) -- (dmp.south east) -- (dmp.south west) -- ($(dmp.north west)+(1.5pt,0)$);
\draw ($(dmp.north)+(0,-1.5pt)$) -- ($(dmp.north)+(0,1.5pt)$);
}
}, decorate]
\tikzstyle{groundflat} = [fill, pattern = north east lines, draw = none, minimum width = 0.75cm, minimum height = 0.3cm]
\draw[fill, pattern = north east lines] (-.5,-.25) cos (0,0) sin (.5,.25) coordinate (top) cos (1,0) sin (1.5,-.25) cos (2,0) sin (2.5,.25) cos (3,0) sin (3.5,-.25) |- (-.5,-1)--cycle;
\node[wheel, above = 0pt of top] (u) {$u$};
\node[mass, above = 2cm of u] (m) {$y$};
\draw [-] (u.north) |- ++(0.5,0.25cm)coordinate (uright);
\draw [-] (u.north) |- ++(-0.5,0.25cm)coordinate (uleft);
\draw [spring] (uleft) -- ($(m.south) +(-0.5,0)$);
\draw [damper] (uright) -- ($(m.south) +(0.5,0)$);
\end{tikzpicture}
\end{document}
Agora, como posso animar esse st u=sin(a*t)com diferentes aao longo de um determinado tempo t. ydeveria, é claro, seguir a solução das EDOs. E a uroda deve “andar” ao longo da sinusóide do solo.
Responder1
Atualizar:Roda “andando” na estrada e observador em movimento, conforme solicitado no OP e comentários.
(Clique na imagem para executar o SVG interativo [4,2 MiB] no navegador.)
O sistema ODE governante foi resolvido com \pstODEsolveo pacotepst-ode.
Para a roda “andar” no perfil da estradavocê(x) em velocidade constante, uma terceira equação diferencial, dx/dt, precisa ser resolvido. Descreve o movimento doxcoordenada do ponto de contato da roda com a estrada. O sistema de equações diferenciais é agora
dx/dt= ω roda Rroda / sqrt(1+(você'(x)) 2 )
dsimmassa /dt=vmassa
dvmassa /dt= −ω 0 2 (simmassa -simdeslocamento -simroda ) − 2ζω 0 (vmassa -vroda )
O parâmetro de integraçãotrepresenta o tempo, ω roda a velocidade angular da roda eRroda seu raio. (Portanto, roda ω Rroda té a distância percorrida ao longo do perfil curvilíneo da estrada.) ω 0 é a frequência natural não amortecida e ζ é a taxa de amortecimento. Baseado em dx/dt,x(t), raio da rodaRroda , elevação da estrada localvocê(x), declivevocê'(x) e 2ª derivadavocê''(x) as coordenadas do eixo da rodaxroda ,simroda e seu vert. velocidadevroda são calculados.simroda evroda agora serve como entrada para o sistema massa-mola-amortecedor.
O perfil da estradavocê(x) é modelado como uma explosão de forma de onda:
você(x) = cos(x-xdesvio )e-λ(x-xdeslocamento ) 2
Observe a trajetória não tão harmônica da roda onde ela bate nos “buracos”, bem como a curva de resposta resultante.
A janela móvel do enredo animado é obtida definindo dinamicamente as opções xmine xmaxdo axisambiente.
O solucionador EDO (RKF45método) em pkg pst-odefoi codificado em PostScript. Recentemente, Marcel Krüger implementou umIntérprete PostScriptem Lua tal que ps2pdfnão é mais necessário.
Composto como PDF executando lualatex3 vezes. Para saída GIF e SVG, remova o comentário de uma das linhas no início da fonte e siga as instruções nela contidas.
%\PassOptionsToPackage{export}{animate} % lualatex <file> ; convert -density 300 -alpha remove -delay 4 <file>.pdf <file>.gif
%\PassOptionsToPackage{dvisvgm,autoplay}{animate} % dvilualatex <file> ; dvisvgm --zoom=-1 --font-format=woff2 --bbox=papersize <file>.dvi
\documentclass[margin=3pt]{standalone}
\usepackage{pst-ode}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\usepgfplotslibrary{fillbetween}
\usetikzlibrary{arrows.meta, calc, decorations}
\usepackage{listofitems} % read space separated items
\usepackage{animate}
\usepackage{xsavebox} %\xsbox
% adjustable parameters & definitions
\def\tEnd{80} % max. time (integration parameter)
\def\animationframes{201} % number of output points for animated objects
\def\outputpoints{1001} % number of output points for curves
\def\omegaWheel{1.0} % wheel angular velocity
\def\rWheel{1.0} % wheel radius
\def\rCenterOfMass{0.16em} % center-of-mass symbol size
\def\springlength{6} % length of spring in neutral state (and of damper cylinder)
\def\springturns{10} % (integer!) number of spring windings
\def\windowsizeright{5} % moving window x-size around moving object
\def\windowsizeleft{25}
\def\ymin{-1} % y axis limits
\def\ymax{17}
\pstVerb{
tx@Dict begin
% mass-spring-damper properties
/w0 1.0 def
/zeta 0.1 def
% initial conditions
/x0 0.0 def % x(t0)=0 starting x coordinate
/yMass0 0.0 def % y coord of mass centre (/wo vert. offset)
/vMass0 0.0 def % y velocity of mass centre
% road profile U(x) with wave-form burst, and its first and second derivatives U'(x), U''(x)
/A 1 def % amplitude
/a 1 def % mode
/lambda (1/40) AlgParser cvx exec def % burst growth/decay rate
/xOff 20.0 def % burst offset
/U (A*Euler^(-lambda*(x[0]-xOff)^2)*cos(a*(x[0]-xOff))) AlgParser cvx bind def
/dUdx (-A*Euler^(-lambda*(x[0]-xOff)^2)*a*sin(a*(x[0]-xOff))-2*lambda*(x[0]-xOff)*U) AlgParser cvx bind def
/d2Udx2 (-(a^2+2*lambda+4*lambda^2*(x[0]-xOff)^2)*U-4*lambda*(x[0]-xOff)*dUdx) AlgParser cvx bind def
% do not change anything below
/rWheel \rWheel\space def
/omegaWheel \omegaWheel\space def
/massWheelOff (rWheel+1.5*\springlength+2) AlgParser cvx exec def
/dxdt (rWheel*omegaWheel/sqrt(1+dUdx^2)) AlgParser cvx bind def
% local wheel hub coordinates
/xWheel (x[0]-rWheel*dUdx/sqrt(1+dUdx^2)) AlgParser cvx bind def
/yWheel (U+rWheel/sqrt(1+dUdx^2)) AlgParser cvx bind def
/dyWheeldx (dUdx*(1-d2Udx2/(1+dUdx^2)*(yWheel-U))) AlgParser cvx bind def
% angular position of wheel in terms of t
/thetaWheel {omegaWheel t mul 2 Pi mul div dup cvi sub neg 360 mul} bind def
end
}
% solve equations of motion with many output points for smooth curves
\pstODEsolve[algebraicAll,saveData]{curves}{% PS variable and file that take result table
% table format to be saved in `curves' and `curves.dat':
x[0] | % x coordinate of contact point
U | % y coordinate U(x) of contact point
xWheel | yWheel | % wheel hub x (same as that of mass) and y coordinates
x[1]+massWheelOff+rWheel % y coordinate of mass (response with const. vertical offset)
}{0}{\tEnd}{\outputpoints}{% t_0, t_end, number of t steps plus 1
x0 | yMass0 | vMass0 % initial conditions
}{ % ODE system's RHS:
dxdt | % dx/dt, movement of contact point (x coordinate)
x[2] | % centre of mass vert. velocity
-w0^2*(x[1]-yWheel+rWheel)% centre of mass vert. acceleration
-2*zeta*w0*(x[2]-dyWheeldx*dxdt)
}
% solve equations with fewer output points for mass-spring-damper animation; save solution
% table (wheel and mass centre coordinates, valve angular position) in `objects' and `objects.dat'
\pstODEsolve[algebraicAll,saveData]{objects}{
xWheel | yWheel | x[1]+massWheelOff+rWheel | thetaWheel-90
}{0}{\tEnd}{\animationframes}{x0 | yMass0 | vMass0}{
dxdt | x[2] | -w0^2*(x[1]-yWheel+rWheel)-2*zeta*w0*(x[2]-dyWheeldx*dxdt)
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \fileopenr{<file stream>}{<file name>}, opens file for reading
\newcommand\fileopenr[2]{%
\newread#1%
\immediate\openin#1=#2%
}
% \readtolist[<sep char>]{<file stream>}{\list}
% reads a line from file stream and splits at <sep char> into \list[1], \list[2], ...
\newcommand\readtolist[3][,]{{%
\setsepchar{#1}%
\immediate\read#2 to \inputline%
\ifeof#2
\immediate\closein#2%
\ifdefined\multiframebreak\multiframebreak\fi%
\else%
\greadlist*#3\inputline%
\fi%
}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\tikzset{spring/.style 2 args= {decorate, decoration = {zigzag, segment length = \dimexpr(#1-1pt)/\springturns\relax, amplitude=#2}}}%
\begin{document}
\IfFileExists{curves.dat}{}{dummy text\end{document}}%
\xsbox{centreofmass}{%
\tikz[radius=\rCenterOfMass] {% from https://tex.stackexchange.com/questions/23453
\fill (0,0) -- ++(\rCenterOfMass,0) arc [start angle=0,end angle=90] -- ++(0,{-2*\rCenterOfMass}) arc [start angle=270, end angle=180];%
\draw [thin] (0,0) circle;%
}%
}%
\begin{animateinline}[controls]{25}
\fileopenr{\data}{objects.dat}
\readtolist[ ]{\data}{\table}
\multiframe{10000}{}{
\begin{tikzpicture}[line cap=rect, line join=bevel]
\begin{axis}[
unit vector ratio = 1, axis on top,
xmin={\table[1]-\windowsizeleft}, xmax={\table[1]+\windowsizeright},
ymin=\ymin, ymax=\ymax,
xtick distance=5,
ytick distance=2,
xlabel={$x$}, ylabel={$y$}, ylabel style={rotate=-90},
extra x tick labels={\phantom{\strut00}},extra x ticks={\pgfkeysvalueof{/pgfplots/xmax}}, % prevents wobbling of axis
]
% road profile; off-window coordinates filtered away for smaller file size
\addplot [name path=road,no markers,unbounded coords=jump,x filter/.expression={x<\pgfkeysvalueof{/pgfplots/xmin}-0.1 || x>\pgfkeysvalueof{/pgfplots/xmax}+0.1 ? nan : x}]
table [x index=0, y index=1] {curves.dat} -- (\pgfkeysvalueof{/pgfplots/xmax},0);
\path [name path=bottom] (0,\pgfkeysvalueof{/pgfplots/ymin}) -- (\pgfkeysvalueof{/pgfplots/xmax},\pgfkeysvalueof{/pgfplots/ymin});
\addplot [fill=black!20!white] fill between[of=road and bottom];
% wheel hub trajectory
\addplot [no markers,blue,unbounded coords=jump,x filter/.expression={x<\pgfkeysvalueof{/pgfplots/xmin}-0.1 || x>\table[1] ? nan : x}]
table [x index=2, y index=3] {curves.dat};
% centre of mass trajectory
\addplot [no markers,red,unbounded coords=jump,x filter/.expression={x<\pgfkeysvalueof{/pgfplots/xmin}-0.1 || x>\table[1] ? nan : x}]
table [x index=2, y index=4] {curves.dat};
\coordinate (w) at (\table[1],\table[2]); % wheel hub
\coordinate (m) at (\table[1],\table[3]); % centre of mass
% wheel with valve
\fill[even odd rule] (w) circle [radius={0.7*\rWheel}] (w) circle [radius=\rWheel];
\draw[thick] ($(w)+(axis direction cs:{0.6*\rWheel*cos(\table[4])},{0.6*\rWheel*sin(\table[4])})$)
-- ++(axis direction cs:{0.1*\rWheel*cos(\table[4])},{0.1*\rWheel*sin(\table[4])});
% mass
\node at (m) {\thecentreofmass};
\draw (m) ++(axis direction cs:-1,-1) rectangle ++(axis direction cs:2,2);
% spring/damper support
\draw (m) ++(axis direction cs:-1,-1) -- ++(axis direction cs:0,{-0.25-0.25*\springlength}) coordinate (sprtop);
\path (m) ++(axis direction cs:1,-1) coordinate (dmptop);
\draw (w) ++(axis direction cs:0,{\rWheel+0.5}) -| ++(axis direction cs:-1,{0.25+0.25*\springlength}) coordinate (sprbot);
\draw let \p1=(axis direction cs:0.35,0) in [{Circle[open, length=\x1]}-, shorten <=-0.5*\x1] (w) -- ++(axis direction cs:0,{\rWheel+0.5}) -| ++(axis direction cs:1,0.5) coordinate (dmpbot);
% spring
\path let \p1=($(sprtop)-(sprbot)$), \p2=(axis direction cs:0.6,0) % length and amplitude (half width)
in [draw,spring={\y1}{\x2}] (sprbot) -- (sprtop);
% damper lower part
\draw (dmpbot) -| +(axis direction cs:0.4,\springlength) (dmpbot) -| +(axis direction cs:-0.4,\springlength);
% damper upper part
\draw (dmptop) -- ++(axis direction cs:0,-\springlength) ++(axis direction cs:-0.3,0) -- ++(axis direction cs:0.6,0);
\end{axis}
\end{tikzpicture}%
\readtolist[ ]{\data}{\table}
}
\end{animateinline}
\end{document}