1

いわゆるLorenz95 モデル(1995 年に Ed Lorenz によって発明された) の時間ステップをインクリメントするコードがあります。通常は 40 変数のモデルとして実装され、無秩序な動作を示します。アルゴリズムの時間ステップを次のようにコーディングしました。

class Lorenz:
    '''Lorenz-95 equation'''

    global F, dt, SIZE
    F = 8
    dt = 0.01
    SIZE = 40

    def __init__(self):
        self.x = [random.random() for i in range(SIZE)]

    def euler(self):
        '''Euler time stepping'''
        newvals = [0]*SIZE
        for i in range(SIZE-1):
            newvals[i] = self.x[i] + dt * (self.x[i-1] * (self.x[i+1] - self.x[i-2]) - self.x[i] + F)
        newvals[SIZE-1] = self.x[SIZE-1] + dt * (self.x[SIZE-2] * (self.x[0] - self.x[SIZE-3]) - self.x[SIZE-1] + F)
        self.x = newvals

この関数 euler は遅くはありませんが、残念ながら、私のコードでは非常に多くの呼び出しを行う必要があります。タイムステッピングをコーディングして高速に実行する方法はありますか?

どうもありがとう。

4

1 に答える 1

3

少なくとも 2 種類の最適化が考えられます。よりスマートな方法 (アルゴリズムの改善) で作業することと、より高速に作業することです。

  • アルゴリズム側では、オイラー法を使用しています。これは一次法であり (したがって、グローバル エラーはステップ サイズに比例します)、小さな安定領域があります。つまり、あまり効率的ではありません。

  • 反対に、標準の CPython 実装を使用している場合、この種のコードは非常に遅くなります。これを回避するには、単純にPyPyで実行してみてください。そのジャストインタイム コンパイラは、数値コードをおそらく 100 倍高速に実行できます。カスタム C または Cython 拡張機能を作成することもできます。

しかし、もっと良い方法があります。常微分方程式系を解くことは非常に一般的であるため、Python のコア科学ライブラリの 1 つである scipy は、それらを解くために、高速で実戦テスト済みの Fortran ライブラリをラップしています。scipy を使用することで、アルゴリズムの改善 (インテグレーターはより高い次数を持つため) と高速な実装の両方を得ることができます。

一連の摂動のある初期条件について Lorenz 95 モデルを解くと、次のようになります。

import numpy as np


def lorenz95(x, t):
    return np.roll(x, 1) * (np.roll(x, -1) - np.roll(x, 2)) - x + F

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from scipy.integrate import odeint
    SIZE = 40
    F = 8
    t = np.linspace(0, 10, 1001)
    x0 = np.random.random(SIZE)
    for perturbation in 0.1 * np.random.randn(5):
        x0i = x0.copy()
        x0i[0] += perturbation
        x = odeint(lorenz95, x0i, t)
        plt.plot(t, x[:, 0])
    plt.show()

そして、出力（設定np.random.seed(7)、あなたのものは異なる場合があります）はかなり混沌としています。初期条件 (座標の 1 つだけで!) での小さな摂動は、非常に異なる解を生成します。 Lorenz-95 力学系

しかし、それは本当にオイラー時間ステッピングより速いのでしょうか? ほぼ 3 倍高速にdt = 0.01見えますが、ソリューションは最初の部分を除いて一致しません。 オイラー vs オイント

を減らすとdt、オイラー法によって得られる解はますますodeint解に近づきますが、時間がかかります。dt が小さいほど、後のオイラー解が解を見失うことに注意してくださいodeint。最も正確なオイラー解は、t=10 までの odeint よりも t=6 までの解を計算するのに 600 倍の時間がかかりました。完全なスクリプトはこちらをご覧ください。 オイラー vs オイント

結局、このシステムは非常に不安定であるため、odeint の解でさえ、プロットされたすべての時間で正確ではないと思います。

于 2013-05-08T17:23:40.793 に答える