データ同化｜ローレンツ方程式とPythonによる実装

今回はデータ同化領域の技術評価でよく使われるローレンツ方程式を取り上げます🐜

このローレンツ方程式は、カオス研究¹が始まるきっかけとなったものです。

ローレンツアトラクタの例（Wikipediaより）

本記事では上記のようなグラフを作ることを目的とします。

ローレンツ方程式
ローレンツ63モデル
- ルンゲ＝クッタ法
サンプルプログラム
モデルの振る舞い
おわりに
参考文献

ローレンツ方程式

ローレンツ方程式とは、気象モデルを研究していたEdward N. Lorenzが、「Deterministic Nonperiodic Flow」（1963年）という論文で発表した非線形の常微分方程式で、カオス的なふるまいを示すモデルです。

カオス的なふるまいとは、あるモデルの初期値をほんの少し変えただけで、計算される結果が大きく変わることを示します。

ローレンツ方程式（またはローレンツモデル）には、1963年発表の63モデルと1996年発表の96モデル（学会セミナで発表？）がありますが、今回は次元数の低い63モデルを扱います。

ローレンツ63モデル

ローレンツ63モデルを以下に示します。

$\displaystyle \begin{align} \frac{dx}{dt} &= - \sigma (x - y) \\ \frac{dy}{dt} &= -xz + rx - y \\ \frac{dz}{dt} &= xy - bz \end{align}$

モデル自体は単純な形ですね。ちなみに差分方程式で表すと、以下のように書けます。

$\displaystyle \begin{align} x_{t} &= x_{t-1} + \Delta t (- \sigma x_{t-1} + \sigma y_{t-1}) \\ y_{t} &= y_{t-1} + \Delta t (-x_{t-1}z_{t-1} + rx_{t-1} - y_{t-1}) \\ z_{t} &= z_{t-1} + \Delta t (x_{t-1}y_{t-1} - bz_{t-1}) \end{align}$

ルンゲ＝クッタ法

ローレンツ63モデルのふるまいを求めることは、いわゆる数値解析で微分方程式の近似解を求めることです。

できるだけ真値に近しい値を求めるためには、時刻の刻み幅 $\Delta t$ を細かくとる方法がすぐに思いつきます。
しかし、それだけでは時刻が進むにつれて誤差が積算していき、不十分です。

より精度の高い近似解を得る手法として、オイラー(Euler)法などがありますが、今回はよく使われるルンゲ＝クッタ法(Runge-Kutta method)の4次(RK4)を扱います。RK4は実装が簡単で十分な精度が得られるため、よく使われる手法です。

ルンゲ＝クッタ法の詳細については、以下のWebページを参考にしてください²。

サンプルプログラム

Pythonでサンプルプログラムを作りました。環境さえ整っていれば、コピペで動きます。

開発環境

開発環境は以下の通りです。

項目	内容
OS	Windows 10
Python	3.6.8
NumPy	1.18.5
Pandas	1.0.4
SciPy	1.4.1
Matplotlib	3.2.1

Lorenz63.py

スクリプトは以下の通りです。動かすと、計算結果が.csvファイルで保存されて、グラフが3枚出てきます。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Lorenz63シミュレーション

Lorenz63シミュレーション

