この記事の前の話はこちら

波動の力学モデル

記事のタイトルで「多体問題を解く」と言っているのにこれまで1体問題ばかり解いてきた。
いよいよ本題である多体問題を解いてみる。

テーマは「ばね質点系」である。
1次元の波動の力学モデルであり、この系に振動を入力すれば、波が発生するはずである。
理論モデルとしては単純明快で高校生でもわかるレベルのものだが、計算機で「実験」して、果たしてそのとおりになるのか。

ばね質点系

$N$個の質量 $m$ の質点が $x \hspace{0.1em}$軸上に配置されていて、 $i \hspace{0.1em}$番目の質点の位置が $x_{i} = i \hspace{0.1em} \Delta x$ （ $\Delta x$ は定数）であるとする。

隣り合う質点の間は、すべて、ばね定数 $k$ のばねで繋がれている。

横波の力学モデル

この系において、 $i \hspace{0.1em}$番目の質点の $y \hspace{0.1em}$座標を $y_{i}$ とする。
すべての質点は $x \hspace{0.1em}$方向と $z \hspace{0.1em}$方向には変位しないとする。
ばねの自然長を $L$ とする。

$i \hspace{0.1em}$番目の質点と $( \hspace{0.1em} i - 1 \hspace{0.1em} ) \hspace{0.1em}$番目の質点の間のばねの長さを $r_{i}$ とする。
$r_{i}$ は、 $y_{i} - y_{i - 1}$ を $\Delta y_{i}$ とおけば
$$
r_{i} = \sqrt{\Delta x^{2} + {\Delta y_{i}}^{2}}
$$
である。

このばねの弾性力の大きさ $f_{i}$ は
$$
f_{i} = k \hspace{0.1em} ( \hspace{0.1em} r_{i} - L \hspace{0.1em} )
$$
である。

このばねから $i \hspace{0.1em}$番目の質点が受ける力の $y \hspace{0.1em}$成分 $F_{iy-}$ は、 $\Delta y_{i}$ が正である場合に負の向きであることに注意すると
$$
F_{iy-} = - \frac{\Delta y_{i}}{r_{i}} \hspace{0.1em} f_{i}
$$
である。

いっぽう、 $( \hspace{0.1em} i + 1 \hspace{0.1em} ) \hspace{0.1em}$番目の質点との間のばねから $i \hspace{0.1em}$番目の質点が受ける力の $y \hspace{0.1em}$成分 $F_{iy+}$ は、 $\Delta y_{i+1}$ が正である場合に正の向きであることに注意すると
$$
F_{iy+} = + \frac{\Delta y_{i+1}}{r_{i+1}} \hspace{0.1em} f_{i+1}
$$
である。

$i \hspace{0.1em}$番目の質点が受ける合力の $y \hspace{0.1em}$成分 $F_{i}$ は
$$
F_{i} = F_{iy-} + F_{iy+}
$$
である。

$i \hspace{0.1em}$番目の質点の運動方程式の $y \hspace{0.1em}$成分は、時刻を $t$ として
$$
m \hspace{0.1em} \frac{d^{2} y_{i}}{dt^{2}} = F_{i}
$$
である。

以上をすべて合わせると
$$
m \hspace{0.1em} \frac{d^{2} y_{i}}{dt^{2}} = k \hspace{0.1em} \left( - \Delta y_{i} \hspace{0.1em} \left( \hspace{0.1em} 1 - \frac{L}{\sqrt{\Delta x^{2} + {\Delta y_{i}}^{2}}} \hspace{0.1em} \right) + \Delta y_{i+1} \hspace{0.1em} \left( \hspace{0.1em} 1 - \frac{L}{\sqrt{\Delta x^{2} + {\Delta y_{i+1}}^{2}}} \hspace{0.1em} \right) \hspace{0.1em} \right)
$$
が得られる。
これが、 $y \hspace{0.1em}$方向にだけ運動できるばね質点系の運動方程式である。

