» 詳しい使用方法や、エラーで展開できない際の対応方法などはこちら
円形波のアニメーション表示
振幅・波長・周期などのパラメータに基づいて、面の上を伝わる円形波をアニメーションで描画するVCSSLプログラムです。
前回 と 前々回 は、線の上、つまり1次元的な媒質上を伝わる 正弦波 のアニメーションを扱いました。 今回はその延長で、面の上、つまり2次元的な媒質上を伝わる波を扱ってみましょう。 具体的には、波源から同心円状に広がる円形波をアニメーションで表示してみます。
[ 関連記事 ]
- 前回: 波の干渉(線上の正弦波)のアニメーション表示
- 前々回: 正弦波のアニメーション表示
スポンサーリンク
使用方法
ダウンロードと展開(解凍)
まず、PC(スマホは未対応)で上の画面の「 ダウンロード 」ボタンを押してください。 するとZIP形式で圧縮されたファイルがダウンロードされます。
その後、ZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開(解凍)してください。 展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。
» 展開がエラーで止まってしまう場合や、ファイル名が文字化けしてしまう場合は…
プログラムの起動
Windows をご使用の場合
上記でZIPファイルを展開したフォルダ内にある、以下のバッチファイルをダブルクリック実行してください:
もしプログラムを書き変えながら使いたい場合は、代わりに「 VCSSL_Editor__プログラム編集はこちら.bat 」を実行してください。
正常に起動できると、初回のみ、Java実行環境を入手するか等を尋ねられるので、適時答えて済ませると、プログラムが起動します。 2回目以降はすぐに起動します。
Linux 等をご使用の場合
ZIPファイルを展開したフォルダ内へコマンドライン端末で cd して、以下の通り入力して実行してください:
(プログラムの内容を書き変えながら使いたい場合は、代わりに VCSSL_Editor.jar を実行)
» javaコマンドが使用できない等のエラーが表示される場合は…
起動後の画面
起動すると、まず波のパラメータを数値入力するか尋ねられます。 特に「ぴったりこの値の波長にしたい」などの必要が無ければ、「いいえ」でスキップして構いません。
その後、グラフ画面とパラメータ設定画面が立ち上がります。

グラフ画面では、円形波がアニメーション表示されます。 パラメータ設定画面では、波の波長や振幅、周期などの値をスライダーで操作できます。
ウィンドウを閉じると、プログラムの実行が終了します。
題材解説
円形波は、波源と呼ばれる振動源の動きが、周囲の媒質へ同心円状に伝わっていく事で起こります。

お風呂などで、水面をそっと指で触れると、実際にそこから円形に波が広がっていくのが見られるかと思います。 さらに、指を水面に対して垂直に振動させると、同心円状の波が連続的に広がっていきます(この場合は指が波源です)。
この運動を単純化して表すために、まず「波の各点は、波源と全く同じ動きを、波源からの距離に応じて遅れながら行っている(スポーツ応援のウェーブのように)と仮定」しましょう。後述の通り厳密にはちょっとおかしい仮定なのですが(エネルギー保存則的に)、しかしこうすると一気に話が簡単になります。
そしてこの条件下で、波源が 単振動(調和振動) の動きをするとしましょう。するとその動きが、距離が離れた領域に、その距離分だけ遅れて伝わっていきます。 結果、波を放射方向に切ると、その断面はまさしく前々回に扱った正弦波となります。

ここで波の振動方向(=波源が上下に振動する方向)を \(z\) 軸にとり、波源からの水平距離を \(r\) としています。
この運動を数式で表すために、波源の位置を座標原点にとると、面上のある点 \((x,y)\) における波源からの水平距離 \(r\) は以下のように表されます :
\[ r = \sqrt{x^2 + y^2} \tag{1}\]図にすると以下の通りです:

