弦を考えます。弦は重さ\(M\), 長さ\(L\), ばね定数\(K\)を持ちます。

この弦を波動方程式に頼らないで定式化しましょう。
波動方程式は、2次元の波動のうち、1次元分の方向を\(y=y(x, t)\)として表現するため、x軸にとって多価関数となる弦を表現することが出来ません。物性物理学等で格子の振動を考える場合は、結晶を構成する原子が平衡位置から殆ど動かないのでこういった問題は起こりません（起こったとしても無視される位少ないので、大多数の振る舞いを考える場合には問題になりません）。
よって多次元の問題を考えるには波動方程式を元にして考えるのは不適切です。そのため、弦を単なる連成振動として考えて、完全に2次元の問題として扱いたいと思います。

図のような連成振動のモデルを考えます。

質量\(m\)を持つ\(i\)番目の質点が、位置ベクトル\(\mathbf{x}_i\)で指定され、隣り合う質点との間が、ばね定数\(k\)、自然長\(l\)を持つ理想的なばねで繋がれているとします\((i=1,2,\cdots, N)\)。この時、運動方程式を立てると
\begin{eqnarray}
m\frac{d^2 \mathbf{x}_i}{dt^2}&=&
+k(|\mathbf{x}_{i+1}-\mathbf{x}_{i}|-l) \frac{\mathbf{x}_{i+1}-\mathbf{x}_{i}}{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|}
-k(|\mathbf{x}_{i}-\mathbf{x}_{i-1}|-l)\frac{\mathbf{x}_i-\mathbf{x}_{i-1}}{|\mathbf{x}_i-\mathbf{x}_{i-1}|}
\end{eqnarray}
となります。全く同じですが、見通しをよくするため式変形すれば、
\begin{eqnarray}
m\frac{d^2 \mathbf{x}_i}{dt^2}&=& k\Bigl[A_i \mathbf{x}_{i+1}-\bigl(A_i+A_{i-1}\bigr)\mathbf{x}_{i}+A_{i-1} \mathbf{x}_{i-1}\Bigr]
\end{eqnarray}
と表せます。ここで、
\begin{equation}
A_i = \frac{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|-l}{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|},~~(i=1,\cdots, N, A_N = 0)
\end{equation}
と置きました。

さて、こうしたモデルを考えた時に、連成振動として考えた\(N\)分割した1つ1つの小さいばね（”ミクロなばね”と名づけます）とその弦を1つにみなした時の大きなばね（”マクロなばね”と名づけます）の間の関係式を考えましょう。
仮定として、等間隔に分割し、弦の密度は一様であるとします。

この条件の下で、
マクロなばねの重さ\(M\), 自然長\(L\), ばね定数\(K\)
と
ミクロなばねの重さ\(m\), 自然長\(l\), ばね定数\(k\)
の間の関係式を導出します。

結論として、マクロとミクロなばねの間に成り立つ関係は
\begin{eqnarray}
L&=&N l \\
M&=&N m \\
K &=& k/N
\end{eqnarray}
となります。

自然長

自然長\(L\)の弦を\(N\)分割するので、
\begin{equation}
L=N l
\end{equation}
の関係があります。

重さ

重さ\(M\)の弦を\(N\)分割するので、
\begin{equation}
M=N m
\end{equation}
と書けます。

ばね定数

最後にばね定数の関係を考えましょう。この関係は、弾性エネルギーを考えることで導くことができます。
マクロなばねとして全体を見た場合に長さ\(D\)だけ縮んだ時に蓄えられる弾性エネルギーと
同じだけ縮んだ時にミクロな\(N\)個のばねに蓄えられる弾性エネルギーが同じであってほしいと要請します。すなわち、ミクロなばねの縮みの合計がマクロなばねの縮み\(D\)に等しいとすると、
\begin{equation}
D=N d
\end{equation}
と書けます。
実は、大胆な過程をしないとマクロなばねとミクロなばねの関係を求めることが出来ません。
なぜなら、マクロなばねはいつでも線形応答をするかのように現在考えていますが、ミクロなばねでは場所ごとに違っても良いという様に定式化しています。なのでもともと無理な話なのです。
ではどうやって妥当な関係式を導出すればよいかといえば、ばねが非常にゆっくり運動が行われるときを仮定するのです。ゆっくり動くときの振る舞いは全てのばねが同じ動きをするとして差し支えないでしょう。そうして定式化します。
よって、
\begin{eqnarray}
\frac{1}{2}K D^2 &=& N\cdot \frac{1}{2} k d^2 \nonumber \\
\to ~~ K &=& k/N
\end{eqnarray}
という関係式が導けます。

