はじめに：Python で数値解析

葛西 真寿　弘前大学大学院理工学研究科

Python で数値解析のために必要なプログラミングと，いくつかの応用例を示します。 ここでは主に SymPy を使った例を示し，SciPy の利用についても解説します。

Python の文末は，単に改行するか，; で終わって改行するかです。; で終わる場合は結果が表示されません。

1 + 2

3

1 + 2;

条件分岐

条件分岐の if 文の書式は以下のようになっています。

if 条件式 1:
    条件式 1 が満たされた場合に実行する文
elif 条件式 2: 
    条件式 2 が満たされた場合に実行する文
else:
    それ以外の場合に実行する文

elseif かな？と思わせて elif と書きます。

import random
x = random.random()
if x <= 0.5:
    print(x, " は 0.5 以下")
else:
    print(x, " は 0.5 より大きい")

0.4432030865092451  は 0.5 以下

繰り返し処理

$\displaystyle \sum_{i=1}^{5} i\ $ を計算する例

for 文

range() 関数を使った for 文の書式は以下の通りです。開始値を省略すると 0 とみなし，増分を省略すると 1 とみなします。

for i in range(開始値, 終端値, 増分):
    i が開始値以上，終端値未満の場合に実行される文

sum = 0
for i in range(1, 6, 1):
    sum = sum + i
print(sum)

15

while 文

while 文の書式は以下の通りです。

while 条件式:
    条件式が成り立つ場合に実行する式

sum = 0
i = 1
while i <=5:
    sum = sum + i
    i = i + 1
print(sum)

15

再帰的定義関数

$\displaystyle \mbox{sum}(n) = \sum_{i=1}^{n} i , \ (n \geq 1)\ $ を以下のように考えます。

もし，$n = 1$ なら $\mbox{sum}(n) = \mbox{sum}(1) = 1$

もし，$n > 1$ なら $\mbox{sum}(n) = \mbox{sum}(n-1) + n$

def sum(n):
    if n == 1:
        return 1
    else:
        return sum(n-1) + n

print(sum(5))

15

参考：Sum()

上の例で，sum(n) という関数を再帰的に定義しましたが，定義などしなくても，実は SymPy を import すると，数列の和を計算する関数 Sum() が使えます。

from sympy import *
from sympy.abc import *

Sum(k, [k, 1, 5])

$\displaystyle \sum_{k=1}^{5} k$

あれ，計算してくれませんね。以下のように .doit() をつけてみます。

Sum(k, [k, 1, 5]).doit()

$\displaystyle 15$

一般の n までの和も...

Sum(k, [k, 1, n]).doit()

$\displaystyle \frac{n^{2}}{2} + \frac{n}{2}$

Sum(k**2, [k, 1, n]).doit()

$\displaystyle \frac{n^{3}}{3} + \frac{n^{2}}{2} + \frac{n}{6}$

数値微分

解析的な微分の定義は（前方差分の極限として） $$\frac{df}{dx} \equiv \lim_{h \rightarrow 0} \frac{f(x+h) - f(x)}{h}$$ でした。滑らかな関数であれば，（後方差分の極限として） $$\frac{df}{dx} = \lim_{h \rightarrow 0} \frac{f(x) - f(x-h)}{h}$$ としても良いです。

現実世界では，$h\rightarrow 0$ の極限はとれませんから十分小さい値として h を定義し，近似的な数値微分を（前方差分と後方差分の平均である中心差分として）以下のように定義しましょう。

$$\begin{eqnarray} \frac{df}{dx} &\simeq& \frac{1}{2} \left\{\frac{f(x+h) - f(x)}{h} + \frac{f(x) - f(x-h)}{h} \right\} \\ &=& \frac{f(x+h) - f(x-h)}{2h} \end{eqnarray}$$

SymPy の import

# SymPy を使うときのおまじない
from sympy import *
from sympy.abc import *
from sympy import pi

# 以下はグラフを SVG で Notebook にインライン表示させる設定
from IPython.display import set_matplotlib_formats
%matplotlib inline
set_matplotlib_formats('svg')

