$$
\def\DIS{\displaystyle}
\def\fraction(#1/#2){\frac{#1}{#2}}
\def\ff(#1/#2){\frac{#1}{#2}}
\def\dd#1#2#3{\frac{d ^{#1}#2}{d #3^{#1}}}
\def\pdd#1#2{\frac{\partial ^{#1}#2}{\partial #3^{#1}}}
\def\bi#1#2{{}{#1}! , \mathrm C{#2}}
\def\bm{\boldsymbol}
$$

調和振動子の方程式

ばね定数 $k$, 自然長 $\ell$ のバネにつながった質点の運動方程式は

$$ m \dd{2}{x}{t} = -k(x-\ell) $$

$m$ は質点の質量。

この方程式の右辺は次のように書き換えます。

$$ - \ff(k/m) \ff((x-x_0)/|x-x_0|) (|x-x_0| - \ell) $$

$x_0$はばねのもう一方の端点の座標です。このように書き換えておくと，ばねの右側か左側かによらず同じ式が使えます。この関数はTeXでは次のように書きます。

\def\funcAccelPX#1#2#3#4#5{
    (-(#1-#3)/abs(#1-#3)*\paraK*(abs(#1-#3) - \Nlength)/(#5))}%

この調和振動子の微分方程式をTeXで数値的に解きます。

https://youtu.be/ErUNg6D16js

積分の手法はこちら
https://qiita.com/hironarin/items/cc4e1c8f961155a8e511

TeX

tikzpictureの中身以外は以下のようになります。

%計算回数
\def\stepB{60}%1秒の刻み幅
\def\stepA{20}%時間
%カウンター定義
\newcount\countA%
\newcount\countB%

%%%%%%%%%%%%%%%%%%%%%%%%%
%物理パラメータ 定数
%%%%%%%%%%%%%%%%%%%%%%%%%
        \tikzmath{
        int \ParNum ;
        \ParNum=1; %粒子数
        \Nlength = 6;%ばねの自然
        \paraK=3;%ばね定数
        \ceiling = 8;%ばねの端点の座標（固定）
        \mass{1} = 0.5;%質点の質量
		}

\tikzmath{
    %時間刻み幅の設定
    real \dT;
    \dT = 1 / \stepB;
}


%微分方程式の右辺
%#1 x座標	 #2 速度x成分 #3ばねのもう片方の端点の座標 #4ばねのもう片方の端点の速度 #5質点の質量
\def\funcAccelX#1#2{#2}%dH/dp_x
\def\funcAccelPX#1#2#3#4#5{
    (
    -(#1-#3)/abs(#1-#3)*\paraK*(abs(#1-#3) - \Nlength)/(#5))}%-dH/dx
%%%%%%%%%%%%%%%%%%%%%%%%%





%ぬこの初期値を代入
\tikzmath{
    \iqXNuko{1} = 4;
    \ivXNuko{1} = 0;
}
%グローバル変数に代入する
\countA = 0 \loop\ifnum\countA < \ParNum{
    \tikzmath{
        int \intCount;
        \intCount = \countA + 1;
    }
    \global\expandafter\edef\csname variXCdNuko[\intCount] \endcsname{\iqXNuko{\intCount}}%初期位置を代入
    \global\expandafter\edef\csname variXVelocityNuko[\intCount] \endcsname{\ivXNuko{\intCount}}%初速度を代入
}
\advance\countA by1\repeat %







\countA = 0 \loop\ifnum\countA < \stepA
{
\countB = 0 \loop\ifnum\countB < \stepB
{



\begin{tikzpicture}

%%%%ここで描画

\end{tikzpicture}





%%%数値計算ルンゲクッタ
\tikzmath{
    \NqtXTako = \csname variXCdNuko[1] \endcsname;
    \NvtXTako = \csname variXVelocityNuko[1] \endcsname;
}


\tikzmath{
    real \variFuji;
    \dummy{1}=\NqtXTako;
    \dummy{2}=\NvtXTako;
    \NatXTako = \funcAccelX{\dummy{1}}{\dummy{2}};
    \NatPXTako =  \funcAccelPX{\dummy{1}}{\dummy{2}}{\ceiling}{0}{\mass{1}};
%
%
    \dummy{1}=\NqtXTako + 0.5 * \NatXTako *\dT;
    \dummy{2}=\NvtXTako + 0.5 * \NatPXTako * \dT;
    \NatXTako = \funcAccelX{\dummy{1}}{\dummy{2}};
    \NatPXTako =  \funcAccelPX{\dummy{1}}{\dummy{2}}{\ceiling}{0}{\mass{1}};
%
%
    \NvtXTako = \NvtXTako +  \dT  * \NatPXTako;
    \NqtXTako = \NqtXTako + \NatXTako * \dT;
    %
}
%計算結果を代入
\global\expandafter\edef\csname variXCdNuko[1] \endcsname{\NqtXTako}
\global\expandafter\edef\csname variXVelocityNuko[1] \endcsname{\NvtXTako}


%計算終了，次のページへ
\newpage 
}%loop \countB
\advance\countB by1\repeat

}%loop \countA
\advance\countA by1\repeat