あとは数値計算を行うだけです。陽的ルンゲクッタ法で解いたプログラムを以下に示します。

https://slpr.sakura.ne.jp/qp/supplement_data/wave_particles.tar.gz

実行は、

でokだと思います。初期状態や、重さや自然長はinputファイルを見てください。

実行後、動画をgnuplot上で出力したければ、

としてください。100は表示する時間ステップを表します。この引数の値の最大はinputファイルの中のNtまでです。

数値計算結果

実際に数値計算を行う条件として、自然長や重力があると波動方程式と若干異なる振る舞いになってしまうので、対応を見る上ではなくしています。すなわち、自然長\(L=0\)、重力\(g=0\)として計算します。

三角パルス

三角形の波がパルス状に存在するとします。
そして初期状態として、初速度がゼロであることを仮定しましょう。

この場合、波動方程式の考察から、
同じ関数であらわされる右向きの波\(f(x-vt)\)と左向きの波\(f(x+vt)\)の重ね合わせとして書かれることが分かっています。すなわち、解\(y=y(x,t)\)は
\begin{equation}
y(x,t)=f(x-vt)+f(x+vt)
\end{equation}
と書かれます。図で表せば以下のように時間発展していくことが期待されます。

これが本稿の連成振動モデルで実際にそうなっているのかを確かめるのは良い問いかけでしょう。
実際に連成振動を数値的に解くと以下の通りになります。

ゆっくりと時間発展を見ていくと

となります。明らかに右向きに進む波と左向きに進む波として分離します。分離した結果から、初期状態がそれらの波の重ね合わせで表現されていたことが分かるでしょう。

矩形パルス

矩形パルスを考えてみます。実際に数値計算をするとこんな感じです。

結局右向きと左向きの重ね合わせとしてあらわされることも確認できますが、決して単純なそれらの足し合わせで書かれているわけではないことがわかります。この振る舞いは、波動方程式に対応するように連成振動モデルを作っていないことから由来します。波動方程式としてあらわした二階微分の偏微分方程式の場合は多価関数をとることができません。なので縦に弦が伸びる場合、間に質点を置くことができないので質点の位置がいきなり飛ぶことになります。しかし、連成振動モデルではその制約をして初期状態を準備しているわけではありませんのでその違いが出て来ます。図で表せば、以下のようになります。

両者の違いがあるので、波動方程式の結論である
\begin{equation}
y(x,t)=f(x-vt)+f(x+vt)
\end{equation}
が成り立っていない、と解釈することができます。

少しずつ時間発展をみるとこの通りになります。

三角形の波

端まで三角形のパルスであった場合にどのように振る舞うかを見てみましょう。特に面白い振る舞いがあるわけではなく、上のパルスの、右向きと左向きの波の重ね合わせで説明ができるので、単にシミュレーション結果を載せるだけにとどめます。実際に計算してみると以下のようになります。

星形

初期状態が星型である場合も見てみましょう。振る舞いも同じです。

注意点としては、実際にこの振る舞いは起こりえないことに注意しておきましょう。なぜなら、本当の弦で行った場合、重なる場所があると弦の衝突が起こります。しかし、この連成振動モデルでは計算上、衝突は起こりません。

偏りがある三角形のパルス

偏りがある場合も波動方程式では記述することができません。こんな感じの初期状態です。

これを実際にシミュレートすると以下の通りになります。

3次元の場合

3次元の場合も載せておきます。ふるまい自体はそんなに変わるわけではないので、シミュレーション結果を載せるだけにとどめておきます。

動かない壁に対する束縛運動と反射を考えます。
例えば、初め跳ねてた運動が、壁に沿って動く運動に変化する、という状況です。
あんまり見たことが無いので、面白そうだと思いました。

束縛されている時と反発するときは、式(1),(2)によってあらわされます。それらは
壁に束縛されている場合

壁と反発する場合

です。壁と反発する場合、反発後の速度は式(3)に沿って動きます。

※式(1), (2)では壁が時間依存しても良いように定式化しています。この定式化は恐らく正しいです。また、本稿に載せているプログラムも壁が時間依存しても良いように作成していますが、動く壁の場合、プログラムではうまく計算が出来ません。

ここで、\(e\)は反発係数、\(\mathbf{n}\)は壁の法線ベクトルであり、