以下の例では，$f(x) = \sin x$ に対して $f'(x) = \cos x$ と，上記の数値微分との誤差を描いています。

def f(x):
    return sin(x)

def diff_f(x, h):
    return (f(x+h)-f(x-h))/(2.0*h)

plot(cos(x) - diff_f(x, 1.0E-6), (x, -2*pi, 2*pi));

以下の例からわかるように，h をやみくもに小さくすればするほど，数値微分の精度があがる，というわけでもなさそうです。（縦軸の目盛りに注目。）

plot(cos(x) - diff_f(x, 1.0E-7), (x, -2*pi, 2*pi));

参考：diff()

実は数値微分などしなくても，Python は sympy を import すると，関数の微分をしてくれます。

def f(x):
    return sin(x)

diff(f(x), x)

$\displaystyle \cos{\left(x \right)}$

SciPy による数値微分

scipy.misc.derivative()

SciPy には数値微分してくれる関数 scipy.misc.derivative() が用意されています。それを使ってみましょう。（すでに定義されている函数を使うほうが楽ちん。）

import scipy.misc
# scipy.misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3)

ans = scipy.misc.derivative(f, 3, 1.0E-6)

print(cos(3.0))
print(diff_f(3.0, 1.0E-6))
print(ans)

-0.989992496600445
-0.989992496730485
-0.989992496730485

数値積分

シンプソン法で $\displaystyle \int_a^b f(x)\, dx$ を求める。

積分区間 $[a, b]$ を $N (=2n)$ 等分（偶数等分）し， $$h = \frac{b-a}{N}, \quad f_i = f(a + i\,h)$$ とする。シンプソン法の公式は $$\begin{eqnarray} \int_a^b f(x)\, dx &\simeq& \frac{h}{3}( f_0 + 4 f_1 + f_2) + \frac{h}{3}( f_2 + 4 f_3 + f_4) + \cdots\\ &&\quad \cdots + \frac{h}{3}( f_{N-2} + 4 f_{N-1} + f_N)\\ &=&\frac{h}{3} \left(f_0 + f_{2n} + 2 \sum_{i=1}^{n-1} f_{2i} + 4 \sum_{i=1}^{n} f_{2n-1}\right) \end{eqnarray}$$

for 文による解法

$\displaystyle \int_0^{\pi/2} \sin x \,dx$ の定積分を，積分区間を $N=20$ 分割してシンプソン法で解く例です。

for 文の繰り返し処理で和をとっています。

def f(x):
    return sin(x)

a = 0.0         # 積分の下端。
b = float(pi/2) # 積分の上端。
N = 20          # 分割数 N は偶数とすること（シンプソン法の決め事）

# ==== 
h = (b-a)/N    

def fi(i):
    return f(a + i*h)

simpson_f = 0

for i in range(2, N+1, 2):
    simpson_f  = simpson_f  + h/3*(fi(i-2) + 4*fi(i-1) + fi(i))

print(N, " 分割の数値解は ", simpson_f)

# ==== 精度確認のため，分割数を倍にして計算
N = 2*N
h = (b-a)/N    

def fi(i):
    return f(a + i*h)

simpson_f = 0

for i in range(2, N+1, 2):
    simpson_f  = simpson_f  + h/3*(fi(i-2) + 4*fi(i-1) + fi(i))

print(N, " 分割の数値解は ", simpson_f)

20  分割の数値解は  1.00000021154659
40  分割の数値解は  1.00000001321438

再帰定義関数による解法

$\displaystyle \int_0^{\pi/2} \sin x \,dx$ の定積分を，積分区間を $N=20$ 分割してシンプソン法で解く例です。

再帰的定義関数で和をとっています。

1行が長すぎるような場合，継続行を示す \ で改行して次の行に続けます。

a, b, N = symbols('a b N')

def f(x):
    return sin(x)

