はじめに

微分方程式は、物体の運動を叙述する場合やベクトル場を描写したい場合などに用いられる。そこで、今回は$xy$平面において$y'=f(x,y)$で与えられる微分方程式を数値計算してベクトル場を描写する。ただし、$y'=f(x,y)$は、変数分離型で解くことのできるものを扱い、変数分離型の解法を用いて解いた解と計算機での数値計算結果を比較して齟齬が無いことを確認する。

例題1(同心円状のベクトル場)

$\frac{dy}{dx}=-\frac{x}{y}$を変数分離で解き、数値計算での結果と比較する。

変数分離型の解法

$ydy=-xdx$と変形することができるので、両辺を以下のように積分する。

\int ydy=-\int xdx

\frac{1}{2}y^2+C_1=-\frac{1}{2}x^2+C_2

ただし、$C_1,C_2$は積分定数とする。上式を整理すると以下のことが分かる。

x^2+y^2=r^2

ただし、$r\ge 0$とする。

プログラム

以下のようなプログラムを作成する。

python differential_equation.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
# ドット数
n=1000
#xy座標系における描写範囲
L=1
# x軸の刻み幅
dx= 2*L/(n-1)

x=np.linspace(-L, L, n)
y=np.linspace(-L, L, n)
X, Y = np.meshgrid(x, y)
U = np.zeros((n, n))
V = np.zeros((n, n))
for p in range(n):
    for q in range(n):
        xx=X[p][q]
        yy=Y[p][q]
        #微分方程式
        f_y=-xx/yy
        U[p][q]=dx
        V[p][q]=f_y*dx        

plt.streamplot(X, Y, U, V, color='blue', density=1.5, linewidth=0.5, arrowsize=1.5)
#plt.quiver(X, Y, U, V, color='blue', angles='xy', scale_units='xy', scale=1)
plt.savefig('流れ.png')
plt.show()

これを実行すると以下のようなグラフが出力される。

同心円状の流れの場を描写することができた。
したがって、以上の結果から理論解と数値計算による結果は一致する。

例題2(原点を通る直線状のベクトル場)

$\frac{dy}{dx}=\frac{y}{x}$を変数分離で解き、数値計算での結果と比較する。

変数分離型の解法

 \frac{1}{y}dy=\frac{1}{x} dx

と変形できるので、両辺を以下のように積分する。

 \int \frac{1}{y}dy= \int \frac{1}{x} dx

 \ln |y|+C_1 = \ln |x| +C_2

 \ln |\frac{y}{x}| = C

  |\frac{y}{x}| = \exp(C)

y=ax

ただし、$C_1,C_2,C,a$は実数の定数とする。

プログラム

以下のようなプログラムを作成する。

python differential_equation2.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
# ドット数
n=1000
#xy座標系における描写範囲
L=1
# x軸の刻み幅
dx= 2*L/(n-1)

x=np.linspace(-L, L, n)
y=np.linspace(-L, L, n)
X, Y = np.meshgrid(x, y)
U = np.zeros((n, n))
V = np.zeros((n, n))
for p in range(n):
    for q in range(n):
        xx=X[p][q]
        yy=Y[p][q]
        #微分方程式
        f_y=yy/xx
        U[p][q]=dx
        V[p][q]=f_y*dx        

plt.streamplot(X, Y, U, V, color='blue', density=1.5, linewidth=0.5, arrowsize=1.5)
#plt.quiver(X, Y, U, V, color='blue', angles='xy', scale_units='xy', scale=1)
plt.savefig('流れ2.png')
plt.show()

原点を通る直線状の流れの場を描写することができた。
したがって、以上の結果から理論解と数値計算による結果は一致する。

まとめ

今回は、$xy$座標系において微分方程式を解くことで、流れの場を描写するということに挑戦した。具体的には、まず、変数分離型で解くことのできる微分方程式を2題扱い、それを理論的に解いた。次に、その理論解と比較するために、計算機を用いた数値計算を実行して流れの場を描写した。最後に理論解と数値計算解に差がないことを確認した。微分方程式は、理論的に解くことができない場合が圧倒的に多い。したがって、数値計算が重要になる場合も多い。

参考文献