\(\nabla\)はナブラ演算子であり

と与えられます。また、\(\hat{H}f\)は関数\(f\)のヘッセ行列であり,

と与えられます。

定式化や数値計算手法の詳細は以下のページを参照してください。

壁との反発と束縛運動の定式化

質点と壁との反発を表す運動方程式
束縛条件下の運動 – 自由度がうまく落とせない運動

数値計算手法

ルンゲ=クッタ法の説明と刻み幅制御
Hyper-dual numbersによる二階偏微分の計算
ゼロ点を探す（二分法、挟み撃ち法、Anderson-Björk法、Brent法、Newton法、Steffensen法）

解放・束縛判定

ここでいう”解放”とは、束縛されていなく、反射を繰り返している状態を表します。また、”束縛”は壁に沿って運動している状態です。

前提として、壁を通り抜けることは無いと考えます。

すなわち、時刻\(t=t_0\)で位置\(\mathbf{r}=\mathbf{r}_0\)の時、任意の時刻\(t\)について

が満たされるとします。
もし、\(f(\mathbf{r}_0,t_0) =0\)ならば判断がつかないため、プラスかマイナスはこちらから与えます。

解放→解放判定、解法→束縛判定

解放状態から壁によって単に反発する場合

関数\(f(\mathbf{r},t)\)の符号が変化した時、壁との反発を考えます。
壁の法線方向の速度が十分に大きい場合、壁と反発し、そうでない場合、壁に束縛されると考えます。
その前提の元、

を満たす\(t=t_i, \mathbf{r}=\mathbf{r}_i\)を見つけます。
その後、式(6)に従い、速度を変化させます。

数値計算的には、衝突直前の時刻を採用します。すなわち、
ゼロ点を探す際にある範囲\(t_a \leq t_i \leq t_b\)で挟み込んでいくのですが、\(t_b\)は壁を越えてしまうので採用しません。

解放状態から壁沿いに束縛される場合

もし、壁の法線方向の速度が十分に小さい場合（ある閾値を下回った場合）、壁に束縛されると考えます。
この時、壁の法線方向の速度はゼロに変更します。
すなわち、速度は時刻\(t=t_i\)で

を持ちますが、\(v_{\parallel}=0\)にしてから、束縛運動に移行するということです。
これは、\(\mathbf{v}\)と壁の法線方向のベクトル\(\mathbf{e}_n\)の内積を取ることで得られます。
また、束縛された瞬間(\(t=t_i\))の束縛力\(C(t_i)\)を計算し、その符号を記録しておきます。
束縛力\(C(t)\)は、

です。
この符号が変化した瞬間が壁からの束縛が無くなる時なので、そのために記録します。

束縛→解放判定

束縛力\(C(t)\)の符号と\(C(t_i)\)の符号が変わるまで、式(1)に従い、時間発展させます。
すなわち、

を満たす\(t=t’, \mathbf{r}=\mathbf{r}’\)を見つけます。
その後、式(2)に従い運動します。
式(2)の運動では束縛力は働かないので、符号は自然と初期条件の符号と同一になる（※）。

※この条件はあまり良くありません。この判別方法のせいで、壁が時間依存している場合、束縛力が働いていない一瞬で質点が壁を超えてしまいます。プログラム自体は壁は時間依存しても良いことになっていますが、この条件分岐は上手く動きません。以下に示す計算プログラムは、質点が束縛されている場合に壁が時間依存しなければ正しいです。

プログラム

プログラムは以下のリンク先においておきます。

https://slpr.sakura.ne.jp/qp/supplement_data/lag_ver1.0.tar.gz

適当なディレクトリで展開し、lag_ver1.0というディレクトリに移動してから以下のコマンドで実行できます。

動かすと以下のような動画が描画されます。

デフォルトでは
ポテンシャル\(V=mgy\)（サブルーチンfp2d）
壁\(f(x,y)=x^2+y^2-4\)（サブルーチンfw2d）
に書かれています。
摩擦、空気抵抗は入っていません。唯一、反発係数(els)だけがinputファイルの中に書かれています。

初期条件の

だけを上の通り変更すると以下のような振る舞いをします。

確かめ

確かめを行います。
重力\(g\)の下で、質量\(m\)の物体が、半径\(r\)の円形の壁の内側に沿ってボールが進み、壁からの抗力が無くなり、壁から離れる状況を考えます(下の図を参照)。

円形の壁に沿っている時、垂直抗力\(N\)は

と書けます。エネルギー保存則より、

