Introduction

** The ordinary differential equation for $ y (x) $ is linearly linear, and the parameters ([eigenvalue](https://ja.wikipedia.org/wiki/%E5%9B%BA%E6] are as follows. % 9C% 89% E5% 80% A4)) If $ \ lambda $ is also linear, then the Boundary Value Problem (https://ja.wikipedia.org/wiki/%E5%A2%83%E7) % 95% 8C% E5% 80% A4% E5% 95% 8F% E9% A1% 8C) is reduced to the eigenvalue question of the matrix, and the solution of the matrix equation solves the solution (eigenfunction) and eigenvalue of the ordinary differential equation. It is possible to determine [1]. ** **

** It is easier to implement than the method of integrating and solving ordinary differential equations as an initial value problem (shooting method [2], etc.). ** **

Here, let us consider a second-order ordinary differential equation.

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}

As a simultaneous boundary condition

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

think of.

Divide the interval $ x = [x_0, x_1] $ into N equal parts

\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}

And.

Now, when the derivatives $ y'' (x) $ and $ y'(x) $ are evaluated by the center difference approximation **

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}

Can be expressed as. The eigenvalues and eigenfunctions when solving (1) by these approximations include an error of $ O (\ delta x ^ 2) $.

By substituting Eqs. (5, 6) into Eq. (1) and rearranging them, the problem can be expressed as a ** general eigenvalue problem ** in matrix format.

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

Where $ \ mathbf {y} $ is the solution vector. \mathbf{y}=(y_1, y_2, y_3, ...., y_{N-1}) {\tag 9}

$ A = (a_ {ij}) $ and $ B = (b_ {ij}) $ are $ N-1 $ dimensional matrices, which are triple diagonal matrices. The matrix elements are $ 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}}

It is expressed as.

** Equation (8) can be solved using Python libraries (numpy, scipy, etc.) [3]. ** **

In this paper, I will try to solve a simple example using this method.

This method is different from the function expansion method ( spectral method ) in which the solution is expanded by an appropriate functional system and the matrix equation for the coefficient is obtained.

** Boundary value problem of second-order ordinary differential equations (simple vibration) **

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

Solve ** by the matrix method. ** **

In this case, each function in equation (1) is

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

And

In $ A \ \ mathbf {y} = \ lambda \ B \ \ mathbf {y} $, $ A $ is

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}} And

$ B $ is ** identity matrix * * By $ E $

B = -E{\tag {16}} It is expressed as.

** And for these matrices ** $ A \ \mathbf{y_n} = - \lambda_n \ E \ \mathbf{y} {}\tag {17} $

Solve ** to find the eigenvalue $ \ lambda_n $ and the corresponding eigenfunction $ y_n (x) $. ** **

code

Solve the general eigenvalue problem (Equation 17) using ** scipy.linalg.eig ** [3].

"""
Solving the boundary value problem by the matrix method
"""
import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt

delta_x = 0.025
x0, x1 = 0, 1 #Boundary x coordinate
N=int((x1-x0)/delta_x) #Number of data grids
print("N=",N)

y = np.zeros([N-1,N+1])

#Simultaneous boundary conditions
y[:,0] = 0
y[:,-1] = 0

A=np.zeros([N-1,N-1])
B=-np.identity(N-1) # B= -E

#Construction of matrix A
for i in range(N-1):  #Be careful of the index position because it is a triple diagonal matrix.
    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) #Solve the general eigenvalue problem and find the eigenvalues"eigen"， Eigenfunction"vec"Store in the object

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') #x-axis label
plt.ylabel('Y') #y-axis label

plt.show()

result

The first three are displayed in ascending order of eigenvalues. From the left, it is $ \ lambda_1 $, $ \ lambda_2 $, $ \ lambda_3 $. eigen values= [ 9.86453205+0.j 39.39731010+0.j 88.41625473+0.j]

The eigenfunctions $ y_1 (x) $, $ y_2 (x) $, $ y_3 (x) $ corresponding to these eigenvalues are shown in the figure.

References

[1] Ryo Natori, ["Numerical Analysis and Its Applied Computer Mathematics Series 15"](https://www.amazon.co.jp/%E6%95%B0%E5%80%A4%E8%A7%A3 % E6% 9E% 90% E3% 81% A8% E3% 81% 9D% E3% 81% AE% E5% BF% 9C% E7% 94% A8-% E3% 82% B3% E3% 83% B3% E3% 83% 94% E3% 83% A5% E3% 83% BC% E3% 82% BF% E6% 95% B0% E5% AD% A6% E3% 82% B7% E3% 83% AA% E3% 83% BC% E3% 82% BA-% E5% 90% 8D% E5% 8F% 96-% E4% BA% AE / dp / 4339025488), Corona Publishing, 1990.

[2] Example of solving ordinary differential equations using shooting method: [[Science / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (2), harmonic oscillator potential, quantum mechanics](http //qiita.com/sci_Haru/items/edb475a6d2eb7e901905)

[3] How to find the eigenvalues and eigenvectors (functions) of a matrix using a library: [[Scientific and technical calculation by Python] Solving (generalized) eigenvalue problems using numpy and scipy, using the library](http: / /qiita.com/sci_Haru/items/034c6f74d415c1c10d0b)

[Scientific / technical calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation

Introduction

This method is different from the function expansion method ( spectral method ) in which the solution is expanded by an appropriate functional system and the matrix equation for the coefficient is obtained.

Contents

code

result

References

[Scientific / technical calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation

Introduction

This method is different from the function expansion method (** spectral method **) in which the solution is expanded by an appropriate functional system and the matrix equation for the coefficient is obtained.

Contents

code

result

References

This method is different from the function expansion method ( spectral method ) in which the solution is expanded by an appropriate functional system and the matrix equation for the coefficient is obtained.