質点の $y \hspace{0.1em}$方向の変位が $\Delta x$ に比べてじゅうぶん小さい場合について考えると、 $\frac{\Delta y_{i}}{\Delta x}$ について1次近似し
$$
\frac{L}{\sqrt{\Delta x^{2} + {\Delta y_{i}}^{2}}} = \frac{1}{\Delta x} \hspace{0.1em} \frac{L}{\sqrt{1 + \frac{{\Delta y_{1}}^{2}}{\Delta x^{2}}}} = \frac{L}{\Delta x} \hspace{0.1em} \left( \hspace{0.1em} 1 - \frac{1}{2} \hspace{0.1em} \frac{{\Delta y_{i}}^{2}}{\Delta x^{2}} + \cdots \hspace{0.1em} \right) \simeq \frac{L}{\Delta x}
$$
$$
1 - \frac{L}{\sqrt{\Delta x^{2} + {\Delta y_{i}}^{2}}} \simeq \frac{\Delta x - L}{\Delta x}
$$
同様に
$$
1 - \frac{L}{\sqrt{\Delta x^{2} + {\Delta y_{i+1}}^{2}}} \simeq \frac{\Delta x - L}{\Delta x}
$$
とすることにより、運動方程式を
$$
m \hspace{0.1em} \frac{d^{2} y_{i}}{dt^{2}} = k \hspace{0.1em} ( - \Delta y_{i} + \Delta y_{i+1} ) \hspace{0.1em} \frac{\Delta x - L}{\Delta x}
$$
$$
\frac{m}{\Delta x} \hspace{0.1em} \frac{d^{2} y_{i}}{dt^{2}} = \Bigl( \hspace{0.1em} k \hspace{0.1em} ( \hspace{0.1em} \Delta x - L \hspace{0.1em} ) \hspace{0.1em} \Bigr) \hspace{0.1em} \frac{- \Delta y_{i} + \Delta y_{i+1}}{\Delta x^{2}}
$$
と変形することができる。

ここで、 $\frac{m}{\Delta x}$ を $\rho$ とおき、 $k \hspace{0.1em} ( \hspace{0.1em} \Delta x - L \hspace{0.1em} )$ を $S$ とおくと
$$
\rho \hspace{0.1em} \frac{d^{2} y_{i}}{dt^{2}} = S \hspace{0.1em} \frac{\frac{\Delta y_{i+1}}{\Delta x} - \frac{\Delta y_{i}}{\Delta x}}{\Delta x}
$$
となる。
$\Delta x \to 0$ の極限で、 $\frac{\Delta y_{i}}{\Delta x}$ は $y$ の $x$ に対する偏導関数 $\frac{\partial y}{\partial x}$ になり、 $\frac{\frac{\Delta y_{i+1}}{\Delta x} - \frac{\Delta y_{i}}{\Delta x}}{\Delta x}$ は $\frac{\partial y}{\partial x}$ の $x$ に対する偏導関数 $\frac{\partial^{2} y}{\partial x^{2}}$ になる。
したがって、運動方程式は
$$
\rho \hspace{0.1em} \frac{\partial^{2} y}{\partial t^{2}} = S \hspace{0.1em} \frac{\partial^{2} y}{\partial x^{2}}
$$
なる偏微分方程式となる。

これは、速さ $v = \sqrt{\frac{S}{\rho}}$ で伝わる波を解とする波動方程式である。

ばね質点系で $\Delta x \to 0$ とする極限は、質点の列の間隔が無限に小さくなり連続とみなせるようになった状態、つまり連続体としての弦である。
それを踏まえると、 $\rho$ は弦の線密度、 $S$ は弦の張力であると考えられる。

理屈

$N \hspace{0.1em}$個の質点から成るばね質点系において、 $i \hspace{0.1em}$番目の質点の質量を $m_{i}$ とする。

$i \hspace{0.1em}$番目の質点の位置の座標を ${\mathbf r_{i}} = ( \hspace{0.1em} x_{i} \hspace{0.1em} , \hspace{0.1em} y_{i} \hspace{0.1em} , \hspace{0.1em} z_{i} \hspace{0.1em} )$ とする。

$i \hspace{0.1em}$番目の質点と $j \hspace{0.1em}$番目の質点との間がばね定数 $k_{ij}$ 、自然長 $L_{ij}$ のばねで繋がれているとする。

質点の運動方程式