が成り立っています。ここで、\(v_0\equiv v(t=0)\)と置きました。
垂直抗力\(N\)がゼロになる点が壁から離れる条件ですので、\(v_0\)を用いて

と書けます。初速度が分かっている時、壁から離れる角度は

です。もしくは、壁から離れる角度が分かっている場合、初速度は

と与えられます。重力加速度, 半径を\(g=1, r=2\)とし、\(y=1\)と壁との交点、すなわち\(\theta=2\pi/3\)の場合、初速度は

です。実際に本稿のプログラムを動かしますと、\(y=1\)でちょうど離れていることが確認できます。

ここで、青線は壁に沿って動いて運動している状態であり、赤線は壁から離れている運動している状態です。

剛体に対する散乱

ポテンシャルを無くし(g=0)、弾性散乱(els=1)を考えると
古典的な散乱問題的なものもできます。

質点が壁に衝突し、反発することを数式で表現します。

1. 壁の定義
2. 運動方程式の導出
   • 1次元の壁との衝突を表す運動方程式
   • 多次元の壁との衝突を表す運動方程式
   • 多次元の反発
3. プログラム
4. 実行結果
   • 重力無しの運動、動かない壁
   • 重力有りの運動、動かない壁
   • 重力有りの運動、動く壁
5. 参考文献

壁の定義

壁とは、壁の法線方向に平行な質点の速度成分を反転させるデバイスである。
と定義します。

運動方程式の導出

まず議論を簡単にする為に一次元の運動について考え、その後多次元の運動について定式化をしていきます。

1次元の壁との衝突を表す運動方程式

ポテンシャル\(V(x)\)保存力下の質点の運動を考えます。
壁が
\(f(x,t)=0\)
で表現されているとします。
壁との衝突では位置は連続、速度は不連続な振る舞いをすると考えると、壁に衝突した時、力はデルタ関数の振る舞いをしていることが予想できます。
よって、未知の定数\(c\)を用いて、

と書くことが出来ます。
ここで、時刻\(t＝t_i\)は、
\(
f(x(t),t)=0~~\to~~t=t_i
\)
を満たす時刻(壁との衝突時刻)です。\(f(x,t)\)の時間微分は、

と書くことが出来ます。
\(c\)を定めるために式(2)の両辺を、時刻\(t=t_j\)周りを微小時間\(\Delta\)で積分します。すると

を得ます。ここで、保存力の時間変化は連続であることを仮定します。すなわち、

がいかなる時刻\(t\)で成立するとします。計算を進めると、

という結果が導かれます。ここで、\(\delta_{i,j}\)はクロネッカーのデルタを表します。式(6)を式(2)の右辺に代入し、\(c\)を消去すると、運動方程式

を得ます。式(7)はこのままでは解くことが出来ません。なぜなら、未来の時刻の速度\(v(t_i+0)\)が含まれているからです。これをどうにかするには、新たな条件式、すなわち反発前後の条件式が必要になります。

壁と衝突する場合、質点の衝突前後の速度\(v(t_i-0), v(t_i+0)\)と、壁の衝突前後の速度\(v_\text{w}(t_i-0), v_\text{w}(t_i+0)\)の間には、反発係数\(e\)を用いて

の関係があると予想します。今、壁の質量が無限大であり、衝突前後で速度変化がない場合を考えると、壁の速度は、質点との衝突によって変化しないと考えます。すなわち

が成立していると考えます。すると、衝突後の質点の速度は、\(e\)を用いて

と書けます。よって、式(7)に代入して、

を得ます。

最後に、位置\(x(t_i)\)における壁の速度\(v_\text{w}(t_i)\)を考えましょう。衝突の前後のごく短い時間では、質点の動きは壁に追従すると考えます。すなわち、

が成り立っているとします。衝突時には質点の位置\(x\)と衝突位置は同じであることを注記します。式(12)を書き換えると、衝突が起こる時刻\(t=t_i\)の周りで

が成立しています。このことから、衝突前後のごく短い時間では

が成立します。よって、

から、壁の速度

を得ます。よって、壁との衝突を記述する運動方程式は

となります。

多次元の壁との衝突を表す運動方程式

壁

と質点が衝突することを考えます。運動方程式は

であり、1次元の場合と同様に時刻\(t=t_j\)周りの微小時間で積分して

を得ます。式(20)に代入すれば運動方程式

を得ます。次の節で述べる結果(多次元の反発)を先取りすると、式(22)は