def simp_f(a, b, N, i):
# a  積分の下端
# b  積分の上端
# N  分割数 N は偶数とすること（シンプソン法の決め事）
# 実際の計算の際は，simp_f(a, b, N, N) と同じ N を2回書く
    h = float((b-a)/N)
    if i == 2:
        return h/3*(f(a+(i-2)*h)+4*f(a+(i-1)*h)+f(a+i*h))
    else:
        return h/3*(f(a+(i-2)*h)+4*f(a+(i-1)*h)+f(a+i*h)) \
               + simp_f(a,b,N,i-2)

N = 20; print(N, " 分割の数値解は ", simp_f(0, pi/2, N, N))
N = 40; print(N, " 分割の数値解は ", simp_f(0, pi/2, N, N))

20  分割の数値解は  1.00000021154659
40  分割の数値解は  1.00000001321438

参考：integrate(), evalf()

実は数値微分などしなくても，Python は sympy を import すると，関数の積分をしてくれます。

integrate(f(x), (x, 0, pi/2))

$\displaystyle 1$

解析的に積分できない場合は，数値積分によって近似値を求めることができます。

$\displaystyle \int_0^{\pi} \sin(\sin x)\, dx$ の計算をします。解析的に解けないと，そのままで返します。

integrate(sin(sin(x)), (x, 0, pi))

$\displaystyle \int\limits_{0}^{\pi} \sin{\left(\sin{\left(x \right)} \right)}\, dx$

直前の結果を表す _ に .evalf() すると数値解を返します。

_.evalf()

$\displaystyle 1.78648748195005$

SciPy による数値積分

scipy.integrate.quad()

SciPy には数値積分してくれる関数 scipy.integrate.quad() が用意されています。それを使ってみましょう。（すでに定義されている函数を使うほうが楽ちん。）結果にはデフォルトで積分値と推定誤差が出力されます。

import scipy.integrate
import numpy as np

def sinsin(x):
    return np.sin(np.sin(x))

x0 = 0
x1 = np.pi # numpy の円周率 π
scipy.integrate.quad(sinsin, x0, x1)

(1.786487481950052, 2.4687274069628184e-11)

方程式の数値解

$f(x) = 0$ の解を求める。

まず，何らかの方法で（例えば，plot f(x) としてグラフを描いてみて $x$ 軸と交差する点のあたりをつけるなどして）$f(x) = 0$ の解である $x$ に近いと思われる値 $x_k$ がわかったとする。

この $x_k$ のまわりに $f(x)$ をテイラー展開すると， $$f(x) \simeq f(x_k) + f'(x_k)\,(x - x_k)$$

$f(x) = 0$ であるから，上式の左辺をゼロとおいて，$x$ について解くと， $$x \equiv x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$$

これを $|f(x_{k+1})| < \epsilon$ （ここで $\epsilon$ 十分小さい数値，たとえば $1.0\times 10^{-8}$ とか？）になるまで反復計算を行い，このときの $x_{k+1}$ をもって，$f(x) = 0$ の近似解とする。

$f'(x_k)$ を求めるのが面倒な場合は，以下のようにして数値的な微分で代用する。 $$f'(x_k) \simeq \frac{f(x_k + h) - f(x_k - h)}{2 h}, \quad |h| \ll 1$$

つまり， $$x_{k+1} = x_k - \frac{f(x_k)\times 2 h}{f(x_k + h) - f(x_k - h)}$$

まず，解を求める関数 $f(x)$ を定義し，$f(x) = 0$ となる $x$ のあたりをつけるためにグラフを描いてみます。

以下では例としてベッセル関数 besselj(0,x) のゼロ点を求めています。ベッセル関数 besselj() は from sympy import * すれば使えます。

def f(x):
    return besselj(0,x)

plot(f(x), (x, 0, 10));

while 文による解法

$f(x) = 0$ となる $x$ の1つは，$2$ から $3$ の間にあることがわかります。

探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件 h は，数値解の精度に関係しますので，いくつか値を変えて出力させます。

# x には探索の初期値を入れる
# h は解の収束条件である微小量と数値微分をするときに使う微小量を兼ねる。
def f(x):
    return float(besselj(0,x))

x = 2.5
h = 1E-6
while abs(f(x)) > h:
    x = x - float(f(x)*2*h/(f(x + h)-f(x - h)))

