[ 前へ | 目次 | 次へ ]
Now Loading...
ダウンロード
PC (※スマートフォンでは動きません) でダウンロードし、ZIPファイルを右クリックメニューから展開して、できたフォルダ内の「 VCSSL.bat(バッチファイル) 」または「 VCSSL.jar 」をダブルクリックなどで実行すると、プログラムが起動します。
» 詳しい使用方法や、エラーで展開できない際の対応方法などはこちら

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

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

使用方法

ダウンロードと展開(解凍)

まず、PC(スマホは未対応)で上の画面の「 ダウンロード 」ボタンを押してください。 するとZIP形式で圧縮されたファイルがダウンロードされるので、そのZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開(解凍)してください。 展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。

» 展開がエラーで止まってしまう場合は…

なお、Linux® 等をご使用で、右クリックメニューから展開するとファイル名が文字化けしてしまう場合は、 コマンドライン端末でZIPファイルのある場所まで cd した上で「 unzip -O cp932 ZIPファイル名 」で展開してみてください。

プログラムの起動

Microsoft® Windows® をご使用の場合

上記の通りにZIPファイルを展開したフォルダ内にある、 「 VCSSL.bat(種類はバッチファイル) 」をダブルクリック実行してください。 もしプログラムの内容を書き変えながら使いたい場合は、代わりに「 VCSSL_Editor.bat 」を実行してください。

実行すると、最初にメモリー使用量や、(必要な場合のみ)Java®実行環境を自動で入手するか 等を尋ねられるので、適時答えると、プログラムが起動します。2回目以降はすぐに起動します。

※ ここで入手したJava®実行環境は、ZIPファイルを展開した中の「 jre 」フォルダ内にダウンロードされ、このプログラムの実行のみに使用されます。PC全体に影響する形でインストールされる事はありません。

Linux® 等やその他のOSをご使用の場合

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

java -jar VCSSL.jar
(プログラムの内容を書き変えながら使いたい場合は、代わりに VCSSL_Editor.jar を実行)

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

操作方法

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

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

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

題材解説

線上を伝わる波

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

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

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

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

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

バネで繋がった質点のモデル
モデルの概要

伸び縮みする連続的な媒質を、バネと質点(質量のある点)で置き換えて扱うモデルです。弾性体の張力はバネの張力で、密度は質点の質量で表現します。

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

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

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

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

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

質点間の X 方向の距離を dx、Y 方向の距離を dy と置きます。さらにバネの張力を T , バネの仰角を θ とすると、質点は水平に T cos θ, 垂直に T sin θ の大きさの力を受けます。

この通り、バネの張力を 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のプログラミングガイド(無料)はこちらへ!

上記のコードはプログラミング言語VCSSLで記述されており、VCSSLのプログラミングガイドは下記で無料公開しています。 上記のコードを改造したい方や、新しいコードを書いてみたい方はぜひご活用ください!

ブラウザで読めるWeb版だけでなく、PDF版も無料で配布しています!

スタートアップガイド( プログラミングがはじめての方向け )
プログラミングの入門書に相当する内容です。プログラミングが初めての方はこちらがおすすめです。
即席ガイド( C系言語ユーザー向け )
C言語や C++ などのC系の言語を扱われている方が、即席でVCSSLを扱うための簡易ガイドです。
文法ガイド
VCSSLの文法や基本的な機能を淡々とまとめた、リファレンスマニュアル的な位置づけのガイドです。
GUI開発ガイド
ボタンや入力項目などのGUI部品が並ぶ、画面を備えたVCSSLプログラムを開発するためのガイドです。
2DCG開発ガイド
画面上や画像ファイルなどに、2次元的な描画を行うVCSSLプログラムを開発するためのガイドです。
3DCG開発ガイド
画面上や画像ファイルなどに、3次元的な描画を行うVCSSLプログラムを開発するためのガイドです。
標準ライブラリ 仕様書
コード内で呼び出される関数は、大半が標準ライブラリのものです。その詳細仕様を掲載しています。

ライセンス