$i \hspace{0.1em}$番目の質点の運動方程式は
$$
m_{i} \hspace{0.1em} \frac{d^{2} {\mathbf r_{i}}}{dt^{2}} = - \sum_{j \ne i} k_{ij} \hspace{0.1em} \left( \hspace{0.1em} {\mathbf r_{i}} - {\mathbf r_{j}} \hspace{0.1em} \right)
$$
である。

隣の質点とだけ繋がれている状況

上の式は $N \hspace{0.1em}$個の質点の中のあらゆる $2 \hspace{0.1em}$個の組合せの間にばねが存在するという状況を表すが、今回は、上の『横波の力学モデル』で説明した状況をシミュレートするため、 $i$ と $j$ の差が 1 以外である組については $k_{ij} = 0$ とする。

自由端・固定端

すると、端の $2 \hspace{0.1em}$個を除くすべての質点について、両側の $2 \hspace{0.1em}$つのばねから力を受ける形の運動方程式になるが、 $1 \hspace{0.1em}$番目と $N \hspace{0.1em}$番目の質点だけについては、 $１ \hspace{0.1em}$つのばねからしか力を受けない形の運動方程式になる。

これらの端の質点について、制約をなにも与えず、ただこの運動方程式に従うとするとき、これらの端は 自由端 である。

いっぽう、これらの端の質点が一切変位しないように固定されているとき、すなわち、たとえば $N \hspace{0.1em}$番目の質点について ${\mathbf r_{N}} = {\mathbf r_{N,0}}$ （ ${\mathbf r_{N,0}}$ は初期状態における $N \hspace{0.1em}$番目の質点の位置の座標）がつねに満たされるようになっているとき、その端は 固定端 である。

外力

特定の質点（ $p \hspace{0.1em}$番目とする）に外力 ${\mathbf F_{ex}}$ を加えることにより、それが波源となって系に波が発生すると考えられる。
その場合、 $p \hspace{0.1em}$番目の質点の運動方程式が
$$
m_{p} \hspace{0.1em} \frac{d^{2} {\mathbf r_{p}}}{dt^{2}} = - \sum_{j \ne p} k_{pj} \hspace{0.1em} \left( \hspace{0.1em} {\mathbf r_{p}} - {\mathbf r_{j}} \hspace{0.1em} \right) + {\mathbf F_{ex}}
$$
となる。

さらに、 ${\mathbf F_{ex}}$ が正弦波形で振動する場合、すなわち、角振動数を $\omega$ とし、 ${\mathbf F_{0}}$ を定数のベクトルとして
$$
{\mathbf F_{ex}} = {\mathbf F_{0}} \hspace{0.1em} \sin \bigl[ \hspace{0.1em} \omega \hspace{0.1em} t \hspace{0.1em} \bigr]
$$
である場合には、系に発生する波が正弦波となるはずである。

固有振動

両端が自由端の場合と両端が固定端の場合には、この系の端から端までが半波長の整数倍になるような波は、自由端を腹とし固定端を節とするような定常波になり、振幅が大きくなる。

両端が自由端と固定端の組み合わせである場合は、この系の端か
ら端までが4分の1波長の奇数倍であるような波は、自由端を腹とし固定端を節とするような定常波になり、振幅が大きくなる。

上の 2つに該当しないような波は、この系ではすぐに減衰してしまい、振幅が大きくならない。

上に述べた、振幅が大きくなる（系の固有振動に共鳴する）波の角振動数は、次のように定まる。
例として、両端が自由端と固定端の組み合わせの場合で、系の $M \hspace{0.1em}$倍振動（ $M$ は奇数）に共鳴する波を考える。

系の端から端までの長さを $D$ とすると、波の波長 $\lambda$ は
$$
\lambda = \frac{4}{M} \hspace{0.1em} D
$$
である。

系を波が伝わる速さを $v$ とすると、この波の振動数 $\nu$ は
$$
\nu = \frac{v}{\lambda}
$$
である。

この波の角振動数 $\omega$ は
$$
\omega = 2 \hspace{0.1em} \pi \hspace{0.1em} \nu
$$
である。