print(x)

2.404824591791648

print(x)
# ==== 精度確認のため，h を変えて計算
x = 2.5
h = 1E-12
while abs(f(x)) > h:
    x = x - f(x)*2*h/(f(x + h)-f(x - h))

print(x)

2.404824591791648
2.4048255576957636

再帰定義関数による解法

以下のように，再帰定義により $f(x) = 0$ の解を求める関数 newton_f(x, h) を定義します。

探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件 h は，数値解の精度に関係しますので，いくつか値を変えて出力させます。

# 再帰定義による解
# x には探索の初期値を入れる
# h は解の収束条件である微小量と数値微分をするときに使う微小量を兼ねる。
def f(x):
    return float(besselj(0,x))

def newton_f(x, h):
    if abs(f(x)) > h:
        return newton_f(x-f(x)*2*h/(f(x+h)-f(x-h)), h)
    else:
        return x

print(newton_f(2.5, 1.E-6))
print(newton_f(2.5, 1.E-8))
print(newton_f(2.5, 1.E-10))
print(newton_f(2.5, 1.E-12))
print(newton_f(2.5, 1.E-14))

2.404824591791648
2.4048255576955846
2.404825557695499
2.4048255576957636
2.404825557695766

解の収束条件 h を変えても変わらない桁までが，精度のよい数値解ということになります。上の例では，

2.40482555769576

くらいでしょうか。

besselj(0, 2.40482555769576).evalf()

$\displaystyle 6.62479860154324 \cdot 10^{-15}$

参考：nsolve

実は Python は sympy を import すると，関数の積分をしてくれる関数 nsolve() が使えます。

x = Symbol('x')

nsolve(besselj(0,x), x, 2.5, prec=15)

$\displaystyle 2.40482555769577$

besselj(0,_).evalf()

$\displaystyle -6.10876525973673 \cdot 10^{-17}$

微分方程式の数値解法

4次のルンゲ・クッタ法で常微分方程式を解きます。

1階常微分方程式

1階の常微分方程式 $$\frac{dx}{dt} = f(x, t)$$ を4次のルンゲ・クッタ法で解く例です。

初期条件を $t=t_0$ で $x(t_0) = x_0$ とし，$t = t_1$ までを $N$ 分割して解きます。

$$h = \frac{t_1 - t_0}{N}, \quad t_i = t_0 + i\, h, \quad x_i = x(t_i)$$

とし，以下の式から次のステップの値 $x_{i+1}$ を決めます。

$$\begin{eqnarray} k_1 &=& h\, f(x_i, t_i) \\ k_2 &=& h\, f(x_i + k_1/2, t_i + h/2) \\ k_3 &=& h\, f(x_i + k_2/2, t_i + h/2) \\ k_4 &=& h\, f(x_i + k_3, t_i + h) \\ x_{i+1} &=& x_i + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) \end{eqnarray}$$

例：ロジスティック方程式

例として，簡単なロジスティック方程式

$$\frac{dx}{dt} = (1-x)\,x$$

を，初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 0.1$，$t = t_0 = 0$ から $t_1 = 10$ までを $N = 100$ 分割して解きます。

def f(x, t):
    return (1-x)*x

t0 = 0
x0 = 0.1
t1 = 10
N = 100

h = (t1 - t0)/N

t = t0
x = x0
# リスト T, X を初期化し，初期条件を append
T = []
X = []
T.append(t)
X.append(x)

for i in range(1, N+1):
    k1 = h*f(x, t)
    k2 = h*f(x+0.5*k1, t+0.5*h)
    k3 = h*f(x+0.5*k2, t+0.5*h)
    k4 = h*f(x+k3, t+h)
    x = x + (k1 + 2*k2 + 2*k3 + k4)/6
    t = t + h
    T.append(t)
    X.append(x)

リストデータを plot.plot

Python でリスト形式の離散的な数値データをグラフに描く際，SymPy には適当な関数がないようなので，matplotlib.pyplot を import して使ってみます。

