葛西 真寿 弘前大学大学院理工学研究科
Python で数値解析のために必要なプログラミングと，いくつかの応用例を示します。
ここでは主に SymPy を使った例を示し，SciPy の利用についても解説します。
Python の文末は，単に改行するか，; で終わって改行するかです。; で終わる場合は結果が表示されません。
1 + 2
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 より大きい")
$\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)
sum = 0
i = 1
while i <=5:
sum = sum + i
i = i + 1
print(sum)
$\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))
Sum()¶上の例で，sum(n) という関数を再帰的に定義しましたが，定義などしなくても，実は SymPy を import すると，数列の和を計算する関数 Sum() が使えます。
from sympy import *
from sympy.abc import *
Sum(k, [k, 1, 5])
あれ，計算してくれませんね。以下のように .doit() をつけてみます。
Sum(k, [k, 1, 5]).doit()
一般の n までの和も...
Sum(k, [k, 1, n]).doit()
Sum(k**2, [k, 1, n]).doit()
解析的な微分の定義は（前方差分の極限として） $$\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 を定義し，近似的な数値微分を（前方差分と後方差分の平均である中心差分として）以下のように定義しましょう。
# 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)
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)
シンプソン法で $\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}
$\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)
$\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))
integrate(), evalf()¶実は数値微分などしなくても，Python は sympy を import すると，関数の積分をしてくれます。
integrate(f(x), (x, 0, pi/2))
解析的に積分できない場合は，数値積分によって近似値を求めることができます。
$\displaystyle \int_0^{\pi} \sin(\sin x)\, dx$ の計算をします。解析的に解けないと，そのままで返します。
integrate(sin(sin(x)), (x, 0, pi))
直前の結果を表す _ に .evalf() すると数値解を返します。
_.evalf()
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)
$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));
$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)
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)
以下のように，再帰定義により $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))
解の収束条件 h を変えても変わらない桁までが，精度のよい数値解ということになります。上の例では，
2.40482555769576
くらいでしょうか。
besselj(0, 2.40482555769576).evalf()
nsolve¶実は Python は sympy を import すると，関数の積分をしてくれる関数 nsolve() が使えます。
x = Symbol('x')
nsolve(besselj(0,x), x, 2.5, prec=15)
besselj(0,_).evalf()
4次のルンゲ・クッタ法で常微分方程式を解きます。
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])
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
初期条件（$x(t=0) = 0.1$）を ics に与えて，dsolve() 関数を使って微分方程式 eq を解きます。
a1 = dsolve(eq, x, ics={x.subs(t,0):1/10})
a1
解析的に解いた $x(t)$ に $t=10$ を入れて値を確かめます。
def xan(s):
return a1.rhs.subs(t,s)
xan(10.0)
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])
あれれ？
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
print("解析解 ", xan(10.0))
print("自作解 ", X[N])
print("odeint ", ans[N, 0])
print("solve_ivp ", ansivp.y[0,-1])
次のような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) $$を，
して解きます。
$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 = 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)
まずは $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
次に，$a = \frac{1}{2}, b = \frac{2}{10}$ の場合。以下のように Rational() 関数を使うと，実数値に変換しないで（約分された）分数のままで扱います。
5/10, Rational(5,10), 2/10, Rational(2,10)
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
求まった解析解を関数 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())
# 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])