$v$ は、
質点の質量がすべて等しく $m$ であり、
ばねのばね定数と自然長がすべてのばねで等しくそれぞれ $k$ , $L$ であれば
$$
v = \sqrt{\frac{k \hspace{0.1em} ( \hspace{0.1em} \Delta x - L \hspace{0.1em} ) \hspace{0.1em} \Delta x}{m}}
$$
である。

以上をすべて合わせると
$$
\omega = \frac{\pi \hspace{0.1em} M}{2 \hspace{0.1em} D} \hspace{0.1em} \sqrt{\frac{k \hspace{0.1em} ( \hspace{0.1em} \Delta x - L \hspace{0.1em} )
\hspace{0.1em} \Delta x}{m}}
$$
が得られる。

プログラム

運動方程式の関数 `f`

f <- function(f.t,f.rv) {

f.deltaxij <- matrix(rep(f.rv[,"x"],times=n),nrow=n,ncol=n,byrow=T) - matrix(rep(f.rv[,"x"],times=n),nrow=n,ncol=n,byrow=F)
f.deltayij <- matrix(rep(f.rv[,"y"],times=n),nrow=n,ncol=n,byrow=T) - matrix(rep(f.rv[,"y"],times=n),nrow=n,ncol=n,byrow=F)
f.deltazij <- matrix(rep(f.rv[,"z"],times=n),nrow=n,ncol=n,byrow=T) - matrix(rep(f.rv[,"z"],times=n),nrow=n,ncol=n,byrow=F)

f.deltarij <- sqrt( f.deltaxij^2 + f.deltayij^2 + f.deltazij^2 )

f.expandij <- f.deltarij - natural.length

f.Fij <- spring.constant * f.expandij

f.Fijx <- f.Fij * ( - f.deltaxij ) / f.deltarij
f.Fijy <- f.Fij * ( - f.deltayij ) / f.deltarij
f.Fijz <- f.Fij * ( - f.deltazij ) / f.deltarij

f.Fx <- colSums(f.Fijx,na.rm=TRUE)
f.Fy <- colSums(f.Fijy,na.rm=TRUE)
f.Fy[point.of.action] <- f.Fy[point.of.action] + F0 * sin( omega * f.t )
f.Fz <- colSums(f.Fijz,na.rm=TRUE)

return( cbind( x=f.rv[,"vx"] , y=f.rv[,"vy"] , z=f.rv[,"vz"] , vx=f.Fx/m , vy=f.Fy/m , vz=f.Fz/m ) )

}

