我创建了一个质量弹簧阻尼器系统的 TikZ,该系统由以下 ODE 描述
dx1/dt = x2
dx2/dt = -w0^2*x1 -2*zeta*w0*x2 + u
y = K*w0*x1 + 2*zeta*w0*x2
其中u
是输入,y
是输出。我们假设w0=K=1
和zeta=0.1
以及初始条件x1(0)=x2(0)=0
。
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}
现在,我该如何为这段时间u=sin(a*t)
的不同时间a
点的 st 制作动画t
。y
当然,应该遵循 ODE 的解。车轮u
应该沿着地面正弦曲线“行驶”。
答案1
更新:正如 OP 和评论中所要求的那样,车轮在道路上“行驶”并移动观察者。
(点击图片运行交互式 SVG [4.2 MB]。
控制 ODE 系统已使用\pstODEsolve
包解决pst-ode
。
为了让车轮“行驶”在路面上你(X)以恒定速度,第三个微分方程,dX/天吨,需要解决。它描述了X车轮与道路接触点的坐标。微分方程组现在为
dX/天吨= ω轮 r轮子/ sqrt(1+(你'(X))2)
d是质量/d吨=五大量的
d五质量/d吨= −ω 0 2(是质量−是偏移量−是轮) − 2ζω 0 (五质量−五车轮)
积分参数吨表示时间,ω wheel表示车轮的角速度,r车轮的半径。(因此,ω车轮 r车轮 吨是沿弯曲道路轮廓行进的距离。)ω 0是无阻尼固有频率,ζ 是阻尼比。基于 dX/天吨,X(吨)、车轮半径r. 轮, 地方 道路 高程 ;你(X)、斜率你'(X)和二阶导数你''(X)轮轴坐标X车轮,是车轮及其垂直速度五车轮经过计算。是车轮和五车轮现在作为质量弹簧阻尼器的输入。
道路概况你(X)被建模为波形突发:
你(X) = 余弦(X−X抵消)埃−λ(X−X偏移)2
请注意车轮撞到“坑洼”时不太和谐的轨迹以及由此产生的响应曲线。
动画图的移动窗口是通过动态设置环境的xmin
和xmax
选项来实现的axis
。
ODE 求解器(RKF45pkg 中的算法)pst-ode
是用 PostScript 编码的。最近,Marcel Krüger 实现了一个PostScript 解释器在 Lua 中,ps2pdf
不再需要这样做。
运行lualatex
3 次即可排版为 PDF。对于 GIF 和 SVG 输出,请取消注释源代码开头的其中一行,然后按照其中的说明进行操作。
%\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}