と変形することが出来ます。ここで、\(e\)は反発係数、\(\mathbf{v}_\text{w}(t)\)は壁の速度ベクトル、\(\mathbf{n}\)は壁の法線ベクトルであり、\(\mathbf{n}\)は

と表されます。
実際に数値計算を行う際には、衝突時刻\(t_i\)と位置\(\mathbf{r}=\mathbf{r}_i\)を求めた後、

に従って速度ベクトルを変更すると良いと思います（根拠は特にありません）。

多次元の反発

衝突時刻\(t=t_i\), 点\(\mathbf{r}(t_i)\)で質点が\(f(\mathbf{r},t)=0\)で表される壁に衝突することを考えます。
衝突前、後の質点の速度ベクトルをそれぞれ\(\mathbf{v}(t_i-0), \mathbf{v}(t_i+0),\)と置きます。
位置\(\mathbf{r}_i\)における壁の単位法線ベクトルを\(\mathbf{n}, (\mathbf{n}^2=1)\)と書くと、

と書くことが出来ます。ここで、\(c\)は未知の定数でこれから定めていきます。
衝突では壁の法線方向の速度成分のみが変化すると仮定しているので、

が成立します。ここで、\(e\)は反発係数で\(0\leq e\leq 1\)である(\(e=0\):完全非弾性衝突、\(e=1\):完全弾性衝突。
また、\(v_\perp, v_{\text{w}\perp}\)は\(\mathbf{v}, \mathbf{v}_\text{w}\)の、壁の法線方向の成分であり、

と書きあらわすことが出来ます。
式(27)の\(c\)を求めましょう。式(27)の両辺に\(\mathbf{n}^\mathsf{T}\)を掛けて

のように得ます。よって衝突後の質点の速度は

と表すことが出来ます。壁が\(f(\mathbf{r},t)=0\)を満たす線と表現されていれば、時刻\(t=t_i\)で\(\mathbf{n}\)は

と求める事が出来ます。
続いて壁面の速度\(\mathbf{v}_\text{w}(t)\)を考えます。1次元の場合と同様に衝突前後のごく短い時間では質点の動きは壁に追従すると予想します。すなわち、

が成立すると考えます。ここで、式(35)の\(\mathbf{r}(t)\)は質点の位置ベクトルです。すなわち、衝突が起こる時刻\(t=t_i\)の周りで

が成立しており、衝突のごく短い時間では

が得られます。よって

を得ますので、壁の速度

を得ることが出来ます。

プログラム

2次元平面の運動に対するFortran90のプログラムはこちらです。
時間発展は、刻み幅制御陽的ルンゲクッタ法
衝突時刻を求める際の根の探索には、Anderson-Björk’s法を用いています。

壁の関数(サブルーチン fw)とその偏微分(サブルーチン pwf)は手で入力しています。

壁との衝突は、関数の”符号が変わった時”で判定しているので、法線の方向はどちらでも構いません。
すなわち、
\(
f(x,y,t)=\pm g(x,y,t)
\)
は同じです。

プログラム

gnuplot用のスクリプト

実行結果

プログラムを実行した結果です。
適当な壁を定義して実行しています。

重力なしの運動、動かない壁

楕円の焦点から放たれた質点の軌跡

\(
\begin{align}
x(0)&=4,~~x'(0)=(\text{const}) \\
y(0)&=0,~~y'(0)=(\text{const}’)
\end{align}
\)

壁
\(
\begin{gather}
f(x,y,t)=\left(\frac{x}{l_x}\right)^2+\left(\frac{y}{l_y}\right)^2-1=0,\\
l_x=5,~l_y=3
\end{gather}
\)

重力下の運動、動かない壁

\(
\begin{align}
x(0)&=0,~~x'(0)=2 \\
y(0)&=0,~~y'(0)=7
\end{align}
\)

壁
\(
\begin{align}
f(x,y,t)=y – \sin(x)=0
\end{align}
\)

重力下の運動、動く壁

\(
\begin{align}
x(0)&=4,~~x'(0)=0 \\
y(0)&=4,~~y'(0)=5
\end{align}
\)

壁
\(
\begin{gather}
f(x,y,t)=y – \sin(kx+\phi_x)\sin(\omega t+\phi_t)=0,\\
k=1,~w=1,~\phi_x=-\pi/2,~\phi_t=0
\end{gather}
\)

参考文献

3. 壁との衝突 -物理学の見つけ方

4. 動く壁との衝突 -物理学の見つけ方