n : 質点の個数 $\hspace{0.5em} N$
natural.length : $i \hspace{0.1em}$番目の質点と $j \hspace{0.1em}$番目の質点の間のばねの自然長 $\hspace{0.5em} L_{ij} \hspace{0.5em}$ を格納するマトリックス
spring.constant : $i \hspace{0.1em}$番目の質点と $j \hspace{0.1em}$番目の質点の間のばねのばね定数 $\hspace{0.5em} k_{ij} \hspace{0.5em}$ を格納するマトリックス
point.of.action : 外力を加える質点の番号
F0 : 正弦波形の外力の振幅
omega : 正弦波形の外力の角振動数
m : 質点の質量を格納するベクトル
f.t : 時刻 $\hspace{0.5em} t$
f.rv : 質点の位置 ${\mathbf r}$ と速度 ${\mathbf v}$ を格納するデータフレーム
$\hspace{1em}$ 列 x y z : ${\mathbf r}$ の $x$成分 , $y$成分 , $z$成分
$\hspace{1em}$ 列 vx vy vz : ${\mathbf v}$ の $x$成分 , $y$成分 , $z$成分
f.deltaxij : $i \hspace{0.1em}$番目の質点の $x \hspace{0.1em}$座標 $x_{i}$ と $j \hspace{0.1em}$番目の質点の $x \hspace{0.1em}$座標 $x_{j}$ との差を格納するマトリックス
$\hspace{2em}$ $x_{j} - x_{i}$ の値が $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に入る。
f.deltayij : $i \hspace{0.1em}$番目の質点の $y \hspace{0.1em}$座標 $y_{i}$ と $j \hspace{0.1em}$番目の質点の $y \hspace{0.1em}$座標 $y_{j}$ との差を格納するマトリックス
$\hspace{2em}$ $y_{j} - y_{i}$ の値が $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に入る。
f.deltazij : $i \hspace{0.1em}$番目の質点の $z \hspace{0.1em}$座標 $z_{i}$ と $j \hspace{0.1em}$番目の質点の $z \hspace{0.1em}$座標 $z_{j}$ との差を格納するマトリックス
$\hspace{2em}$ $z_{j} - z_{i}$ の値が $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に入る。
f.deltarij : $i \hspace{0.1em}$番目の質点の座標 ${\mathbf r_{i}}$ と $j \hspace{0.1em}$番目の質点の座標 ${\mathbf r_{j}}$ との差の大きさを格納するマトリックス
$\hspace{2em}$ $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に値が入る。
f.expandij : $i \hspace{0.1em}$番目の質点と $j \hspace{0.1em}$番目の質点の間のばねの伸びを格納するマトリックス
$\hspace{2em}$ $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に値が入る。
f.Fij : $i \hspace{0.1em}$番目の質点と $j \hspace{0.1em}$番目の質点の間のばねの弾性力の大きさを格納するマトリックス
$\hspace{2em}$ $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に値が入る。
f.Fijx : $i \hspace{0.1em}$番目の質点との間のばねが $j \hspace{0.1em}$番目の質点におよぼす力の $x \hspace{0.1em}$成分を格納するマトリックス
$\hspace{2em}$ $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に値が入る。
f.Fijy : $i \hspace{0.1em}$番目の質点との間のばねが $j \hspace{0.1em}$番目の質点におよぼす力の $y \hspace{0.1em}$成分を格納するマトリックス
$\hspace{2em}$ $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に値が入る。
f.Fijz : $i \hspace{0.1em}$番目の質点との間のばねが $j \hspace{0.1em}$番目の質点におよぼす力の $z \hspace{0.1em}$成分を格納するマトリックス
$\hspace{2em}$ $i \hspace{0.1em}$行 $j \hspace{0.1em}$列に値が入る。
f.Fx : 質点が受ける合力の $x \hspace{0.1em}$成分を格納するベクトル
f.Fy : 質点が受ける合力の $y \hspace{0.1em}$成分を格納するベクトル
f.Fz : 質点が受ける合力の $z \hspace{0.1em}$成分を格納するベクトル

境界条件の関数 `boundary.condition`

boundary.condition <- function(boundary.condition.t,boundary.condition.rv) {

boundary.condition.rv[,"x"] <- x0
boundary.condition.rv[n,"y"] <- 0.0
boundary.condition.rv[n,"z"] <- 0.0

return( boundary.condition.rv )

}

x0 : 質点の初期位置の $x \hspace{0.1em}$座標を格納するベクトル

すべての質点について、 $x \hspace{0.1em}$方向には変位しないという条件を与えている。

また、 $N \hspace{0.1em}$番目の質点について、つねに初期位置にあるという条件を与えている。
つまり、 $N \hspace{0.1em}$番目の質点の側の端を、固定端としている。

なお、 $1 \hspace{0.1em}$番目の質点には、 $x \hspace{0.1em}$方向に変位しないという以外の
条件を特に与えていない。
つまり、 $1 \hspace{0.1em}$番目の質点の側の端を、自由端としている。

パラメータ

n <- 100
m0 <- 0.01
m <- rep(m0,times=n)
spring.constant.0 <- 100.0
spring.constant <- rbind(numeric(n),spring.constant.0*diag(n))[-(n+1),] + cbind(numeric(n),spring.constant.0*diag(n))[,-(n+1)]
diag(spring.constant) <- NA
natural.length.0 <- 0.1
natural.length <- matrix(rep(natural.length.0,times=n*n),nrow=n,ncol=n)
diag(natural.length) <- NA

point.of.action <- 80
omega <- 2 * pi * 1.1616
F0 <- 0.1

t0 <- 0.0
x0 <- seq(1.0,length=n)
y0 <- rep(0.0,times=n)
z0 <- rep(0.0,times=n)
vx0 <- rep(0.0,times=n)
vy0 <- rep(0.0,times=n)
vz0 <- rep(0.0,times=n)

