
我創建了質量彈簧阻尼器系統的 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)
以不同的方式動畫這個 st 。當然,應該遵循 ODE 的解決方案。車輪應該沿著地面正弦曲線「行駛」。a
t
y
u
答案1
更新:按照OP和評論中的要求,車輪在道路上「騎行」並移動觀察者。
(點擊圖片即可運行互動式 SVG [4.2 米巴] 在瀏覽器中。
控制 ODE 系統已\pstODEsolve
透過套件解決pst-ode
。
讓車輪在路面上“行駛”你(X)以恆定速度,第三個微分方程 dX/dt,需要解決。它描述了X車輪與路面接觸點的座標。現在的微分方程組是
dX/dt= ω輪 r輪/ sqrt(1+(你'(X)) 2 )
dy品質/天t=v大量的
dv品質/天t= - ω 0 2 (y品質-y偏移量-y輪) − 2ζω 0 (v品質-v車輪)
積分參數t代表時間,ω輪子的角速度和r輪子的半徑。 (因此,ω輪 r車輪 t是沿著彎曲道路輪廓行駛的距離。基於dX/dt,X(t), 車輪半徑r. 車輪, 當地 道路 高程 ;你(X), 坡度你'(X) 和二階導數你''(X) 輪軸座標X車輪,y輪子及其垂直度。速度v輪進行計算。y輪子和v車輪現在用作質量彈簧阻尼系統的輸入。
道路概況你(X) 被建模為波形突發:
你(X) = 餘弦(X-X抵消)e-λ(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}