PostgreSQL – 基本コマンド

概要

ターミナルやコマンドプロンプトから実行するPostgreSQLのツールコマンド。

PostgreSQLの開始

psql -U ユーザー名 -d データベース名

または

psql --username ユーザー名 --dbname データベース名

以下はユーザー名:postgres、データベース名:testで開始した例。

  • ユーザーpostgresのパスワードを入力してログイン
  • ログイン後のプロンプトは接続しているDB名となる

データベースのバックアップ

pg_dump -U ユーザー名 データベース名 > バックアップファイル名

以下はユーザー名:postgresでデータベースtestをバックアップファイルtest.dumpに保存する例。

  • ユーザーpostgresのパスワードを入力してログイン
  • バックアップ後は元のプロンプトに戻る
  • 場所を指定しないと、カレントディレクトリーにバックアップファイルが作成される

バックアップファイルのリストア

SQLスクリプトタイプの場合

予めデータベースを作成しておいた上で次のコマンドを実行。

psql -U ユーザー名 -d データベース名 < バックアップファイル名

以下はユーザーpostgresでバックアップファイルtest.dumpからデータベースtestにリストアした例。

カスタムフォーマットの場合

カスタムフォーマットのダンプ(アーカイブファイル)の場合はpg_restoreを使う。

pg_restore -U ユーザー名 -d データベース名 バックアップファイル名

 

Python – 図形を含むアニメーション

ArtistAnimationの場合

ArtistAnimationにグラフを追加するには、たとえばplotメソッドの戻り値をそのままリストに追加していけばよい。

ところが、matplotlib.patchesの図形の場合に同じようにすると、アニメーションにならず全フレームの図形が重ね書きされてしまう。

これを解決するには、appendの際に[]で囲んでリスト化する必要がある。

以下は円の図形のアニメーションのコード例。

実行結果は以下の通り。

グラフと図形の重ね合わせ

グラフと図形を重ね合わせるときも、patchesの図形を明示的にリストとしてappendする必要がある。

以下は、plotによる線とpatches.Circleによる円を重ねて描く例。

実行結果は以下の通り。

FuncAnimationの場合

FuncAnimationの場合は描画関数中でグラフを描くだけで、ArtistAnimationのようなコレクションを用いない。

この場合はArtistAnimationより簡明で、グラフはAxesオブジェクトなどに描画、図形はadd_artistメソッドで追加するだけ。

以下はArtistAnimationと同じく、plotによる線とpatches.Circleによる円を重ね牡蠣する例。

 

Python – FuncAnimationによる動画表示

概要

Pythonのmatplotlibにはanimationパッケージがあり、ArtistAnimationFuncAnimationの2つのクラスが準備されている。

  • ArtistAnimationは、予め動画に表示する各フレームを生成しておいて、これらを動画として表示する
  • FuncAnimationは、全結果を保存せず、フレームごとに描画を実行する関数を実行して、動画して表示する

ここでは、FuncAnimationによる動画の表示方法を整理する(ArtistAnimationについてはこちら)。

1つのグラフの表示

流れ

  1. animationパッケージ又はArtistAnimationクラスをインポート
  2. グラフを描画する関数を定義
  3. FuncAnimationクラスのコンストラクターに描画関数を渡してインスタンスを生成
  4. 表示

処理ステップ

インポート

パッケージをインポートする

または

フレームリストの準備

動画の各フレームを保存するリストを準備する。animationパッケージでは、個々のグラフフレームをartistと表現している。

FigureとAxesオブジェクトを準備

フレームを描画するオブジェクトを準備しておく

1フレームを描画する関数を定義

4周期分の正弦関数を描画する関数を定義している。

  • ここではループごとにxを計算しているが、今回は全フレーム共通なので、ループの外で準備してもよい
  • plot命令のcolor指定をしないと、描画するごとに違うグラフとして色を変えて描画されてしまう

アニメーション設定して表示させる