tend <- 20.0

Deltat <- 2.0e-5
draw.step <- 1000
xmin <- 1
xmax <- 100
ymin <- - 0.25
ymax <- 0.25

n : 質点数今回は 100
m0 : 今回はすべて等しい質量の質点とする。その値。 0.01 とした
m : 質点の質量 $\hspace{0.5em} m_{i} \hspace{0.5em}$ を格納するベクトル。すべての要素を m0 に等しくした
spring.constant.0 : 今回はすべて等しいばね定数のばねとする。その値。 100 とした
spring.constant : $i \hspace{0.1em}$番目の質点と $j \hspace{0.1em}$番目の質点の間のばねのばね定数 $\hspace{0.5em} k_{ij} \hspace{0.5em}$ を格納するマトリックス。 $i$ と $j$ の差が 1 の要素は spring.constant.0 に等しく、 $i$ と $j$ の差が 1 より大きい要素は 0 。 $i$ と $j$ が等しい要素（対角成分）には 0 や他の数値を入れると後でエラーの元となりかねないため NA を入れた
natural.length.0 : 今回はすべて等しい自然長のばねとする。その値。 0.1 とした
natural.length : $i \hspace{0.1em}$番目の質点と $j \hspace{0.1em}$番目の質点の間のばねの自然長 $\hspace{0.5em} L_{ij} \hspace{0.5em}$ を格納するマトリックス。 $i$ と $j$ が異なる要素（対角成分以外）はすべて natural.length.0 に等しくした。 $i$ と $j$ が等しい要素（対角成分）には 0 や他の数値を入れると後でエラーの元となりかねないため NA を入れた

point.of.action : 外力を加える質点の番号。 80 とした
omega : 正弦波形の外力の角振動数
F0 : 正弦波形の外力の振幅。 0.1 とした。

t0 : 初期状態の時刻
x0 y0 z0 : 初期状態での質点の座標。 x0 は [ 1 , 2 , 3 , $\cdots$ ] とした（$\hspace{0.1em}$間隔 1 ずつ離して 1 から 100 まで）。 y0 と z0 はすべての要素が 0
vx0 vy0 vz0 : 初期状態での質点の速度の各成分。いずれもすべての要素が 0
tend : 終状態の時刻

Deltat : Runge-Kutta法の計算において1ステップで進める時間
draw.step : Runge-Kutta法で何ステップか計算するたびに画像を出力する、そのステップ数

結果を出力するコード

output <- function() {

plot(rv[,"x"],rv[,"y"],xlim=c(xmin,xmax),ylim=c(ymin,ymax),xlab="x",ylab="y",tcl=0.5,las=1,mgp=c(2.5,0.3,0),text(85,0.23,labels=paste("t =",as.character(round(t,digits=3)))))

}

今回はエネルギーの出力はなし。
それ以外は振り子や天体のときと内容的には同じ。

全リスト

################################################################
## parameters
## modify this section to change condition
################################################################

n <- 100
m0 <- 0.01
m <- rep(m0,times=n)
spring.constant.0 <- 100.0
spring.constant <- rbind(numeric(n),spring.constant.0*diag(n))[-(n+1),] + cbind(numeric(n),spring.constant.0*diag(n))[,-(n+1)]
diag(spring.constant) <- NA
natural.length.0 <- 0.1
natural.length <- matrix(rep(natural.length.0,times=n*n),nrow=n,ncol=n)
diag(natural.length) <- NA

A <- 0.1
point.of.action <- 80
omega <- 2 * pi * 1.3939
F0 <- 0.1

t0 <- 0.0
x0 <- seq(1.0,length=n)
y0 <- rep(0.0,times=n)
z0 <- rep(0.0,times=n)
vx0 <- rep(0.0,times=n)
vy0 <- rep(0.0,times=n)
vz0 <- rep(0.0,times=n)

tend <- 20.0

Deltat <- 2.0e-5
draw.step <- 1000
xmin <- 1
xmax <- 100
ymin <- - 2.5 * A
ymax <- 2.5 * A


################################################################
## output
## modify this section to change the method to output
################################################################