このVCSSLコード( 拡張子が「.vcssl」のファイル )は実質的な著作権フリー(パブリックドメイン) である CC0 の状態で公開しています。 そのままでのご利用はもちろん、言語の種類を問わず、改造や流用などもご自由に行ってください。

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

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

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


スポンサーリンク



[ 前へ | 目次 | 次へ ]
波の干渉(面上の円形波)のアニメーション表示

面上の円形波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。
円形波のアニメーション表示

振幅・波長・周期をスライダ―で操作しながら、円形波のグラフをアニメーションで見られるプログラムです。
波の干渉(線上の正弦波)のアニメーション表示

線上(1次元の)の正弦波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。
正弦波のアニメーション表示

振幅・波長・周期をスライダ―で操作しながら、正弦波のグラフをアニメーションで見られるプログラムです。
凹レンズを通過する波のシミュレーション

凹レンズ形状の高密度媒質を通過する、波のシミュレーションです。
凸レンズを通過する波のシミュレーション

凸レンズ形状の高密度媒質を通過する、波のシミュレーションです。
乱雑な密度分布における波のシミュレーション

密度分布が乱雑な媒質中における、波の伝播のシミュレーションです。
ローレンツアトラクタ(ファイル出力版)

4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。
波の屈折のシミュレーション

密度の異なる領域を、波が屈折しながら通過するシミュレーションです。
力学アルゴリズムによる波のシミュレーション(面上の波)

媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。
手動で波を発生させるシミュレーション

スライダーをマウスで動かす事により、波を発生させるシミュレーションです。
力学アルゴリズムによる波のシミュレーション(線上の波)

媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。
二重振り子のシミュレーション

ラグランジュ方程式を用いた、二重振り子のシミュレーションです。
ローレンツアトラクタ(GUI版)

4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。
スポンサーリンク

この階層の目次
[ 前へ | 目次 | 次へ ]
お知らせ

Vnanoのスクリプトエンジンアーキテクチャ解説1: 全体像
2019年05月28日 - RINEARNでは現在、 アプリケーション組み込み用スクリプトエンジン「 Vnano 」を、オープンソースで開発中です。今回は、このスクリプトエンジンのアーキテクチャ面を掘り下げて解説します。

リニアンプロセッサー nano の先行開発版やソースコードリポジトリを公開
2019年04月16日 - オープンソースで開発中の小型プログラム関数電卓ソフト、「 リニアンプロセッサー nano 」の先行開発版やソースコードリポジトリを公開しました。概要と使用方法、ビルド方法などについて解説します。

各ソフトウェアの最新版を一括でリリース、OpenJDKのJava実行環境(JRE)に対応
2019年03月06日 - RINEARNでは3月6日に、主要なソフトウェアの最新版を一括でリリースしました。今回のアップデートには、以前お知らせした、OpenJDKで生成したJREへの対応が含まれています。その概要等をお知らせします。

新着
連番ファイルから3Dグラフをアニメーション描画するツール

フォルダ内の連番データファイルを読み込み、3Dグラフを高速で連続描画して、アニメーションさせるツールです。
2019年06月03日
連番ファイルから2Dグラフをアニメーション描画するツール

フォルダ内の連番データファイルを読み込み、2Dグラフを高速で連続描画して、アニメーションさせるツールです。
2019年05月24日
[公式ガイドサンプル] 立体モデルを生成して3D空間に配置する

「VCSSL 3DCG開発ガイド」内のサンプルコードです。立体モデルを生成し、3D空間に配置します。
2019年05月21日
[公式ガイドサンプル] ポリゴンを生成して3D空間に配置する

「VCSSL 3DCG開発ガイド」内のサンプルコードです。立体の基本的な構成要素となるポリゴンを生成し、3D空間に配置します。
2019年05月20日
[公式ガイドサンプル] CSVファイルにデータを書き出し&読み込んで、複雑な3次元曲面のグラフを描く(魔法陣形)

「VCSSLスタートアップガイド」内のサンプルコードです。CSVファイルにデータを書き出し、さらにそれを読み込んで、魔法陣のような3次元曲面のグラフを描画します。
2019年05月17日
開発元Twitterアカウント

スポンサーリンク