FuncAnimationクラスのコンストラクターに必要な引数を与えてインスタンスを生成し、表示させる。intervalはフレーム間の待機時間で、ミリ秒で与える。

GIF画像として保存する

結果をGIF画像として保存するには、plt.show()を実行せず、以下のようにimagemagicでファイルとして書き出す。

frames=20でフレーム数を設定しているが、これがないとGIF保存のときにファイル容量が肥大化してしまう。それでも容量は800KB程度になり、ArtistAnimationのときの150KB程度よりも大きい。

実行結果

以下のような動画が表示される。

コード

全体の実行コードは以下のとおり。

複数のグラフの表示

流れ

全体の流れは単一のグラフの場合と同じ。リスト保存時に複数グラフを加算して重ね合わせるところだけが異なる。

  1. animationパッケージ又はArtistAnimationクラスをインポート
  2. グラフを描画する関数を定義し、複数のグラフをプロット
  3. FuncAnimationクラスのコンストラクターに描画関数を渡してインスタンスを生成
  4. 表示

処理ステップ

以下の描画関数が単一グラフの場合と異なる。関数の中で、異なる2つのグラフを描画している

実行結果

以下のように、2つのグラフが同時に表示される。

コード

全体の実行コードは以下のとおり。

描画関数に引数を渡す

描画関数に独自の引数を渡すことができる。

たとえば最初の「1つのグラフの表示」で、frame_countと表示する最大時刻max_xを引数とする場合、まず描画関数でこれらを引数として宣言し、関数内でこれを使った処理を書く。

FuncAnimationインスタンス生成時に、fargsで独自宣言した引数の値をタプルで与える。

 

Python – ArtistAnimationによる動画表示

概要

Pythonのmatplotlibにはanimationパッケージがあり、ArtistAnimationFuncAnimationの2つのクラスが準備されている。

  • ArtistAnimationは、予め動画に表示する各フレームを生成しておいて、これらを動画として表示する
  • FuncAnimationは、全結果を保存せず、フレームごとに描画を実行する関数を実行して、動画して表示する

ここでは、ArtistAnimationによる動画の表示方法を整理する(FuncAnimationについてはこちら)。

1つのグラフの表示

流れ

  1. animationパッケージ又はArtistAnimationクラスをインポート
  2. フレームの数だけグラフを描画してリストに保存
  3. フレームリストを渡してArtistAnimationインスタンスを生成
  4. 表示

処理ステップ

インポート

パッケージをインポートする

または

フレームリストの準備

動画の各フレームを保存するリストを準備する。animationパッケージでは、個々のグラフフレームをartistと表現している。

FigureとAxesオブジェクトを準備

フレームを描画するオブジェクトを準備しておく

フレーム数を指定して各フレームを保存する

以下ではフレーム数を20として、20回分のグラフを少しずつずらして描画して保存している。

ここではグラフを正弦波として、1周期分を20分割して表示させている。

  • ここではループごとにxを計算しているが、今回は全フレーム共通なので、ループの外で準備してもよい
  • plot命令のcolor指定をしないと、描画するごとに違うグラフとして色を変えて描画されてしまう
  • 描画結果はオブジェクトとして、先に準備したartistコレクションに追加している

アニメーション設定して表示させる

ArtistAnimationクラスのコンストラクターに必要な引数を与えてインスタンスを生成し、表示させる。intervalはフレーム間の待機時間で、ミリ秒で与える。

GIF画像として保存する

結果をGIF画像として保存するには、plt.show()を実行せず、以下のようにimagemagicでファイルとして書き出す。

実行結果

以下のような動画が表示される。

コード

全体の実行コードは以下のとおり。

複数のグラフの表示

流れ

全体の流れは単一のグラフの場合と同じ。リスト保存時に複数グラフを加算して重ね合わせるところだけが異なる。

  1. animationパッケージ又はArtistAnimationクラスをインポート
  2. フレームの数だけグラフを描画して、重ね合わせてリストに保存
  3. フレームリストを渡してArtistAnimationインスタンスを生成
  4. 表示