SymPy は確かに便利ですが，解析関数と数値データのグラフの重ね書きができないようで残念です（gnuplot でも Maxima でもできるのに）。できる！と知っている方がいましたら，ぜひお知らせください。

import matplotlib.pyplot as plt

plt.plot(T, X);

分割数 N を変えて結果をみます。分割数を変えても変わらない桁数までをもって，精度の良い数値解とします。

# まず，直前の結果を表示
print("N = ", N, "  ", X[N])

t0 = 0
x0 = 0.1
t1 = 10
N = 400

h = (t1 - t0)/N

t = t0
x = x0
# リスト T, X を初期化し，初期条件を append
T = []
X = []
T.append(t)
X.append(x)

for i in range(1, N+1):
    k1 = h*f(x, t)
    k2 = h*f(x+0.5*k1, t+0.5*h)
    k3 = h*f(x+0.5*k2, t+0.5*h)
    k4 = h*f(x+k3, t+h)
    x = x + (k1 + 2*k2 + 2*k3 + k4)/6
    t = t + h
    T.append(t)
    X.append(x)

print("N = ", N, "  ", X[N])

N =  100    0.9995915652218679
N =  400    0.9995915675088635

参考：dsolve で解析解

数値解析の本来の目的は，解析的に解けない問題を何とかして解きたい，解いたときの精度もちゃんと把握したい，ということです。

なので，上記の数値解がどの程度の精度で精確であるかを調べるために解析解と比較するのは本来は邪道なのですが，たまたま，上記の微分方程式は解析的に解けて，以下のようになります。

$$x(t) = \frac{\exp(t)}{9 + \exp(t)}$$

実は SymPy では dsolve() 関数を使って常微分方程式を解析的に解くことができます。

まず，解くべき微分方程式を eq として定義します。

t = Symbol('t')
x = Function('x')(t)
eq = Eq(Derivative(x,t), (1-x)*x)
eq

$\displaystyle \frac{d}{d t} x{\left(t \right)} = \left(1 - x{\left(t \right)}\right) x{\left(t \right)}$

初期条件（$x(t=0) = 0.1$）を ics に与えて，dsolve() 関数を使って微分方程式 eq を解きます。

a1 = dsolve(eq, x, ics={x.subs(t,0):1/10})
a1

$\displaystyle x{\left(t \right)} = \frac{1}{1 + 9.0 e^{- t}}$

解析的に解いた $x(t)$ に $t=10$ を入れて値を確かめます。

def xan(s):
    return a1.rhs.subs(t,s)
xan(10.0)

$\displaystyle 0.999591567517392$

SciPy による数値解法

SciPy には常微分方程式を数値的に解いてくれる関数 が用意されています。それを使ってみましょう。（すでに定義されている函数を使うほうが楽ちん。）

scipy.integrate.odeint()

import scipy.integrate  # SymPy のと混同しないようにフルネームで
import numpy as np
import matplotlib.pyplot as plt

def f_dxdt(x, t):
    return (1-x)*x

t_ini = 0
x_ini = 0.1
t_end = 10
N = 400

t_list = np.linspace(t_ini, t_end, N+1)
ans = scipy.integrate.odeint(f_dxdt, x_ini, t_list)

print("解析解 ", xan(10.0))
print("自作解 ", X[N])
print("odeint ", ans[N, 0])

解析解  0.999591567517392
自作解  0.9995915675088635
odeint  0.9995915695583809

あれれ？ scipy.integrate.odeint() の精度がイマイチですねぇ。手作りの4次のルンゲ・クッタ法の結果より精度が悪いようです。

scipy.integrate.solve_ivp()

SciPy のマニュアルによれば，odeint() は Old API だそうで，そうであれば，まずマニュアルのトップにある scipy.integrate.solve_ivp() を利用してみましょう。

# scipy.integrate.solve_ivp() は
# dy/dt = fun(t, y) の形を解く

def fun(t, y):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
# 連立でないときは，1次元リストとして返す
    return [(1-y)*y]

t_ini = 0
y_ini = [0.1]
t_end = 10
t_span = [t_ini, t_end]

t_list = np.linspace(t_ini, t_end, 21)

