力学アルゴリズムによる波のシミュレーション（線上の波）

このプログラムは、線上を伝わる波を、単純な力学アルゴリズムと数値積分法によって計算し、 その結果をリアルタイムでアニメーション表示するVCSSLログラムです。

スポンサーリンク

---

使用方法

ダウンロードと展開（解凍）

まず、PC（スマホは未対応）で上の画面の「 ダウンロード 」ボタンを押してください。 するとZIP形式で圧縮されたファイルがダウンロードされます。

その後、ZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開（解凍）してください。 展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。

» 展開がエラーで止まってしまう場合や、ファイル名が文字化けしてしまう場合は…

プログラムの起動

Windows をご使用の場合

上記でZIPファイルを展開したフォルダ内にある、以下のバッチファイルをダブルクリック実行してください：

もしプログラムを書き変えながら使いたい場合は、代わりに「 VCSSL_Editor__プログラム編集はこちら.bat 」を実行してください。

正常に起動できると、初回のみ、Java実行環境を入手するか等を尋ねられるので、適時答えて済ませると、プログラムが起動します。 2回目以降はすぐに起動します。

» うまく起動できずにエラーになってしまう場合は…

Linux 等をご使用の場合

ZIPファイルを展開したフォルダ内へコマンドライン端末で cd して、以下の通り入力して実行してください：

» javaコマンドが使用できない等のエラーが表示される場合は…

操作方法

プログラムを起動すると、線上を伝わる波の様子が、グラフ上にアニメーションで表示されます。

このプログラムはアルゴリズムの解説用であるため、波のパラメータを操作するような設定画面はありません。

波のパラメータを変えたい場合は、プログラム内の先頭領域に定義されている数値を直接書き換えてください。

題材解説

線上を伝わる波

「波」というと、一般的には水面の波など、面上を伝わる波（2次元波）を想像します。

しかし、波の中で最も基本的と言えるのは、線上を伝わる波(1次元波)です。これは例えば、ギターの弦の振動や、長く伸ばされたコイルばねの振動などに見られます。

シミュレーションでは、バネと質点でモデル化

波を伝えるのは、弾性のある、つまり伸び縮みする媒質です。

このような弾性体をシミュレーションするにはいくつかの方法があります。簡単なのは、媒質を、バネで繋がった質点（質量のある点）で置き換える事です。そして弾性体の張力はバネの張力で、密度は質点の質量で表現します。

このようにモデルを設定すれば、媒質全体の運動は、左右のバネから力を受ける質点の運動が、横に連なったものとして扱えます。

右のバネから質点が受ける力

個々の質点の運動は、力さえ分かれば、ニュートンの運動方程式を使ってシミュレーションできます。なので、あと必要なのは、個々の質点が受ける力の値です。

個々の質点は、左右のバネから受ける力の和です。簡単のために、右のバネからのみ受ける力を考えましょう。下図のように、質点間の X 方向の距離を dx、Y 方向の距離を dy と取ります。

この通り、バネの張力を T , バネの仰角（X軸とのなす角）を θ とすると、質点は水平に T cos θ, 垂直に T sin θ の大きさの力を受けます。各方向の力を Fx , Fy と書くと：

\( F_{x, 右のバネ} = T \cos{ \theta } \)　　　… (1)
\( F_{y, 右のバネ} = T \sin{ \theta } \)　　　… (2)

張力 T については、厳密にはバネが伸びる影響も聞いてくるので、区間に依存するように思えます。しかし、今回はこの影響は考えません。 即ち、バネはある程度の強さの初期張力を持っていて、波でバネが伸び縮みする影響は、初期張力に比べてほぼ無視できるオーダーという近似を行います。 つまるところ、T は全域で一様な定数と見なします。

さらに、θがあまり大きくない（微小振動）という近似も用います。すると、sin と cos も以下のように近似されます。

\( \cos{ \theta } = 1 + O(\theta^2) \)
\( \sin{ \theta } = \tan{ \theta } + O(\theta^2)\)

O(θ^2)というのは、θについて2次以上の形で聞いてくる量という事です。θが微小なとき、この量は 1 や tanθ に比べて、ずっと小さいのでほぼ無視できます。

tanθ は dy / dx で置き換えられます。2次以上の項を無視すると、質点が受ける力は、

\( F_{x, 右のバネ} \simeq T \)　　　… (3)
\( F_{y, 右のバネ} \simeq T dy/dx = T (y_{i+1} - y_{i})/dx \)　　　… (4)

と近似的に求まります。

左右のバネから受ける力の合計

さて、質点は左右のバネから引っ張られます。左右のバネでは引っ張られる力の向きが反対なので、符号を反転させて足し合わせると、合計の力が求まります。