output <- function() {

plot(rv[,"x"],rv[,"y"],xlim=c(xmin,xmax),ylim=c(ymin,ymax),xlab="x",ylab="y",tcl=0.5,las=1,mgp=c(2.5,0.3,0),text(85,0.23,labels=paste("t =",as.character(round(t,digits=3)))))

}


################################################################
## models
## modify this section to change physical model
################################################################

f <- function(f.t,f.rv) {

f.deltaxij <- matrix(rep(f.rv[,"x"],times=n),nrow=n,ncol=n,byrow=T) - matrix(rep(f.rv[,"x"],times=n),nrow=n,ncol=n,byrow=F)
f.deltayij <- matrix(rep(f.rv[,"y"],times=n),nrow=n,ncol=n,byrow=T) - matrix(rep(f.rv[,"y"],times=n),nrow=n,ncol=n,byrow=F)
f.deltazij <- matrix(rep(f.rv[,"z"],times=n),nrow=n,ncol=n,byrow=T) - matrix(rep(f.rv[,"z"],times=n),nrow=n,ncol=n,byrow=F)

f.deltarij <- sqrt( f.deltaxij^2 + f.deltayij^2 + f.deltazij^2 )

f.expandij <- f.deltarij - natural.length

f.Fij <- spring.constant * f.expandij

f.Fijx <- f.Fij * ( - f.deltaxij ) / f.deltarij
f.Fijy <- f.Fij * ( - f.deltayij ) / f.deltarij
f.Fijz <- f.Fij * ( - f.deltazij ) / f.deltarij

f.Fx <- colSums(f.Fijx,na.rm=TRUE)
f.Fy <- colSums(f.Fijy,na.rm=TRUE)
f.Fy[point.of.action] <- f.Fy[point.of.action] + F0 * sin( omega * f.t )
f.Fz <- colSums(f.Fijz,na.rm=TRUE)

return( cbind( x=f.rv[,"vx"] , y=f.rv[,"vy"] , z=f.rv[,"vz"] , vx=f.Fx/m , vy=f.Fy/m , vz=f.Fz/m ) )

}


boundary.condition <- function(boundary.condition.t,boundary.condition.rv) {

boundary.condition.rv[,"x"] <- x0
boundary.condition.rv[n,"y"] <- 0.0
boundary.condition.rv[n,"z"] <- 0.0

return( boundary.condition.rv )

}


################################################################
## engines
## do not modify this section
################################################################

runge.kutta.4th <- function(runge.kutta.4th.t,runge.kutta.4th.rv) {

runge.kutta.4th.k0 <- f( runge.kutta.4th.t , runge.kutta.4th.rv ) * Deltat
runge.kutta.4th.k1 <- f( runge.kutta.4th.t + Deltat / 2.0 , runge.kutta.4th.rv + runge.kutta.4th.k0 / 2.0 ) * Deltat
runge.kutta.4th.k2 <- f( runge.kutta.4th.t + Deltat / 2.0 , runge.kutta.4th.rv + runge.kutta.4th.k1 / 2.0 ) * Deltat
runge.kutta.4th.k3 <- f( runge.kutta.4th.t + Deltat , runge.kutta.4th.rv + runge.kutta.4th.k2 ) * Deltat

runge.kutta.4th.Deltarv <- ( runge.kutta.4th.k0 + 2.0 * ( runge.kutta.4th.k1 + runge.kutta.4th.k2 ) + runge.kutta.4th.k3 ) / 6.0

return( runge.kutta.4th.rv + runge.kutta.4th.Deltarv )

}


### preparation to output

png(filename="%04d.png",width=960,height=960,units="px",res=192)

###

t <- t0
rv <- cbind( x=x0 , y=y0 , z=z0 , vx=vx0 , vy=vy0 , vz=vz0 )

while(t<=tend) {

output()

for(i in 1:draw.step) {

rv <- runge.kutta.4th(t,rv)

rv <- boundary.condition(t,rv)

t <- t + Deltat

}

}


### post-processing to outuput

graphics.off()

###

実行結果

外力を加える質点（波源）の位置は $x = 80$ 。
波源の振動数を2通りに変えてやってみた。