ansivp = scipy.integrate.solve_ivp(fun, t_span, y_ini, 
                    t_eval = t_list, 
                    rtol = 1.e-12, 
                    atol = 1.e-14)

# odeint では刻み幅と精度は連動しているが，
# solve_ivp では精度は rtol, atol で指定

# 計算結果は .y として取り出します。
ansivp.y

array([[0.1       , 0.1548281 , 0.23196932, 0.33242786, 0.45085306,
        0.57512085, 0.69056786, 0.78630171, 0.85848645, 0.90910664,
        0.94282562, 0.9645239 , 0.97817805, 0.98664969, 0.99185987,
        0.9950469 , 0.99698992, 0.99817213, 0.99889054, 0.99932679,
        0.99959157]])

print("解析解    ", xan(10.0))
print("自作解    ", X[N])
print("odeint    ", ans[N, 0])
print("solve_ivp ", ansivp.y[0,-1])

解析解     0.999591567517392
自作解     0.9995915675088635
odeint     0.9995915695583809
solve_ivp  0.9995915675171966

2階常微分方程式

次のような2階常微分方程式をルンゲ・クッタ法で解く。 $$\frac{d^2 x}{dt^2 } = f\left(x, \frac{dx}{dt}, t\right)$$

この場合には， $\displaystyle v \equiv \frac{dx}{dt}$ とおけば，次のような連立1階常微分方程式の形に帰着できる。 $$\begin{eqnarray} \frac{dx}{dt} &=& F_1(x, v, t) = v \\ \frac{dv}{dt} &=& F_2(x, v, t) = f(x, v, t) \end{eqnarray}$$

初期条件を $t=t_0$ で $x(t_0) = x_0, v(t_0) = v_0$ とし，$t = t_1$ までを $N$ 分割して解きます。

$$h = \frac{t_1 - t_0}{N}, \quad t_i = t_0 + i\, h, \quad x_i = x(t_i), \quad v_i = v(t_i)$$

とし，以下の式から次のステップの値 $x_{i+1}, v_{i+1}$ を決めます。

$$\begin{eqnarray} k_1 &=& h\, F_1(x_i, v_i, t_i) \\ m_1 &=& h\, F_2(x_i, v_i, t_i) \\ k_2 &=& h\, F_1(x_i + k_1/2, v_i + m_1/2, t_i + h/2) \\ m_2 &=& h\, F_2(x_i + k_1/2, v_i + m_1/2, t_i + h/2) \\ k_3 &=& h\, F_1(x_i + k_2/2, v_i + m_2/2, t_i + h/2) \\ m_3 &=& h\, F_2(x_i + k_2/2, v_i + m_2/2, t_i + h/2) \\ k_4 &=& h\, F_1(x_i + k_3, v_i + m_3, t_i + h) \\ m_4 &=& h\, F_2(x_i + k_3, v_i + m_3, t_i + h) \\ x_{i+1} &=& x_i + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4)\\ v_{i+1} &=& v_i + \frac{1}{6} (m_1 + 2 m_2 + 2 m_3 + m_4) \end{eqnarray}$$

例：減衰＋強制振動

例として，簡単な減衰＋強制振動の方程式

$$\frac{d^2 x}{dt^2} = -x - a \frac{dx}{dt} + b \cos(t)$$

を，

• $a , b$ にいくつか値を入れて
• 初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 3, v_0 = v(t_0) = 0$，
• $t = t_0 = 0$ から $t_1 = 20$ までを $N$ 分割

して解きます。

$N$ の値をいくつか変えてみて，精度を確認します。

a = 0.5
b = 0.2

def F1(x, v, t):
    return v
def F2(x, v, t):
    return -x - a*v + b*cos(t)

t0 = 0
x0 = 3.0
v0 = 0.0
t1 = 20