ここで水平方向の力 Fx は打ち消しあってしまい、θの2次以上の項しか残りません。つまり質点は、微小振動であれば水平方向にはほとんど動かないと近似できます。なので、これ以降はdx を全域でほぼ一様な定数と近似します：

\( dx \simeq L / (N-1) \)

ここで L は連続体媒質の全長、N はそれを近似するバネの個数を意味します。

垂直方向の力 Fy は、

\( F_y = F_{y, 右のバネ} - F_{y, 左のバネ} \)
\( \simeq T (y_{i+1} - y_{i})/dx - T (y_{i-1} - y_{i})/dx \)
\( = T (y_{i+1} + y_{i-1} - 2 y_i)/dx \)　　　… (5)

と求まります。

端の処理（固定端、自由端）

上で求めた式は、左右が質点ある事を前提としているので、両端の質点には適用できません。この、両端の質点をどのように処理するかについては、2通りの方式があります。

一つは、固定端と言って、両端を固定してしまう方式です。これは、両端の質点に働く力を、常に 0 にする事によってシミュレーションします。 固定端では、波は逆位相となって反射します（いわゆる固定端反射）。

\( F_{y, 右端, 固定端} = 0 \)
\( F_{y, 左端, 固定端} = 0 \)

もう一つは、自由端と言って、両端を何も固定しない方式です。この場合、両端の質点に働く力は、普通に一個手前の質点と繋がるバネから求めます。例えば右端の質点なら、左側のバネからの力のみを求めます（というより右にバネは無い）。 自由端では、波は同位相で反射します（いわゆる自由端反射）。

\( F_{y, 右端, 自由端} = F_{y, 左のバネ} = T (y_{右端の隣} - y_{右端})/dx \)
\( F_{y, 左端, 自由端} = F_{y, 右のバネ} = T (y_{左端の隣} - y_{左端})/dx \)

速度に比例する摩擦力で、減衰効果を加える

このままでは、全体でエネルギーが保存しているため、波がいつまでたっても消えません。しかし現実では、波はいずれ減衰して消えてしまいます。

このような減衰効果は、モデルによってメカニズムが異なります。最も簡単なのは以下のように、各質点に、速度に比例する摩擦力を加える事です。

\( F_y = T (y_{i+1} + y_{i-1} - 2 y_i)/dx - f v_i \)　　　… (6)

小文字の v は質点の速度です。小文字の f は摩擦力の係数で、この値が大きいほど波の減衰が強くなります。摩擦力は、速度と逆向きに作用するため、項の符合にマイナスが付いています。

数値積分によって運動を解く

個々の質点の受ける力が求まったので、あとはニュートンの運動方程式を立てて数値積分すれば、質点の運動が求まります。個々の質点の運動方程式は、加速度を a 、質点の質量を m として：

\( m a = F \)

と置けます。ここで加速度 a は速度 v の微分 dx /dt であり、さらに速度 v は位置 y の微分 dy / dt なので、以下の 2 式の形に直します。

\( dv / dt = F / m \) 　　… (7)
\( dx / dt = v \) 　　… (8)

あとはこれらの連立微分方程式を解けば良いのですが、シミュレーションで数値的に解くために、時間変数 t を細かい区間 Δt に区切り、数値積分法を用いて解きます。

最も単純なのは、Euler(オイラー)法と呼ばれる数値積分法です。 これは、時間区間 Δt 内では、v や F がほぼ変わらないと見なす方法です。こうすると、

\( \Delta v / \Delta t = ( v_{t+1} - v_t ) / \Delta t = F_t / m \) 　　… (9)
\( \Delta x / \Delta t = ( x_{t+1} - x_t ) / \Delta t = v_t \) 　　… (10)

となり、現在の時刻の状態から、次の時刻の状態を求める漸化式が得られます：

\( v_{t+1} = v_t + (F_t/m) \Delta t \) 　　… (11)
\( x_{t+1} = x_t + v_t \Delta t \) 　　… (12)

ただし、Euler法はあまり精度が良くありません。科学技術計算で要求される精度を出したければ、4次Runge-Kutta法(RK4)など、もっと高精度の数値積分法を用いる必要があります。

しかし、Euler法を次のように一部改変するだけで、RK4ほどではありませんが、精度を上げる事ができます：

\( v_{t+1} = v_t + (F_t/m) \Delta t \) 　　… (13)
\( x_{t+1} = x_t + v_{t+1} \Delta t \) 　　… (14)

このように、位置の更新に、古い速度ではなく、新しい速度を用いるだけです。こうしておく事で Symplectic数値積分の一種となり、エネルギーが振動的に保存するようになると共に、精度自体もEuler法より向上します。

---

スポンサーリンク

---

コード解説

コード全体

このプログラムのコード全体は、以下のようになっています。

以下では、コードの各部について解説します。

先頭領域