"""
__author__  = 'ari23(Twitter: @ari23ant)'
__version__ = '0.0.1'
__date__    = '2020/06/11'
__status__  = 'Development'

import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


class Lorenz63():

    def __init__(self):
        self.msg = 'Lorenz63 simulation'
        # --- solve_ivp引数 --- #
        # 時刻tの間隔 t0からtfまで
        self.t0 = 0
        self.tf = 100
        self.t_span = (self.t0, self.tf)
        # 関数(Lorenz63)の初期値
        self.init_x = 1
        self.init_y = 1
        self.init_z = 1
        self.y0 = np.array([self.init_x, self.init_y, self.init_z])
        self.y0_ = np.array([self.init_x, self.init_y, self.init_z+0.001])
        # 数値計算手法
        self.method = 'RK45'  # 4次のルンゲクッタ 精度は5次相当
        # 連続値フラグ
        self.dense_output = True  # Falseだと最終値しか戻らない
        # 関数の初期値以外の設定値(Lorenz63パラメタ)
        self.s = 10
        self.r = 28
        self.b = 8/3
        self.args = (self.s, self.r, self.b)
        # 時刻tの刻み幅
        self.dt = 0.01
        self.t = np.arange(self.t0, self.tf, self.dt)

        # --- 出力 --- #
        self.fname = 'lorenz63.csv'

    def Process(self):
        print(self.msg)
        # --- 常微分方程式計算 --- #
        # scipy.integrate.solve_ivpでLorenz63モデルを計算
        self.solver = solve_ivp(fun=self.lorenz63, t_span=self.t_span, y0=self.y0, method='RK45',
                                t_eval=self.t, dense_output=self.dense_output, args=self.args)

        # z初期値ちょいずれ
        self.solver_ = solve_ivp(fun=self.lorenz63, t_span=self.t_span, y0=self.y0_, method='RK45',
                                t_eval=self.t, dense_output=self.dense_output, args=self.args)

        # 変数入れ替え
        self.xyz = self.solver.y
        self.xyz_ = self.solver_.y

        # --- csv保存 --- #
        self.arr_txyz = np.vstack((self.t, self.xyz))
        self.df_txyz = pd.DataFrame(self.arr_txyz.T, columns=['t', 'x', 'y', 'z'])  # 行と列を入れ替え
        self.df_txyz.to_csv(self.fname, index=False)

        # --- グラフ --- #
        # 3次元
        fig3d = plt.figure()
        ax3d = fig3d.add_subplot(111, projection='3d')
        ax3d.plot(self.xyz[0], self.xyz[1], self.xyz[2])
        ax3d.set_xlabel('x')
        ax3d.set_ylabel('y')
        ax3d.set_zlabel('z')
        fig3d.tight_layout()
        # 2次元
        fig2d = plt.figure()
        axx = fig2d.add_subplot(311)
        axy = fig2d.add_subplot(312)
        axz = fig2d.add_subplot(313)
        axx.plot(self.t, self.xyz[0])
        axx.set_ylabel('x')
        axy.plot(self.t, self.xyz[1])
        axy.set_ylabel('y')
        axz.plot(self.t, self.xyz[2])
        axz.set_ylabel('z')
        axz.set_xlabel('t')
        fig2d.tight_layout()
        # 比較
        fig2d_ = plt.figure()
        axx_ = fig2d_.add_subplot(311)
        axy_ = fig2d_.add_subplot(312)
        axz_ = fig2d_.add_subplot(313)
        axx_.plot(self.t, self.xyz[0]-self.xyz_[0])
        axx_.set_ylabel('x')
        axy_.plot(self.t, self.xyz[1]-self.xyz_[1])
        axy_.set_ylabel('y')
        axz_.plot(self.t, self.xyz[2]-self.xyz_[2])
        axz_.set_ylabel('z')
        axz_.set_xlabel('t')
        fig2d_.tight_layout()
        # 表示
        plt.show()

    def lorenz63(self, t, arr_xyz, ss, rr, bb):
        # 時刻t-1の値
        x, y, z = arr_xyz
        # Lorenz63モデルのパラメタ
        s, r, b = ss, rr, bb

        # 時刻tの値
        dxdt = -s * (x-y)
        dydt = -x*z + r*x -y
        dzdt = x*y - b*z

        # arrayで戻す
        return np.array([dxdt, dydt, dzdt])


if __name__ == '__main__':
    proc = Lorenz63()
    proc.Process()

解説

丁寧にコメントを書いたので、以下だけ解説します。

# --- 常微分方程式計算 --- #
# scipy.integrate.solve_ivpでLorenz63モデルを計算
self.solver = solve_ivp(fun=self.lorenz63, t_span=self.t_span, y0=self.y0, method='RK45',
                        t_eval=self.t, dense_output=self.dense_output, args=self.args)

Pythonでルンゲ＝クッタ法を実装すると、計算時間が膨大になることは目に見えていた³ので、今回はSciPyの常微分方程式を解く関数scipy.integrate.solve_ivpを使いました。

scipy.integrate.solve_ivpを使う上で注意する点が3つあります⁴。

常微分方程式の第1引数に必ずtを置く
今回実装したLorenz63はtを必要としませんが、書いておかないと動きません。
常微分方程式のパラメタは1変数ずつ用意する
Lorenz63ではsとrとbの3つのパラメタがありますが、これをリストなどでひとまとめにすると動きません。
連続値が欲しいときはdense_outputにTrueを入れる
これをしないと最後の時刻の結果のみしか出力されません。

その他の使い方については、公式のドキュメントを参考にしてください。