パターン1

1つめは
$$
\omega = 2 \hspace{0.1em} \pi \times 1.3939
$$
つまり、 $f = 1.3939$ の場合。

結果は次のようになった。

パターン2

2つめは
$$
\omega = 2 \hspace{0.1em} \pi \times 1.1616
$$
つまり、 $f = 1.1616$ の場合。

結果は次のようになった（ファイルサイズの制限のために途中までで切れている）。

考察

$ t = 0 $ から後、波源から波が発生し左右に伝わっていっている。
このモデルが確かに横波の伝播のモデルとして成立していることがわかる。

波が左端（自由端）まで伝わると、端の近くで反射による定常波が生じるのがわかる。
パターン１・２ともそうであるが、とくに、パターン１で短い間だけ現れてすぐ消えてしまうのが見られるのが面白い。
右端（固定端）でも同じことが言えるはずだが、そうだと思って見ればそう見える、かな…。なんとも言えない。

パターン１では最初の僅かな間だけそのような定常波が見えるが、そのあとは、波源からの波や端からの反射波が入り乱れて、系全体の振動は非常に複雑な変化をする。
全般的にあまり大きな波が現れず、ほとんど波が存在しない時間帯も多く、波が現れても細かい動きをしてすぐに消えてしまったり、はっきりしない。
波源の振動数が系の固有振動数と一致していないからだと考えられる。

パターン２は、はっきりと系の固有振動である。
左端を腹、右端を節とする固有振動で、５倍振動になっている。
この系の５倍振動の振動数を理論的に計算してみる。
『理屈』に書いたように、系の固有角振動数 $\omega$ は
$$
\omega = 2 \hspace{0.1em} \pi \times \frac{M}{4 \hspace{0.1em} D} \hspace{0.1em} \sqrt{\frac{k \hspace{0.1em} ( \hspace{0.1em} \Delta x - L \hspace{0.1em} )
\hspace{0.1em} \Delta x}{m}}
$$
である。これに
$$
M = 5
\hspace{0.5em} , \hspace{0.5em}
D = 99
\hspace{0.5em} , \hspace{0.5em}
k = 100
\hspace{0.5em} , \hspace{0.5em}
\Delta x = 1
\hspace{0.5em} , \hspace{0.5em}
L = 0.1
\hspace{0.5em} , \hspace{0.5em}
m = 0.01
$$
を代入すると、 $\omega = 2 \hspace{0.1em} \pi \times 1.1978$ である。
パターン２は $\omega = 2 \hspace{0.1em} \pi \times 1.1616$ であるから、近いが、ちょっとずれがある。
ずれは 2% 程度である。
一致したと言えなくもない。どうだろうか。

パターン２では波の振幅がかなり大きくなるので、波動方程式の導出のさいに利用した $y_{i} \ll \Delta x$ という近似が成り立たなくなっているはずだ。
だから固有角振動数が上の理論式からずれるのは当然と言えよう。
定性的には次のように考えられる。
$\frac{L}{\sqrt{\Delta x^{2} + {\Delta y_{i}}^{2}}}$ を $\frac{L}{\Delta x}$ に近似していたのができなく（$ \hspace{0.1em} \Delta y_{i}$ の影響を無視できなく）なるわけであるから、近似できていたときに比べて $\left( 1 - \frac{L}{\sqrt{\Delta x^{2} + {\Delta y_{i}}^{2}}} \right)$ が大きくなるということである。
これは、質点の運動方程式の右辺が大きくなることを意味する。
すなわち、実際の張力よりも強い張力の場合の運動方程式に近くなるということである。
よって、波の速さ $v$ は理論式よりも大きくなり、ゆえに、固有角振動数 $\omega$ も理論式より大きくなる。

以上の考察が的を得ているかどうかは、他の固有振動を作ってみることにより確かめられるだろう。
まず理論式にある程度近い結果になるかどうかを調べることが必要だし、近いようなら、この５倍振動と同様の傾向でずれが生じるかどうかを調べたいところである。
とくに、基本振動・３倍振動については、比較的簡単に作れる可能性が高そうであるから、早速に試してみたい課題である。