更新:2024/09/26

SymPyを使用して常微分方程式を解く方法について

はるか
はるか
SymPyで微分方程式の解き方、わかってる?
ふゅか
ふゅか
うん!微分方程式はsp.Eq()を使って定義して、それをdsolveで解けばいいんだよね!

1. 微分方程式とSymPy

PythonのライブラリであるSymPyは、数学的な計算や解析に使える便利なツールです。常微分方程式を解く機能も持っており、手軽に微分方程式の解を求めることができます。

1.1. 基本的な考え方

微分方程式を扱う際には、次の3つのステップが重要です。

  1. 微分の表現・・・関数 \( y(x) \) を変数 \( x \) について微分する場合、y(x).diff(x) と記述します。これで関数 \( y(x) \) の1階微分を表現します。
  2. 微分方程式の表現・・・微分方程式は、sp.Eq()を使って定義します。このメソッドに左辺と右辺の式を渡すことで、等式としての微分方程式を作成できます。
  3. 微分方程式を解く・・・作成した微分方程式は、sp.dsolve() を使って解きます。

これらのステップを理解すれば、SymPyを使って簡単に微分方程式を解くことができます。

2. SymPyで微分方程式を解く方法

2.1. 常微分方程式の解法

まずは、1階常微分方程式の具体例を通して、SymPyの使い方を見ていきましょう。

例: 1階常微分方程式

$$\frac{dy}{dx} + y = x$$

この1階常微分方程式をSymPyで解くには、以下のように記述します。

import sympy as sp
# 変数と関数の定義
x = sp.symbols('x')
y = sp.Function('y')

# 微分方程式の定義
ode = sp.Eq(y(x).diff(x) + y(x), x)

# 微分方程式の解法
solution = sp.dsolve(ode)

print("解:", solution)

このコードは、微分方程式 \(\frac{dy}{dx} + y = x\) を解き、結果を表示します。

2.2. 初期条件付き常微分方程式

次に、初期条件付きの微分方程式を解く例を見てみましょう。初期条件も指定することで、特定の解を求めることができます。

例2: 初期条件付き常微分方程式

\(\frac{dy}{dx} = -2y\) かつ \( y(0) = 1 \)

以下は、上記の常微分方程式を解く例です。

# 微分方程式の定義
ode = sp.Eq(y(x).diff(x), -2*y(x))

# 初期条件の定義
initial_conditions = {y(0): 1}

# 微分方程式の解法(初期条件付き)
solution = sp.dsolve(ode, ics=initial_conditions)

print("解:", solution)

このコードは、微分方程式 \(\frac{dy}{dx} = -2y\) を初期条件 \( y(0) = 1 \) で解き、結果を表示します。

2.3. 2階常微分方程式

次に、2階常微分方程式をSymPyで解く例を示します。2階微分方程式では、関数の2階微分を使って方程式を定義します。

$$\frac{d^2y}{dx^2} - 3\frac{dy}{dx} + 2y = 0$$

以下は、上記の常微分方程式を解く例です。

# 微分方程式の定義
ode = sp.Eq(y(x).diff(x, x) - 3*y(x).diff(x) + 2*y(x), 0)

# 微分方程式の解法
solution = sp.dsolve(ode)

print("解:", solution)

このコードは、微分方程式 \(\frac{d^2y}{dx^2} - 3\frac{dy}{dx} + 2y = 0\) を解き、結果を表示します。出力は次のようになります。

はるか
はるか
2階微分の場合は、y(x).diff(x, x)と書く。

2.4. 連立微分方程式の解法

最後に、連立微分方程式を解く方法を紹介します。複数の微分方程式を同時に解く場合も、SymPyで簡単に扱うことができます。

$$\frac{dy_1}{dx} = y_1 + y_2$$

$$\frac{dy_2}{dx} = y_1 - y_2$$

以下は、上記の連立微分方程式を解く例です。

# 関数の定義
y1, y2 = sp.Function('y1'), sp.Function('y2')

# 連立微分方程式の定義
ode1 = sp.Eq(y1(x).diff(x), y1(x) + y2(x))
ode2 = sp.Eq(y2(x).diff(x), y1(x) - y2(x))

# 連立微分方程式の解法
solution = sp.dsolve((ode1, ode2))

print("解:", solution)

このコードは、上記の連立微分方程式を解き、結果を表示します。

PR