Pythonを用いた微分方程式と行列の解法

Pythonは科学計算を得意として発展してきた歴史があります。そのため、微分方程式の数値解析は得意分野のひとつでもあるのです。

行列形式による常微分方程式の解法

行列形式による常微分方程式の解法は、初期値問題として常微分方程式を積分して解く方法よりも実装が簡単である。具体的には、以下の2階常微分方程式を考えます。

$$p(x)y”(x)+q(x)y'(x)+r(x)y(x)=\lambda(u(x)y”(x)+v(x)y'(x)+w(x)y(x))$$

同時境界条件として、$y(x_0)=0, y(x_N)=0$ を考えます。

この問題を行列形式による一般固有値問題として表現できます。

import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt

delta_x = 0.025
x0, x1 = 0, 1
N=int((x1-x0)/delta_x)
y = np.zeros([N-1,N+1])
y[:,0] = 0
y[:,-1] = 0
A=np.zeros([N-1,N-1])
B=-np.identity(N-1)

for i in range(N-1):
    if i == 0:
        A[i,i] = -2/(delta_x**2)
        A[i,i+1] = 1/(delta_x**2)
    elif i == N-2:
        A[i,i-1] = 1/(delta_x**2)
        A[i,i] = -2/(delta_x**2)
    else:
        A[i,i-1] = 1/(delta_x**2)
        A[i,i] = -2/(delta_x**2)
        A[i,i+1] = 1/(delta_x**2)

eigen, vec= scipy.linalg.eig(A,B)

for j in range(N-1):
    for i in range(1,N):
        y[j, i] = vec[i-1,j]

X= np.linspace(x0,x1, N+1)
plt.plot(X, y[0,:],'-',markersize=5,label='y1')
plt.plot(X, y[1,:],'-',markersize=5,label='y2')
plt.plot(X, y[2,:],'-',markersize=5,label='y3')
plt.legend(loc='upper right')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

このコードは、行列形式による常微分方程式の解法を示しています。

まとめ

Pythonを用いて微分方程式と行列の解法を学ぶことは、科学計算やデータ分析において非常に有用です。Pythonのライブラリを活用することで、複雑な問題も簡単に解くことが可能になります。.