あとは、前々回に登場した、正弦波の式の横軸 \(x\) を \(r\) で置き直せば、「波源から放射状に伝わる正弦波」の式が求まります:
\[ z = A \sin 2 \pi \bigg( \frac{t}{T} - \frac{r}{\lambda} \bigg) \tag{2} \]これが今回アニメーションで表現している円形波の数式表現です。理論上扱いやすいシンプルな円形波で、波の干渉の説明などの際にも、よくこの形の円形波が用いられます。
なお、\( A \) は振幅、\( \lambda \) は波長、\( T \) は周期で、これも 前々回 と同じなので、詳しくはそちらの説明をご参照ください。
ところで、波は広がるのに伴って減衰(振幅が小さくなっていく)したり、壁で波が反射したりと、実際には複雑な現象が色々とあります。 上の (2) 式は、そういうものが一切無視できるという、かなり割り切った仮定をしている事に留意が必要です。
実際に面上の波の動きを、もう少し物理的にシミュレーションしたものは、以下の回で扱っています:
上のシミュレーションは、かなり色々な事を簡略化したシミュレーションモデルを用いていますが、 それでも開始と同時に振幅が急に減少し、さらに反射波が重なり合ってきて、すぐによくわからない動きになってしまっています。
実際の波についても、反射の影響については工夫でどうにかできるにしても、広い領域に広がっていく波が減衰していく事は、実はエネルギー保存則と表裏一体です。 なので、特に波源近くに注目する場合は、 恐らく(相当変な事をしない限り)現実的に減衰を無視できるような状況を作り出す事は難しいと思います。 つまるところ、普通の媒質の上で広がっていく波について、「減衰は無視できるものとする」という仮定を置く事は、 間接的に「エネルギー保存則は無視できるものとする」という仮定を置いているようなものなわけで、 そういう仮定を置く事自体がどうなのかという議論にも当然なり得ます。
一方で、(2) 式のように、(ある意味で非現実的な)単純に距離に関する正弦波で表される円形波を考える事で、 次回で扱う「円形波の干渉」で、独特のパターンが表れる理由を、かなり単純化して理解しやすくなります。 減衰を入れるとそれが一気に難しくなるため、最初から減衰ありで話を進めると、 「 重ね合わせたら結果的にこうなるからこうなる 」といった曖昧な説明になってしまいます。 物理的には減衰がある方が本質的に正しいとしても、 まずは減衰の無い場合について考えて、それから減衰があるとどうなるかを考える、 という方が、思考の流れ的にはワンクッション置けて、わかりやすいのではないでしょうか。 その点で、(2) 式のように割り切った円形波を扱う事にも、一定の利点があると思います。 もちろん、物理的に結構無理な仮定を置いている事に留意しておく必要はあります。
スポンサーリンク
コード解説
このプログラムのコードはVCSSLで記述されています。 ここではそのコード内容について簡単に解説します。 VCSSLはC系の単純な文法の言語なので、C言語などに触れた事のある方なら簡単に読めると思います。
グラフの範囲などを変えたり、色々と改造したいといった場合などは、 プログラムのコード「 CircularWave.vcssl 」をテキストエディタで開いて改造してください。 スクリプト言語なので、コンパイラなどの別ソフトは不要で、コードを書き換えるだけでOKです。
コード全体
まずは、コード全体を見てみましょう。
coding UTF-8; // 文字コードの明示(文字化け予防)
import Math; // 数学関数の使用に必要なライブラリの読み込み
import GUI; // GUI画面の構築に必要なライブラリの読み込み
import tool.Graph3D; // 3次元グラフ描画に必要なライブラリの読み込み
// 以下、波のパラメータを格納する変数
double amplitude = 1.0; // 波の振幅(A)
double wavelength = 2.0; // 波の波長(λ)
double period = 1.0; // 波の周期(T)
// グラフ範囲や点数、時間刻みなどの設定変数
const int N = 100; // メッシュの一方向あたりの格子点数(細かいほど滑らか)
const double X_MIN = -5.0; // グラフのX軸の最小値
const double X_MAX = 5.0; // グラフのX軸の最大値
const double Y_MIN = -5.0; // グラフのY軸の最小値
const double Y_MAX = 5.0; // グラフのY軸の最大値
const double Z_MIN = -1.0; // グラフのZ軸の最小値
const double Z_MAX = 1.0; // グラフのZ軸の最大値
const double DT = 0.1; // 1ループごとの時間進行量(speed=1.0の時の値)
const double DX = (X_MAX - X_MIN) / (N-1); // メッシュのX軸方向の刻み幅
const double DY = (Y_MAX - Y_MIN) / (N-1); // メッシュのY軸方向の刻み幅
// グラフに波のデータを渡すための座標値配列
double waveX[ N ][ N ]; // X座標の配列
double waveY[ N ][ N ]; // Y座標の配列
double waveZ[ N ][ N ]; // Z座標の配列
// GUI(設定画面)用の変数
int window; // ウィンドウのIDを格納する
int periodSlider; // 周期設定スライダーのIDを格納する
int wavelengthSlider; // 波長設定スライダーのIDを格納する
int amplitudeSlider; // 振幅設定スライダーのIDを格納する
int speedSlider; // アニメーション速度スライダーのIDを格納する
int infoLabel; // 設定値を表示するラベルのIDを格納する
int timeLabel; // 時刻表示ラベルのIDを格納する
int resetButton; // 時刻リセットボタンのIDを格納する
// その他の制御用変数
float speed = 0.7; // アニメーション速度
int graph; // グラフソフトのIDを格納する
bool continuesLoop = true; // アニメーションON/OFF制御変数
double t; // 時刻変数
// ==================================================
// 中心的な処理 - この関数は起動時に自動で実行されます
// ==================================================
void main( ){
// グラフを起動
createGraphWindow();
// GUI(設定画面)を構築
createSettingWindow();
t = 0.0; // 時刻変数
// アニメーションのループ
while( continuesLoop ) {
// グラフソフトに渡す点列データの座標を計算
for( int yi=0; yi<N; yi++ ) {
for( int xi=0; xi<N; xi++ ) {
double x = X_MIN + xi * DX; // X軸方向のxi番目のメッシュ線のx座標値
double y = Y_MIN + yi * DY; // Y軸方向のyi番目のメッシュ線のy座標値
double r = sqrt(x*x + y*y); // 原点からの距離
waveX[yi][xi] = x;
waveY[yi][xi] = y;
waveZ[yi][xi] = amplitude * sin( 2.0 * PI * ( t/period - r/wavelength ) );
}
}
// グラフソフトにデータを転送
setGraph3DData( graph, waveX, waveY, waveZ );
// 時刻を加算(アニメーション速度をかけて変化のスピードを調整)
t += DT * speed;
// アニメーションウェイト(30ミリ秒待つ)
sleep( 30 );
}
exit(); // アニメーションループを脱出すると実行終了
}
// ==================================================
// グラフソフトを起動して設定を行う関数
// ==================================================
void createGraphWindow() {
// 3Dグラフソフトを起動
graph = newGraph3D( 0, 0, 800, 600, "Graph Window" );
// 最初に、偏平率や目盛り調整、角度や倍率などの細かい設定をファイルから読み込む
setGraph3DConfigurationFile(graph, "RinearnGraph3DSetting.ini");
// 範囲の自動調整機能を無効化し、手動で範囲を指定
setGraph3DAutoRange( graph, false, false, false );
setGraph3DRange( graph, X_MIN, X_MAX, Y_MIN, Y_MAX, Z_MIN, Z_MAX );
// プロットオプションを選択(点プロットを無効化、曲面プロットを有効化)
setGraph3DOption(graph, "WITH_POINTS", false);
setGraph3DOption(graph, "WITH_MEMBRANES", true);
}
// ==================================================
// 設定画面を構築する関数
// ==================================================
void createSettingWindow() {
// ウィンドウの生成
window = newWindow( 0, 610, 800, 230, "Setting Window" );
// 時刻表示ラベルを生成して配置
timeLabel = newTextLabel( 10, 10, 200, 20, "" );
mountComponent( timeLabel, window );
// スライダーの左に振幅、周期、波長のラベルを生成して配置
int amplitudeLabel = newTextLabel( 45, 40, 100, 15, "A 振幅 (0〜1)" );
mountComponent( amplitudeLabel, window );
int omegaLabel = newTextLabel( 45, 60, 100, 15, "T 周期 (0.2〜10)" );
mountComponent( omegaLabel, window );
int lambdaLabel = newTextLabel( 45, 80, 100, 15, "λ 波長 (1〜8)" );
mountComponent( lambdaLabel, window );
// 振幅スライダー(範囲0.0〜1.0)を生成して配置
amplitudeSlider = newHorizontalSlider( 150, 40, 600, 20, amplitude, 0.0, 1.0 );
mountComponent( amplitudeSlider, window );
// 周期スライダー(範囲0.2〜10.0)を生成して配置
periodSlider = newHorizontalSlider( 150, 60, 600, 20, period, 0.2, 10.0 );
mountComponent( periodSlider, window );
// 波長スライダー(範囲1.0〜8.0)を生成して配置
wavelengthSlider = newHorizontalSlider( 150, 80, 600, 20, wavelength, 1.0, 8.0 );
mountComponent( wavelengthSlider, window );
// スライダー値を表示するラベルを生成して配置
infoLabel = newTextLabel( 150, 100, 400, 20, "" );
mountComponent( infoLabel, window );
updateInfoLabel();
// アニメーション速度のラベルと調整スライダーを生成して配置
int speedLabel = newTextLabel( 10, 140, 140, 15, "アニメーション速度" );
mountComponent( speedLabel, window );
speedSlider = newHorizontalSlider( 150, 140, 400, 20, speed, 0.0, 1.0 );
mountComponent( speedSlider, window );
// 時刻リセットボタンを生成して配置
resetButton = newButton(600, 125, 150, 40, "時刻をリセット");
mountComponent( resetButton, window );
// 配置後にウィンドウを描画して表示内容を更新
paintComponent( window );
}
// ==================================================
// スライダー値を表示するラベルの内容を更新する関数
// ==================================================
void updateInfoLabel() {
string infoText
= "( A="
+ round(amplitude, 2, HALF_UP) // 振幅を小数点以下2桁に丸める
+ ", λ="
+ round(wavelength, 2, HALF_UP) // 波長を小数点以下2桁に丸める
+ ", T="
+ round(period, 2, HALF_UP) // 周期を小数点以下2桁に丸める
+ ")";
setComponentString(infoLabel, infoText);
paintComponent(infoLabel);
paintComponent(window);
}
// ==================================================
// スライダーが動いた際に呼ばれる関数(イベントハンドラ)
// ==================================================
void onSliderMove(int id, float value) {
// 以下、スライダーに対応するグローバル変数に値を設定
if(id == amplitudeSlider) {
amplitude = value;
}
if(id == periodSlider) {
period = value;
}
if(id == wavelengthSlider) {
wavelength = value;
}
if (id == speedSlider) {
speed = value;
}
updateInfoLabel();
}
// ==================================================
// ボタンを押した際に呼ばれる関数(イベントハンドラ)
// ==================================================
void onButtonClick(int id) {
// 時刻リセットボタンが押されたら時刻変数 t をリセット
if (id == resetButton) {
t = 0.0;
}
}
// ==================================================
// ウィンドウを閉じた際に呼ばれる関数(イベントハンドラ)
// ==================================================
void onWindowClose(int id) {
// アニメーションループを脱出して終了
continuesLoop = false;
}
// ==================================================
// グラフ画面を閉じた際に呼ばれる関数(イベントハンドラ)
// ==================================================
void onGraph3DClose(int id) {
// アニメーションループを脱出して終了
continuesLoop = false;
}
CircularWave.vcssl
以上です。大体250行程度の規模の短いコードですね。 各部で行っている処理はコード内のコメントの通りですが、 以下では先頭から順を追って、もう少し詳しく説明します。
基本的には、前々回の「 正弦波のアニメーション表示 」のコードと非常によく似ていて、 波のデータの次元が増えたのと、表示する2Dグラフが3Dグラフになっただけです。 なので、コードの基本的な流れについては、まず前回のコードをご参照ください。
以下では、前々回から変わっている部分を中心に、ピックアップして解説していきます。
変数の宣言
先頭から続く変数の宣言部では、基本的には前々回と同じですが、 Z方向の次元が新たに増えたので、それに関連する変数が増えています。
// 以下、波のパラメータを格納する変数
double amplitude = 1.0; // 波の振幅(A)
double wavelength = 2.0; // 波の波長(λ)
double period = 1.0; // 波の周期(T)
// グラフ範囲や点数、時間刻みなどの設定変数
const int N = 100; // メッシュの一方向あたりの格子点数(細かいほど滑らか)
const double X_MIN = -5.0; // グラフのX軸の最小値
const double X_MAX = 5.0; // グラフのX軸の最大値
const double Y_MIN = -5.0; // グラフのY軸の最小値
const double Y_MAX = 5.0; // グラフのY軸の最大値
const double Z_MIN = -1.0; // グラフのZ軸の最小値
const double Z_MAX = 1.0; // グラフのZ軸の最大値
const double DT = 0.1; // 1ループごとの時間進行量(speed=1.0の時の値)
const double DX = (X_MAX - X_MIN) / (N-1); // メッシュのX軸方向の刻み幅
const double DY = (Y_MAX - Y_MIN) / (N-1); // メッシュのY軸方向の刻み幅
// グラフに波のデータを渡すための座標値配列
double waveX[ N ][ N ]; // X座標の配列
double waveY[ N ][ N ]; // Y座標の配列
double waveZ[ N ][ N ]; // Z座標の配列
var.txt
waveX / waveY / waveZ などは、グラフに波のデータを渡すために、波が起こる面上を細かく刻んだ各点の座標値を格納する配列です。 前々回では、波は線の上(X軸上)で起こっていたので、それを単純に刻んだ1次元配列でした。 今回は、波は面の上(X-Y平面上)で起こるので、それをメッシュ状に刻んだ2次元配列になっています。
メッシュはX軸方向とY軸方向にそれぞれ N 個の格子点があるように刻み(つまり区間数は N-1 )、その区間幅はそれぞれ DX および DY です。
main 関数
続いて main 関数です。
// ==================================================
// 中心的な処理 - この関数は起動時に自動で実行されます
// ==================================================
void main( ){
// グラフを起動
createGraphWindow();
// GUI(設定画面)を構築
createSettingWindow();
t = 0.0; // 時刻変数
// アニメーションのループ
while( continuesLoop ) {
// 時刻を加算(アニメーション速度をかけて変化のスピードを調整)
t += DT * speed;
// グラフソフトに渡す点列データの座標を計算
for( int yi=0; yi<N; yi++ ) {
for( int xi=0; xi<N; xi++ ) {
double x = X_MIN + xi * DX; // X軸方向のxi番目のメッシュ線のx座標値
double y = Y_MIN + yi * DY; // Y軸方向のyi番目のメッシュ線のy座標値
double r = sqrt(x*x + y*y); // 原点からの距離
waveX[yi][xi] = x;
waveY[yi][xi] = y;
waveZ[yi][xi] = amplitude * sin( 2.0 * PI * ( t/period - r/wavelength ) );
}
}
// グラフソフトにデータを転送
setGraph3DData( graph, waveX, waveY, waveZ );
// アニメーションウェイト(30ミリ秒待つ)
sleep( 30 );
}
exit(); // アニメーションループを脱出すると実行終了
}
main.txt
ここも基本的には前々回と同じで、 while( continuesLoop ) { 〜 のアニメーションループ内で、 時刻を小刻みに進めながら、波の起こる面上を細かく刻んだメッシュの格子点上で、波の値を求めて配列に格納し、 それをグラフに渡して描画してもらっています。
重要なのが、波の値を求めて配列に格納している部分:
double x = X_MIN + xi * DX; // X軸方向のxi番目のメッシュ線のx座標値
double y = Y_MIN + yi * DY; // Y軸方向のyi番目のメッシュ線のy座標値
double r = sqrt(x*x + y*y); // 原点からの距離
waveX[yi][xi] = x;
waveY[yi][xi] = y;
waveZ[yi][xi] = amplitude * sin( 2.0 * PI * ( t/period - r/wavelength ) );
array.txt
で、この 72行目と76行目は、このページの冒頭で登場した (1) - (2) 式:
\[ r = \sqrt{x^2 + y^2} \tag{1} \] \[ z = A \sin 2 \pi \bigg( \frac{t}{T} - \frac{r}{\lambda} \bigg) \tag{2} \]をそのままプログラムのコードの形にしたものです。 つまり今回の処理は、パラパラ漫画の要領で時刻 t を小刻みに進めながら、(2) 式を3Dグラフに連続で描かせてアニメーションしている内容になっています。
コード内で本質的に重要なのはここまでで、残りの部分はグラフ周りの設定や、GUI画面の構築などの枝葉の処理です。 グラフ周りでは、VCSSLでは2Dグラフと3Dグラフを同じ感覚で扱えるようになっているため、 呼び出している関数名の「2D」が「3D」になっているくらいで、やっている事や流れは前々回とほとんど変わりません。 GUI構築やイベント処理に関しても、前々回とほぼ同じです。
コード内容は以上です。
ライセンス
この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 氏の米国およびその他の国における商標または登録商標です。
- その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。
Vnano版 | ローレンツ方程式を数値的に解くプログラム |
|
![]() |
ローレンツ方程式を4次ルンゲ=クッタ法によって解き、グラフ描画用のデータを出力するプログラムです。 |
波の干渉(面上の円形波)のアニメーション表示 |
|
![]() |
面上の円形波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。 |
円形波のアニメーション表示 |
|
![]() |
振幅・波長・周期をスライダ―で操作しながら、円形波のグラフをアニメーションで見られるプログラムです。 |
波の干渉(線上の正弦波)のアニメーション表示 |
|
![]() |
線上(1次元の)の正弦波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。 |
正弦波のアニメーション表示 |
|
![]() |
振幅・波長・周期をスライダ―で操作しながら、正弦波のグラフをアニメーションで見られるプログラムです。 |
凹レンズを通過する波のシミュレーション |
|
![]() |
凹レンズ形状の高密度媒質を通過する、波のシミュレーションです。 |
凸レンズを通過する波のシミュレーション |
|
![]() |
凸レンズ形状の高密度媒質を通過する、波のシミュレーションです。 |
乱雑な密度分布における波のシミュレーション |
|
![]() |
密度分布が乱雑な媒質中における、波の伝播のシミュレーションです。 |
ローレンツアトラクタ(ファイル出力版) |
|
![]() |
4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。 |
波の屈折のシミュレーション |
|
![]() |
密度の異なる領域を、波が屈折しながら通過するシミュレーションです。 |
力学アルゴリズムによる波のシミュレーション(面上の波) |
|
![]() |
媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。 |
手動で波を発生させるシミュレーション |
|
![]() |
スライダーをマウスで動かす事により、波を発生させるシミュレーションです。 |
力学アルゴリズムによる波のシミュレーション(線上の波) |
|
![]() |
媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。 |
二重振り子のシミュレーション |
|
![]() |
ラグランジュ方程式を用いた、二重振り子のシミュレーションです。 |
ローレンツアトラクタ(GUI版) |
|
![]() |
4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。 |