処理ステップ

以下の部分だけが単一グラフの場合と異なる。2つのグラフを描画し、そのオブジェクトを加算して1フレームとした上でartistsに追加している

実行結果

以下のように、2つのグラフが同時に表示される。

コード

全体の実行コードは以下のとおり。

図形を含むアニメーション

matplotlib.patchesによる図形を含んでアニメーション表示するには、少し工夫がいる。その方法はこちら

 

単振動の微分方程式の解

概要

以下の単振動に関する方程式の解法を整理した。

(1)    \begin{equation*} \frac{d^2x(t)}{dt^2} + \omega^2 x(t) = 0 \quad \left( t = 0 \Rightarrow \left\{ \begin{array}{l} x = X_0 \\ \dot{x} = 0 \end{array} \right) \end{equation*}

変数変換による方法

式(1)の両辺にdx/dtを掛けて変形する。

(2)    \begin{align*} &\frac{dx}{dt} \frac{d^2x}{dt^2} + \omega^2 x \frac{dx}{dt} = 0 \\ &\frac{1}{2} \frac{d}{dt} \left( \frac{dx}{dt} \right)^2 + \frac{1}{2} \omega^2 \frac{d}{dt} \left( x^2 \right) = 0 \end{align*}

これをtで積分して以下を得る。

(3)    \begin{equation*} \left( \frac{dx}{dt} \right)^2 + \omega^2 x^2 = C \end{equation*}

ここでC = ω2A2とおいて

(4)    \begin{equation*} \frac{dx}{dt} \right = \pm \omega \sqrt {A^2 - x^2} \quad (x \le A) \end{equation*}

さらにx = A cos θとおいて、

(5)    \begin{equation*} (- A \sin \theta) \frac{d \theta}{dt} = \pm \omega A \sin \theta \quad \Rightarrow \quad \frac{d \theta}{dt} = \mp \omega A \end{equation*}

θの解は以下の様になる。

(6)    \begin{equation*} \theta = \mp \omega t + \alpha \end{equation*}

x = A cos θだったので、以下を得る。

(7)    \begin{equation*} x = A \cos (\omega t + \alpha) \end{equation*}

初期条件から

(8)    \begin{align*} x|_{t=0} &= A \cos (\alpha) = X0 \\ \dot{x}|_{t=0} &= -A \sin (\alpha) = 0 \end{align*}

これより以下の解を得る。

(9)    \begin{equation*} x = X_0 \cos \omega t \end{equation*}

オイラーの公式を使う方法

式(1)の解を以下の様に置く。

(10)    \begin{equation*} x = e^{\lambda t} \end{equation*}

これを式(1)に代入して

(11)    \begin{equation*} \lambda^2 e^{\lambda t} + \omega^2 e^{\lambda t} = 0 \quad \Rightarrow \quad \lambda ^2 + \omega^2 = 0 \quad \Rightarrow \quad \lambda = \pm i \omega \end{equation*}

この解を式(10)に代入して2つの関数が得られるが、式(1)の一般解はこれらの線形結合として得られる。

(12)    \begin{equation*} x = A e^{i \omega t} + Be^{-i \omega t} \end{equation*}

ここで以下のオイラーの公式を適用する。

(13)    \begin{equation*} e^{\pm i \theta} = \cos \theta \pm i \sin \theta \end{equation*}

これにより式(12)は以下のように書ける。

(14)    \begin{align*} x &= C_1 (\cos \omega t + i \sin \omega t) + C_2  (\cos \omega t - i \sin \omega t) \\ &= (C_1 + C_2) \cos \omega t + (C_1 - C_2) i \sin \omega t \end{align*}

xが実数解を持つために、cos, sin双方の係数が実数となるよう、複素係数C1, C2を以下の様に置く。

(15)    \begin{align*} C_1 = \frac{A_1 + i A_2}{2} \\ C_2 = \frac{A_1 - i A_2}{2} \end{align*}

これにより、式(16)は以下の様に変形される。

(16)    \begin{equation*} x = A_1 \cos \omega t + A_2 \sin \omega t \end{equation*}

三角関数の和の関係から、上式は以下の様に表せる。

(17)    \begin{equation*} x = A \cos (\omega t + \alpha) \end{equation*}

初期条件(8)より式(9)と同じ解を得る。

 

 

三角関数 – 2つの三角関数の和

概要

以下の関係を確認する。

(1)    \begin{equation*} A_1 \sin x + A_2 \cos x = A \sin (x + \alpha) = A \cos(x + \beta) \end{equation*}

導出

同じ変数のsineとcosineの和について考える。変数は同じだが、互いの係数が異なっている。

(2)    \begin{equation*} s = A_1 \sin x + A_2 \cos x \end{equation*}

ここで以下のような変数を導入する。

(3)    \begin{equation*} A = \sqrt {A_1^2 + A_2^2} \end{equation*}

これによって式(2)を以下の様に変形する。

(4)    \begin{equation*} s = A \left( \frac{A_1}{A} \sin x + \frac{A_2}{A} \cos x \right) \end{equation*}

ここで、以下の様に係数を再定義する。

(5)    \begin{align*} \frac{A_1}{A} = \cos \alpha \\ \frac{A_2}{A} = \sin \alpha \end{align*}

これに加法定理を考慮して以下の様に単一の正弦関数になる。

(6)    \begin{align*} s &= A ( \sin x \cos \alpha + \cos x \sin \alpha) \\ &= A \sin (x + \alpha) \end{align*}

また、係数を以下の様に再定義する。

(7)    \begin{align*} \frac{A_1}{A} = \sin \beta \\ \frac{A_2}{A} = \cos \beta \end{align*}

この場合は結果が余弦関数となる。同じ変数のsineとcosineは位相が違うだけなので、結果としては当然。

(8)    \begin{align*} s &= A ( \sin x \sin \beta + \cos x \cos \beta) \\ &= A \cos (x - \beta) \end{align*}

 

Runge-Kutta法による1階常微分方程式の解

概要

1階常微分方程式について、4次のRunge-Kutta法(RK4)による数値解と解析解を比較する。Euler法による解についても重ねてみた。

  • RK4は、かなり粗い分割間隔でもよい結果になる
  • Euler法は、分割間隔を小さくとらないと精度が上がらない

Logistic関数

方程式と解析解

Logistic関数に関する微分方程式は以下で表される。

(1)    \begin{equation*} \frac{N(t)}{dt} = r \left( 1 - \dfrac{N(t)}{K} \right) N(t) \end{equation*}

この方程式は解析的に解けて、以下の解を得る。

(2)    \begin{equation*} N(t) = \dfrac{K}{1 + \left( \dfrac{K}{N_0} - 1 \right) e^{-rt}} \end{equation*}

数値解

式(1)をRK4、Euler法によって解き、解析解と重ねる(解き方は1階常微分方程式の数値解法を参照)。

r=10, K=1000として計算した結果は以下の通り。

  • RK4は上限に収束するまでの時間(t=1.0)の20等分くらい(dt=0.05)でもよく合っている
  • Euler法では凹な部分で計算値が過少となり、凸の部分で解析解に追いついているが、適合度はよろしくない。

dtを上記の10分の1(dt=0.05)として計算した結果が以下の通り。Euler法もかなり適合度合いがよくなったが、それでもまだ少し右方にずれている。

逆にdt=0.1と大きくしてみた。折れ線がはっきりわかるくらいの粗さだが、RK4はよくフィットしている。Euler法のずれは大きくなっている。

コード

上記の結果は、以下のコードを実行して得ている。

三角関数

方程式と解析解

簡単な正弦関数を考える。

(3)    \begin{equation*} y = \sin t \end{equation*}

この関数は、以下の微分方程式の解である。

(4)    \begin{equation*} y' = \cos t \quad ( t=0 \Rightarrow y=0 ) \end{equation*}

数値解

式(4)をEuler法、Runge-Kutta法によって解き、式(3)のグラフと重ねる。

以下は分割間隔を1/40周期で計算した結果。RK4の精度はいいが、Euler法では上に凸のときに過大となり、下に凸のときに過少となって初期値に戻っている。

分割幅を1/200周期にするとEuler法でも精度は上がってくるが、まだ少し誤差が残っている。

分割数を10とかなり粗くとっても、RK4の適合度合いはよい。

繰り返し数が増えた時に誤差が累積するかどうか確認してみた。

以下は10万周期分計算した最後の2周忌を取り出した結果で、RK4、Euler法とも傾向は変わらない。規則正しい周期関数の場合は、繰り返し数を増やしても誤差は拡散しないと思われる。

コード

以下はこれまでの計算をしたPythonのコード。

計算区間を長くとって一部だけ表示させるため、Logistic関数のときとは異なる処理をしている。

  • ループの各段階で計算結果をリストに保存せず、各ステップで計算された次のステップの結果を元の変数に入れ替える
  • 予め指定した計算区間の間の結果だけ、表示用のリストに保存する

 

2階常微分方程式の数値解法

概要

2階常微分方程式の数値解法、Euler法とRunge-Kutta法をPythonのモジュールとしてまとめた。

1階常微分方程式の解法についてはこちら

基本形

以下のような方程式を考える。

(1)    \begin{equation*} \dfrac{d^2 y}{dx^2} = f\left ( x, y, \dfrac{dy}{dx}, ... \right ) \end{equation*}

(2)    \begin{equation*} y'' = f(x, y, y', ... ) \end{equation*}

y‘(x)をv(x)と置いて以下のように変形する。

(3)    \begin{align*} \frac{dy}{dx} &=y' = v \\ \frac{dv}{dx} &= v' = f(x, y, y', ...) \end{align*}

Euler法

Euler法では、iステップ目のΔv→Δyと計算し、それらを使ってi+1ステップ目のyvを計算する。

(4)    \begin{align*} \Delta v_i &= \Delata x \cdot f(x_i, y_i, y'_i, ...) \\ \Delta y_i &= \Delta x \cdot v_i \\ y_{i+1} &= y_i + \Delta y_i \\ v_{i+1} &= v_i + \Delta v_i \end{align*}

Runge-Kutta法

Runge-Kutta法も同様に、iステップ目のΔv1→Δy1を計算し、それらを使ってΔvi、Δyi (i = 1, 2, 3, 4)を計算していく。

(5)    \begin{align*} \Delta v_{1, i} &= \Delta x \cdot f(x_i, y_i, v_i, ...) \\ \Delta y_{1, i} &= \Delta x \cdot v_i \\ \Delta v_{2, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{1, i}}{2}, v_i + \dfrac{\Delta v_{1, i}}{2}, ... \right ) \\ \Delta y_{2, i} &= \Delta x \cdot \left (v_i + \dfrac{\Delta v_{1, i}}{2} \right ) \\ \Delta v_{3, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{2, i}}{2}, v_i + \dfrac{\Delta v_{2, i}}{2}, ... \right ) \\ \Delta y_{3, i} &= \Delta x \cdot \left (v_i + \dfrac{\Delta v_{2, i}}{2} \right ) \\ \Delta v_{4, i} &= \Delta x \cdot (v_i + \Delta v_{3, i}) \\ \Delta y_{4, i} &= \Delta x \cdot f(x_i + \Delta x, y_i + \Delta y_{3, i}, v_i + \Delta v_{3, i}, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \frac{\Delta y_{1, i} + 2 \Delta y_{2, i} + 2 \Delta y_{3, i} + \Delta y_{4, i}}{6} \\ v_{i+1} &= v_i + \frac{\Delta v_{1, i} + 2 \Delta v_{2, i} + 2 \Delta v_{3, i} + \Delta v_{4, i}}{6} \end{align*}

数値解法モジュール

コード

Euler法とRunge-Kutta法で2階常微分方程式を解くモジュールを、Pythonで以下の様に書いた。

関数euler2rk4_2の使い方はいずれも同じで、以下の通り。

  • x, y, v:現在の計算ステップでの独立変数x、xにおけるyの値、xにおけるyの1次微分値
  • dx:計算の分割幅
  • func:微分方程式の右辺の関数
  • params:関数funcに与えるパラメーター

パラメーターは、呼び出し側の関数定義で設定したパラメーター名をキーとした辞書で与える。その前提で、モジュール内で**paramsとしてアンパックしている。

戻り値は実数型の数値。

使い方

モジュールの使い方の概要は以下のとおり。

  1. 右辺の関数を定義する
  2. 関数に渡すパラメーターを設定する
  3. 計算開始点と終了点のxの値、分割幅dx, yvの初期値を設定する
  4. x, y, vの計算用リストを確保する
  5. 計算開始点のyvのリストに初期値をセットする
  6. 計算実行

単純な放物線の微分方程式を例にしてモジュールの使い方を示す。まず、解の方程式として以下を考える。

(6)    \begin{equation*} y = a x^2 \quad (-1 \le x \le 1) \end{equation*}

この式を2階微分した、以下の常微分方程式を解いていく。

(7)    \begin{equation*} y'' = 2a \quad \left( x = -1 \Rightarrow \left\{ \begin{array}{l} y = a\\ y' =-2a \end{array} \right) \end{equation*}

まず、必要なパッケージをインポート。diffeqは上で示したモジュールを収めている。

次に、微分方程式右辺の関数を定義する(以下のderivative)。この関数はパラメーターとして係数aを1つ持っている。また、パラメーターの値を辞書で定義する。

数値計算に計算開始点x_startと終了点x_end、分割幅dxyvの計算初期値を設定。

次に、予めxの計算点のリストを用意して、これに対応するy, vの値を計算・保存していく。計算点数は植木算で分割数+1とし、ループ計算に用いるリストは計算終了点を除いている。結果を保存しないなら、リストではなく逐次各ステップの計算値を次のステップの入力値にする。

また、計算開始点におけるy, vの初期値をセットしている。

計算自体はシンプルで、定義域リストの値を1つずつ進めながら、現在の値や関数への参照、パラメーターを与えてモジュールを呼び出す。

実行例

元の解析解と併せて表示した実行結果は以下の通りで、1回微分方程式の場合と同じく、Euler法の解は計算が進むごとに乖離が大きくなっている。Runge-Kutta法の解はかなり正確に解析解と一致している。

以下は、解析解も含めてグラフ表示するための、全体の実行コード。

 

1階常微分方程式の数値解法

概要

1階常微分方程式の数値解法、Euler法とRunge-Kutta法をPythonのモジュールとしてまとめた。

2階常微分方程式の解法についてはこちら

基本形

以下のような方程式を考える。

(1)    \begin{equation*} \dfrac{dy}{dx} = y' = f(x, y, ...) \end{equation*}

Euler法

Euler法はシンプルな方法で、ある計算ステップiの時の関数上の点(xi, yi)とその点における微分係数yiを使って、計算幅dxだけ先にある点xi+1 = xi + dxのに対応するyiの値を近似的に計算する。

(2)    \begin{align*} \Delta y_i &= \Delta x \cdot f(x_i, y_i, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \Delta y_i \end{align*}

概念図からもわかるように、たとえば下に凸の関数の場合、微分係数で線形的に予測した値は常に真値より小さくなるため、誤差が拡大していくことが予想される。

Runge-Kutta法

Runge-Kutta法には1次、2次など異なる次数の解法があるが、4次の方法がよく紹介されている。

(3)    \begin{align*} \Delta y_{1, i} &= \Delta x \cdot f(x_i, y_i, ...) \\ \Delta y_{2, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{1, i}}{2}, ... \right ) \\ \Delta y_{3, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{2, i}}{2}, ... \right ) \\ \Delta y_{4, i} &= \Delta x \cdot f(x_i + \Delta x, y_i + \Delta y_{3, i}, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \frac{\Delta y_{1, i} + 2 \Delta y_{2, i} + 2 \Delta y_{3, i} + \Delta y_{4, i}}{6} \end{align*}

これを説明するのに、xy平面上での関数を描き、その微分係数を複数描いているものが多い。

本来は、右辺のz = f(…)上のある曲線があって、それをxy平面に投影した曲線のある点における微分係数と、その点におけるzの値が等しくなるような曲線を求めるという意味となる。実際にzの値と投影曲線の微分係数が等しいというのを直感的に示すのはやっかいだなと思っていた。

いろいろ探してみると、どうやらこの式形は元の関数をTaylor展開して、次数を上げて精度を高めるというものらしい。文献を見てみると、その導出過程でかなり長々とした式が出ていた。時間がある時に試してみたい。

数値解法モジュール

コード

Euler法とRunge-Kutta法で1階常微分方程式を解くモジュールを、Pythonで以下の様に書いた。

関数euler1rk4_1の使い方はいずれも同じで、以下の通り。

  • x, y:現在の計算ステップでの独立変数xと関数値yの値
  • dx:計算の分割幅
  • func:微分方程式の右辺の関数
  • params:関数funcに与えるパラメーター

パラメーターは、呼び出し側の関数定義で設定したパラメーター名をキーとした辞書で与える。その前提で、モジュール内で**paramsとしてアンパックしている。

戻り値は実数型の数値。

使い方

モジュールの使い方の概要は以下のとおり。

  1. 右辺の関数を定義する
  2. 関数に渡すパラメーターを設定する
  3. 計算開始点と終了点のxの値、分割幅dx, yの初期値を設定する
  4. x, yの計算用リストを確保する
  5. 計算開始点のyのリストに初期値をセットする
  6. 計算実行

単純な放物線の微分方程式を例にしてモジュールの使い方を示す。まず、解の方程式として以下を考える。

(4)    \begin{equation*} y = a x^2 \quad (-1 \le x \le 1) \end{equation*}

この式を1階微分した、以下の常微分方程式を解いていく。

(5)    \begin{equation*} y'= 2ax \quad (x = -1 \Rightarrow y = 1) \end{equation*}

まず、必要なパッケージをインポート。diffeqは上で示したモジュールを収めている。

次に、微分方程式右辺の関数を定義する(以下のderivative)。この関数はパラメーターとして係数aを1つ持っている。また、パラメーターの値を辞書で定義する。

数値計算に計算開始点x_startと終了点x_end、分割幅dxyの計算初期値を設定。

次に、予めxの計算点のリストを用意して、これに対応するyの値を計算・保存していく。計算点数は植木算で分割数+1とし、ループ計算に用いるリストは計算終了点を除いている。結果を保存しないなら、リストではなく逐次各ステップの計算値を次のステップの入力値にする。

また、計算開始点におけるyの初期値をセットしている。

計算自体はシンプルで、定義域リストの値を1つずつ進めながら、現在の値や関数への参照、パラメーターを与えてモジュールを呼び出す。

実行例

元の解析解と併せて表示した実行結果は以下の通りで、予想通りEuler法の解は計算が進むごとに乖離が大きくなっている。Runge-Kutta法の解はかなり正確に解析解と一致している。

以下は、解析解も含めてグラフ表示するための、全体の実行コード。