for j in range(0, 4):
    N = 100*2**j
    h = (t1 - t0)/N

    # 計算結果を配列に入れます。
    T = []
    X = []
    V = []
    t = t0; x = x0; v = v0
    T.append(t)
    X.append(x)
    V.append(v)
    
    for i in range(1, N+1):
        k1 = h*F1(x, v, t)
        m1 = h*F2(x, v, t)
        k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
        m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
        k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
        m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
        k4 = h*F1(x+k3, v+m3, t+h)
        m4 = h*F2(x+k3, v+m3, t+h)
        x = x + (k1 + 2*k2 + 2*k3 + k4)/6
        v = v + (m1 + 2*m2 + 2*m3 + m4)/6
        t = t + h 
        T.append(t)
        X.append(x)
        V.append(v)

    print("N = ", N, "  ", X[N]) 

N =  100    0.383972022819058
N =  200    0.383967288337950
N =  400    0.383966889383019
N =  800    0.383966861347239

$N = 400$ で小数点以下7桁くらいの精度が出ているようですので，$N=400$ として続けます。

まずは，$a = 0, b = 0$，すなわち単振動の場合の方程式を解きます。

a = 0
b = 0

def F1(x, v, t):
    return v
def F2(x, v, t):
    return -x - a*v + b*cos(t)

t0 = 0
x0 = 3.0
v0 = 0.0
t1 = 20

N = 400
h = (t1 - t0)/N

# 計算結果を配列に入れます。
T0 = []
X0 = []
V0 = []
t = t0; x = x0; v = v0
T0.append(t)
X0.append(x)
V0.append(v)

for i in range(1, N+1):
    k1 = h*F1(x, v, t)
    m1 = h*F2(x, v, t)
    k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
    m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
    k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
    m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
    k4 = h*F1(x+k3, v+m3, t+h)
    m4 = h*F2(x+k3, v+m3, t+h)
    x = x + (k1 + 2*k2 + 2*k3 + k4)/6
    v = v + (m1 + 2*m2 + 2*m3 + m4)/6
    t = t + h 
    T0.append(t)
    X0.append(x)
    V0.append(v)

plt.plot(T0, X0);

次に，$a = 0.5, b = 0.2$ の場合を解きます。

a = 0.5
b = 0.2

def F1(x, v, t):
    return v
def F2(x, v, t):
    return -x - a*v + b*cos(t)

t0 = 0
x0 = 3.0
v0 = 0.0
t1 = 20

N = 400
h = (t1 - t0)/N

# 計算結果を配列に入れます。
T1 = []
X1 = []
V1 = []
t = t0; x = x0; v = v0
T1.append(t)
X1.append(x)
V1.append(v)

for i in range(1, N+1):
    k1 = h*F1(x, v, t)
    m1 = h*F2(x, v, t)
    k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
    m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
    k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
    m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
    k4 = h*F1(x+k3, v+m3, t+h)
    m4 = h*F2(x+k3, v+m3, t+h)
    x = x + (k1 + 2*k2 + 2*k3 + k4)/6
    v = v + (m1 + 2*m2 + 2*m3 + m4)/6
    t = t + h 
    T1.append(t)
    X1.append(x)
    V1.append(v)

plt.plot(T1, X1);

複数のグラフを重ねて描画

以下のように，plt.plot() を続けて実行すると重ねて描画してくれます。

plt.plot(T0, X0);
plt.plot(T1, X1);

凡例やグリッドの表示

凡例をつけ，グリッドを描く例です。

plt.plot(T0, X0, label='$a = 0, b = 0$');
plt.plot(T1, X1, label='$a = 0.5, b = 0.2$');
plt.grid();
plt.legend();

参考：dsolve で解析解

SymPy には常微分方程式を解く dsolve() 関数があります。これを使って解いてみます。

まず，$a, b$ を変数とした微分方程式を返す関数 eq(a, b) を定義します。

t = Symbol('t')
x = Function('x')(t)
a, b = symbols('a b')

def eq(a, b):
    return Eq(Derivative(x, t, 2), 
              -x -a*Derivative(x, t) + b*cos(t))

eq(a, b)

$\displaystyle \frac{d^{2}}{d t^{2}} x{\left(t \right)} = - a \frac{d}{d t} x{\left(t \right)} + b \cos{\left(t \right)} - x{\left(t \right)}$

