C言語でカオスアトラクタを計算する方法を解説

本稿ではC言語を使ってカオスアトラクタの計算方法について解説します。

カオスアトラクタは、初期条件の微かな違いが結果に大きな影響を与える非線形システムの一例です。

簡単なコード例を通してループ処理や数値計算法を利用し、複雑な挙動をシミュレーションする方法を学べます。

カオスアトラクタの基礎知識

カオスアトラクタの概念

カオスアトラクタは、非線形ダイナミクスにおいて決定論的な法則のもとで生成される複雑な軌道パターンのことです。

システムがカオス状態にある場合、軌道は初期状態から遠く離れた後でも、ある一定の範囲内で収束したような挙動を示します。

すなわち、決定論的な方程式にもかかわらず、長期的な挙動は予測不可能になる性質を持ちます。

また、カオスアトラクタは空間内の幾何学的な構造を持っており、フラクタル次元を備えている場合が多いです。

初期条件が及ぼす影響

カオス系では初期条件に非常に敏感なため、微小な初期値の違いが時間経過とともに大きな違いをもたらすことが知られています。

これはいわゆるバタフライ効果としても言及される現象です。

初期条件の差分$\mathrm{\Delta }x(0)$が、時間$t$が経つと、以下の式に従い急速に拡大する可能性があります。

$$\mathrm{\Delta }x(t)\approx e^{\lambda t}\mathrm{\Delta }x(0)$$

ここで、$\lambda$は最大リャプノフ指数と呼ばれるパラメータで、正の値であればカオス的な振る舞いを示します。

初期条件の僅かな違いが軌道全体に大きな影響を与えるため、数値シミュレーションでは入力値の正確な設定が重要となります。

数値計算モデルの選定

数値積分法の基本

カオスアトラクタの計算では、連続時間の微分方程式を離散化して数値的に解く方法がよく利用されます。

一般的には、以下のような数値積分法が採用されます。

• オイラー法

$$x_{n+1}=x_n+f(x_n)\mathrm{\Delta }t$$

• 4次のルンゲクッタ法

より高精度に解を求めるために、ルンゲクッタ法がよく利用されます。

数値積分法はシステムの安定性や精度に影響するため、具体的なシステムに合わせたアルゴリズムの選定が大切です。

適用可能なシステム例

ローレンツシステムの計算方法

ローレンツシステムはカオスの代表例で、以下の微分方程式で表現されます。

$$\frac{dx}{dt}=\sigma (y–x),\frac{dy}{dt}=x(\rho –z)-y,\frac{dz}{dt}=xy–\beta z$$

ここでは、オイラー法を用いたシンプルな数値シミュレーションのサンプルコードを示します。

#include <stdio.h>
// パラメータの定義
#define SIGMA 10.0
#define RHO 28.0
#define BETA 8.0/3.0
#define DT 0.01   // 時間刻み
#define STEPS 10000  // ループ回数
int main(void) {
    double x = 1.0, y = 1.0, z = 1.0; // 初期条件
    int i;
    // シミュレーション開始
    for(i = 0; i < STEPS; i++) {
        // 現在の値を保持
        double x_now = x, y_now = y, z_now = z;
        // オイラー法による更新
        x = x_now + DT * SIGMA * (y_now - x_now);
        y = y_now + DT * (x_now * (RHO - z_now) - y_now);
        z = z_now + DT * (x_now * y_now - BETA * z_now);
        // 結果を出力
        printf("%.4f %.4f %.4f\n", x, y, z);
    }
    return 0;
}

1.0000 1.0000 1.0000
1.0000 1.0000 1.0000
...

ロジスティック写像による例

ロジスティック写像もカオスの特徴を示すシステムの一例です。

写像は以下の式で表されます。

$$x_{n+1}=r{x}_n(1–x_n)$$

ここで、パラメータ$r$の値がある範囲に入るとカオス的な挙動が観察されます。

サンプルコードでは、オイラー法ではなく、直接写像を適用して反復計算を行います。

#include <stdio.h>
#define R 3.9
#define ITERATIONS 50
int main(void) {
    double x = 0.5;  // 初期条件
    int i;
    // シミュレーション開始
    for(i = 0; i < ITERATIONS; i++) {
        // ロジスティック写像の計算
        x = R * x * (1 - x);
        // 結果を出力
        printf("%.4f\n", x);
    }
    return 0;
}