プログラムの先頭領域では、数学関数を扱う標準ライブラリ「 Math 」と、 2次元グラフを扱う拡張ライブラリ「 tool.Graph2D 」を読み込んでいます。

グローバル変数

その下では、グローバル変数を宣言しています。シミュレーションするモデルのパラメータなどを用意しています。

main関数

関数部分を見ていきましょう。まずはmain関数です。

main関数は、グローバル領域の処理（初期化など）が完了した時点で、システムから自動的にコールされます。つまり、プログラムの具体的な処理はここから始まります。

今回は、以下のようになっています。

このように、ユーザーに開始を伝えて、初期処理を行うinitialize関数（後述）をコールして初期化したあとは、メインループの処理を行うmainLoop関数へ処理を渡しています。

mainLoop関数

続いてmainLoop関数です。この関数はメインループの処理を行います。

メインループとは、プログラムの開始から終了まで、ずっと繰り返し回り続けている部分の事で、アニメーションやゲームでは必須です。毎秒数十回、モデルを少しずつ動かし、その度に画面を更新する事で、アニメーションしているように見せます。

while(true){ } の中は、プログラムの開始から終了までずっと回り続けます。これがメインループ部です。

これがメインループの中では、update関数（後述）でモデルの状態を少しだけ変化させています。

また、plotCounter回に一回、setGraph2DData関数でグラフソフトに座標値配列を転送し、グラフを描画させています。

グラフを描画した後は、グローバル変数 WAIT に指定されたミリ秒数だけプログラムを待機し、そして描画タイミングのカウンタ plotCounter をリセットしています。

メインループを一周するごとに、カウンタ plotCounter の値を 1 ずつ加算し、その値が plotTiming になったタイミングで、グラフを描画するという仕組みです。

initialize関数

その後では、モデルの状態に初期値を設定する initialize 関数を実装しています。

まずは頂点（質点の位置）のX座標配列です。今回のモデルなら、質点はX方向へは動かないのですが、グラフソフトに座標を転送する必要があるので用意しています。値は、モデルの長さをN等分した地点のX座標を入れています。

続いて頂点（質点の位置）のY座標配列です。今回は、ガウス関数による波束を設定しています。Y速度の配列では、静止状態からのスタートとなる 0 に設定しています。

ここまでは前半です。後半では、グラフソフトの起動や設定を行っています。

重要なのは newGraph2D 関数で、この関数は2Dグラフソフトを起動し、それに整数の識別番号「グラフID」を割り振って返します。それをグローバル変数 graph に格納しています。グラフを操作する時は、この graph を第一引数に指定します。

update関数

最後に、モデルの毎時刻の状態変化を処理する update 関数の実装です。この関数が今回のシミュレーションの中核となります。

最初に、定数である質点間水平距離 dx と、質点1個あたりの質量 m を求めています。この値は定数なのでグローバル領域に置いておいてもいいのですが、計算コスト的に大した節約にはならないので、可読性優先でここに置いています。

続いて、頂点（質点）にかかる力を求めています。両端の点については、自由端と固定端で分岐させて求めています。

最後に、位置と速度のデータを更新しています。これは題材解説で求めた数値積分の手続きそのままです。

ライセンス

このVCSSL/Vnanoコード（ 拡張子が「.vcssl」や「.vnano」のファイル ）は実質的な著作権フリー（パブリックドメイン) である CC0 の状態で公開しています※。 記事中にC言語/C++/Java言語などでのサンプルコードが掲載されいてる場合は、それらについても同様です。 そのままでのご利用はもちろん、改造や流用などもご自由に行ってください。

※ ただし、このコードの配布フォルダ内には、ダウンロード後すぐに実行できるように、 VCSSLの実行環境も同梱されており、そのライセンス文書は「 License 」フォルダ内に同梱されています （要約すると、商用・非商用問わず自由に使用できますが、使用の結果に対して開発元は一切の責任を負いません、といった具合の内容です）。 配布フォルダ内の各構成物の一覧やライセンスについては「 ReadMe_使用方法_必ずお読みください.txt 」をご参照ください。

※ Vnano の実行環境については、別途スクリプトエンジンのソースコードも一般公開しており、 何らかのソフトウェア内に組み込んでご利用いただく事も可能です。詳細はこちらをご参照ください。

この記事中の商標などについて

• OracleとJavaは、Oracle Corporation 及びその子会社、関連会社の米国及びその他の国における登録商標です。文中の社名、商品名等は各社の商標または登録商標である場合があります。
• Windows は、米国 Microsoft Corporation の米国およびその他の国における登録商標です。この記事は独立著作物であり、Microsoft Corporation と関連のある、もしくはスポンサーを受けるものではありません。
• Linux は、Linus Torvalds 氏の米国およびその他の国における商標または登録商標です。
• その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。