まずは $a = 0, b = 0$ の場合に，初期条件 $x(0) = 3, v(0) = \dot{x}(0) = 0$ のもとで dsolve() で解析的に解きます。

ans = dsolve(eq(0, 0), x, 
             ics = {x.subs(t, 0):3, diff(x, t).subs(t,0):0})
ans

$\displaystyle x{\left(t \right)} = 3 \cos{\left(t \right)}$

次に，$a = \frac{1}{2}, b = \frac{2}{10}$ の場合。以下のように Rational() 関数を使うと，実数値に変換しないで（約分された）分数のままで扱います。

5/10, Rational(5,10), 2/10, Rational(2,10)

(0.5, 1/2, 0.2, 1/5)

a = Rational(5,10)
b = Rational(2,10)

ans1 = dsolve(eq(a, b), x, 
             ics = {x.subs(t, 0):3, diff(x, t).subs(t,0):0})
ans1

$\displaystyle x{\left(t \right)} = \left(\frac{7 \sqrt{15} \sin{\left(\frac{\sqrt{15} t}{4} \right)}}{75} + 3 \cos{\left(\frac{\sqrt{15} t}{4} \right)}\right) e^{- \frac{t}{4}} + \frac{2 \sin{\left(t \right)}}{5}$

求まった解析解を関数 x1(t) として定義する例です。（ちょっとトリッキーですかね？）

解 ans1 の右辺 rhs の t に t1 を代入して関数 x1(t1) を定義します。

t, u = symbols('t u')

def x1(u):
    return ans1.rhs.subs(t, u)

解析解と数値解を重ねて描画

関数のグラフは SymPy の plot() で以下の通りに描画できます。

plot(x1(t), (t, 0, 20));

一方で，数値的に求めた解は，matplotlib の plt.plot() で描画しました。以下のようにしても，重ねて描画してくれません。

plot(x1(t), (t, 0, 20));
plt.plot(T1, X1);

次善の策として，関数解のほうも以下のように数値データにします。

import numpy as np

t1 = np.arange(0, 20, 0.1)
# 関数 x1(t) は SymPy の関数なので，リスト化するのが面倒
npx1 = []
for i in range(0, len(t1)):
    npx1.append(float(x1(t1[i])))

plt.plot(t1, npx1, label="解析解");
plt.plot(T1[0:N+1:20], X1[0:N+1:20], 'o', label="数値解");
plt.grid();
plt.legend();

念のために，$t=20$ での数値解と解析解を比べてみます。

print("数値解：",X1[N])
print("解析解：",x1(20).evalf())

数値解： 0.383966889383019
解析解： 0.383966859373692

SciPy による数値解法

scipy.integrate.solve_ivp()

SciPy には常微分方程式を数値的に解いてくれる関数 scipy.integrate.solve_ivp() が用意されています。この函数は，1階常微分方程式だけでなく，連立1階常微分方程式も，従って高階の常微分方程式も解いてくれるので，それを使ってみましょう。（すでに定義されている函数を使うほうが楽ちん。）

# scipy.integrate.solve_ivp() は
# dy/dt = Fun(t, y) の形を解く

def Fun(t, y, a, b):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
    x = y[0]
    v = y[1]
    F1 = v
    F2 = -x - a*v + b*np.cos(t)
    return [F1, F2]

t_ini = 0
y_ini = [3, 0]
t_end = 20
t_span = [t_ini, t_end]

t_list = np.linspace(t_ini, t_end, 21)

a = 0.5; b = 0.2
ansivp = scipy.integrate.solve_ivp(Fun, t_span, y_ini, 
                    t_eval = t_list, 
                    args=(a, b), 
                    rtol = 1.e-12, 
                    atol = 1.e-14)

# odeint では刻み幅と精度は連動しているが，
# solve_ivp では精度は rtol, atol で指定

print("数値解：  ",X1[N])
print("解析解：  ",x1(20).evalf())
print("solve_ivp ", ansivp.y[0, -1])

数値解：   0.383966889383019
解析解：   0.383966859373692
solve_ivp  0.3839668593734701