0.9750
0.0956
0.3380
...

C言語での実装方法

プログラム構造と変数設計

C言語による実装では、プログラム全体の構造を明確にして、変数の命名規則やスコープを意識した設計を行うことが大切です。

基本構造としては、main関数の中に初期条件の設定、シミュレーションループ、結果出力を含む形を採用します。

たとえば、システムの状態を格納するために構造体を利用する方法も考えられます。

#include <stdio.h>
// 状態を格納するための構造体定義
typedef struct {
    double x;
    double y;
    double z;
} State;
int main(void) {
    // 初期状態の設定
    State state = {1.0, 1.0, 1.0};
    // 状態更新に必要な変数宣言
    double dt = 0.01;
    int steps = 10000;
    int i;
    // シミュレーションループ
    for(i = 0; i < steps; i++) {
        // 変数の更新処理は後述のループ処理で行ないます
        // ここでは状態出力も行います
        printf("%.4f %.4f %.4f\n", state.x, state.y, state.z);
    }
    return 0;
}

1.0000 1.0000 1.0000
1.0000 1.0000 1.0000
...

ループ処理と数値更新手法

数値計算を行うためのループ処理は、forループやwhileループを用いて実装します。

各ステップで変数の値を更新する際には、計算誤差や丸め誤差に注意しながら、必要に応じて一時変数を利用して値が上書きされないように工夫します。

以下はオイラー法によるシンプルな更新処理の例です。

#include <stdio.h>
#define DT 0.01
#define STEPS 10000
#define PARAM 2.0  // 任意のパラメータ例
int main(void) {
    double x = 1.0;  // 初期条件
    int i;
    for(i = 0; i < STEPS; i++) {
        double x_current = x;
        // 例として簡単な更新式: x = x + dt * PARAM * x
        x = x_current + DT * PARAM * x_current;
        // 更新結果の出力
        printf("%.4f\n", x);
    }
    return 0;
}

1.0200
1.0404
1.0612
...

出力処理の実装

シミュレーション結果は、printf関数などを利用してターミナルに出力されます。

大量のデータを出力する場合には、フォーマット指定子を活用して必要な桁数に調整することで、視認性が向上します。

出力結果は後で検証やグラフ描画などに利用するため、ファイル出力に変更する方法も検討できます。

#include <stdio.h>
#define STEPS 100
#define DT 0.01
int main(void) {
    double x = 1.0;
    int i;
    // ループ処理と出力処理
    for(i = 0; i < STEPS; i++) {
        x = x + DT * (2.0 * x);  // 任意の更新式
        // 結果をコンソールに出力
        printf("Step %d: x = %.4f\n", i, x);
    }
    return 0;
}

Step 0: x = 1.0200
Step 1: x = 1.0404
Step 2: x = 1.0612
...

検証とデバッグ

実行結果の確認方法

シミュレーション実行後は、出力された数値データをグラフ描画ツール(たとえばGnuplotなど)でプロットし、軌道の挙動を視覚的に確認する方法があります。

また、既知のカオス系の特性値やフラクタル次元と比較することで、数値シミュレーションの正確性を検証することが可能です。

シンプルなログ出力を活用し、ステップごとの数値の変化を確認できるようにすることも有効です。

数値不安定性への対処法

カオス系のシミュレーションでは、数値不安定性が発生する可能性があるため、以下の対策を検討するとよいです。

• 時間刻み$\mathrm{\Delta }t$を適切に設定し、小さすぎず大きすぎない値を選ぶ
• ルンゲクッタ法など、より高精度な数値積分法を利用する
• シミュレーション開始前に初期条件やパラメータの妥当性を確認する

数値シミュレーション中に不安定な振る舞いが検出された場合、デバッグ用のログ出力を追加して、どのステップで発散が始まったかなどを特定することが効果的です。

必要に応じて、シミュレーションパラメータを動的に調整する仕組みを取り入れるのも一つの方法です。

まとめ

本記事ではカオスアトラクタの基本的な概念と、初期条件が軌道に及ぼす影響について学びました。

また、数値積分法を基盤としたローレンツシステムやロジスティック写像を例に、カオス現象を数値的に表現する方法を紹介しました。

さらに、C言語による実装方法、ループ処理と数値更新、出力方法、そして実行結果の検証や数値不安定性の対処法についても詳述しています。