소개

$y(x)$에 대한 상미분 방정식이 동차 선형이고, 다음과 같이 파라미터( 고유값 )$\lambda$에 대해서도 선형인 경우는 경계값 문제 을 행렬의 고유치문에 귀착 그런 다음 행렬 방정식을 풀면 상미분 방정식의 해 (고유 함수)와 고유 값을 결정할 수 있습니다 [1].

초기치 문제로서 상미분 방정식을 적분해 풀는 방법(사격법[2] 등)보다 실장이 간단하다.

여기서는 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)) {\tag 1}$

동시 경계 조건으로서,

$y(x_0)=0, y(x_N)=0 {\tag2}$

생각합니다.

구간 $x= [x_0,x_1]$를 N 등분

$\delta x = (x_1-x_0)/N{\tag 3}$

$x_i=x_0+i\\delta x{\tag 4}$

$y_i=y(x_i) {\tag 5}$

한다.

이제 미분 $y ''(x) $ 및 $ y'(x) $를 중심 차분 근사로 평가하면

$y'(x) =\frac{y_{i+1}-y_{i-1}}{2\delta x} +O(\delta x^3){\tag 6}$

$y''(x) =\frac{y_{i+1}-2y_{i}+y_{i-1}}{\delta x ^2}+O(\delta x^3) {\tag 7 }$

라고 표현할 수 있다.
이러한 근사에 의해 (1)을 풀었을 때의 고유치·고유 함수에는 $O(\delta x^2)$의 오차가 포함된다.

식 (5, 6)을 식 (1)에 대입하고 정리하면 문제를 행렬 형식에 따른 일반 고유치 문제로 표현할 수 있다.

$A\\mathbf{y}=\lambda\B\\mathbf{y}\{\tag 8}$

여기서 $\mathbf{y}$는 해 벡터이다.
$\mathbf{y}=(y_1, y_2, y_3, ...., y_{N-1}) {\tag 9}$

$A=(a_{ij})$ 및 $B=(b_{ij})$는 $N-1$차원의 행렬이며, 3중 대각 행렬이다.
이러한 행렬 요소는 $i=1,2,...,N-1$로,

$a_{i,i-1}=\frac{p_i}{\delta x^2}+\frac{q_i}{2\delta x},\a_{i,i}=\frac{-2p_i}{\delta x^2}+r_i,\a_{i,i+1}=\frac{p_i}{\delta x^2}-\frac{q_i}{2\delta x} {\tag {10}} $
$a_{i,j} = 0\( j\neq i-1, i, i+1){\tag {11}}$

$b_{i,i-1}=\frac{u_i}{\delta x^2}+\frac{v_i}{2\delta x},\b_{i,i}=\frac{-2u_i}{\delta x^2}+w_i,\b_{i,i+1}=\frac{u_i}{\delta x^2}-\frac{v_i}{2\delta x} {\tag {12}} $
$b_{i,j} = 0\( j\neq i-1, i, i+1){\tag {13}}$

로 표시된다.

식 (8)은 파이썬 라이브러리 (numpy, scipy 등)를 이용하여 풀 수있다 [3].

이 논문에서는이 방법을 사용하여 간단한 예제를 해결합니다.

#이 방법은 적절한 함수 시스템에서 솔루션을 확장하고 계수에 대한 행렬 방정식을 찾는 함수 확장 방법 (스펙트럼 방법)과 다릅니다.

내용

2층 동차 상미분 방정식의 경계값 문제(단진동)

$y''(x)=-\lambda y,\(y(0)=y(1)=0)$

행렬 방식으로 풀기.

이 경우 식 (1)의 각 함수는

$p(x)=1, q(x)=0, r(x)=0, u(x)=0, v(x)=0, w(x)=1 {\tag {14}}$

되고,

$ A\\mathbf {y} =\lambda\B\\mathbf {y} $에서 $ A $는

$a_{i,i-1}=\frac{1}{\delta x^2},\a_{i,i}=\frac{-2}{\delta x^2},\a_{i, i+1}=\frac{1}{\delta x^2} {\tag {15}}$
그리고,

$B$는 단위 행렬 $E$

$B = -E{\tag {16}}$
로 표시된다.

나머지는이 행렬에
$A\\mathbf{y_n} = -\lambda_n\E\\mathbf{y}{}\tag{17}$

를 풀고 고유치 $\lambda_n$와 그에 대응하는 고유함수 $y_n(x)$를 구한다.

코드

scipy.linalg.eig를 사용하여 일반 고유치 문제 (식 17)를 푸십시오 [3].

"""
行列法による境界値問題の解法
"""
import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt

delta_x = 0.025
x0, x1 = 0, 1 #境界のx座標
N=int((x1-x0)/delta_x) # データグリッド数
print("N=",N)

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) # B= -E

#行列Aの構築
for i in range(N-1):  # 3重対角行列なのでindexの位置を気をつけること
    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) # 一般固有値問題を解き，固有値を"eigen"，固有関数を"vec"というオブジェクトへ格納する

print("eigen values=",eigen)
#print("eigen vectors=",vec)

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

#
# for plot
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') # y軸のラベル

plt.show()

결과

고유치를 작은 순서로 처음 3개를 표시한다. 왼쪽에서 $\lambda_1$, $\lambda_2$, $\lambda_3$이다.
eigen values= [ 9.86453205+0.j 39.39731010+0.j 88.41625473+0.j]

이들 고유치에 대응하는 고유함수 $y_1(x)$, $y_2(x)$, $y_3(x)$를 도면에 나타낸다.

참고문헌

[1] 나토리 료, 「수치 해석과 그 응용 컴퓨터 수학 시리즈 15」, 코로나 사, 1990.

[2] 사격법을 이용한 상미분 방정식의 해법 예 : [Python에 의한 과학·기술 계산] 정상상태의 1차원 슈레딩거 방정식의 사격법에 의한 해법(2), 조화진동자 포텐셜, 양자역학

[3] 라이브러리를 사용하여 행렬의 고유 값 및 고유 벡터 (함수)를 찾는 방법 : [Python에 의한 과학·기술 계산] numpy·scipy를 이용한(일반화) 고유치 문제의 해법, 라이브러리 이용