Python で数値解析のために必要なプログラミングと，いくつかの例を示します。
ここでは数値解析ライブラリである SciPy を使わず，計算機代数システムである SymPy をあえて使った例を示します。なお，SymPy には常微分方程式の数値解法の機能はないようなので，dsolve() で解析的に解ける例のみを紹介します。
# SymPy を使うときのおまじない
from sympy import *
from sympy.abc import *
from sympy import pi
# どうしても matplotlib が必要な場合
import matplotlib.pyplot as plt
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
SymPy による数値微分
解析的な微分の定義は（前方差分の極限として）
$$\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}
例として $f(x) = \sin x$ の $x = 3$ における数値微分を求めてみます。
（答えは $\cos(3)$ です。）
SymPy では以下のように解析的に微分できてしまいますが，ここはあえて数値的に微分するのです。
diff(sin(x), x)
def f(x):
return sin(x)
def diff_f(x, h):
return (f(x+h)-f(x-h))/(2.0*h)
# 数値微分の値
print("数値微分", diff_f(3, 1E-6))
# 解析的微分の値
print("解析解 ", float(cos(3)))
以下の例では，$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*h)
p1 = plot(cos(x) - diff_f(x, 1E-6), (x, -2*pi, 2*pi));
# Web ページ用に画像ファイル保存
p1.save('./sk-sympy1.svg');
以下の例からわかるように，h をやみくもに小さくすればするほど，数値微分の精度があがる，というわけでもなさそうです。（縦軸の目盛りに注目。）
p2 = plot(cos(x) - diff_f(x, 1E-7), (x, -2*pi, 2*pi));
# Web ページ用に画像ファイル保存
p2.save('./sk-sympy2.svg');
SymPy による数値積分
$\displaystyle \int_0^{\pi} \sin(\sin x)\, dx$ の計算をします。解析的に解けないと，そのままで返します。
integrate(sin(sin(x)), (x, 0, pi))
.evalf() すると数値積分した値を返します。
integrate(sin(sin(x)), (x, 0, pi)).evalf()
SymPy による方程式の数値解
まず，解を求める関数 $f(x)$ を定義し，$f(x) = 0$ となる $x$ のあたりをつけるためにグラフを描いてみます。
以下では例としてベッセル関数 besselj(0,x) のゼロ点を求めています。ベッセル関数 besselj() は from sympy import * すれば使えます。
p3 = plot(besselj(0,x), (x, 0, 10),
title = "ベッセル関数 J0(x)");
# Web ページ用に画像ファイル保存
p3.save('./sk-sympy3.svg');
方程式の数値解を求めてくれる関数は nsolve() です。$x=2.5$ のあたりで $x$ 軸と交わっていますから，このあたりの解を探します。
# besselj(0,x) = 0 となる解を
# x が 2.5 のあたりから始めて 15 桁の精度で求める
nsolve(besselj(0,x), x, 2.5, prec=15)
besselj(0,_).evalf()
SymPy による1階常微分方程式の解法
常微分方程式の数値解法をしてくれる関数は SymPy には用意されていないようです。簡単な場合は dsolve() で解析的に解くことができます。
例として，簡単なロジスティック方程式
$$ \frac{dx}{dt} = (1-x)\,x$$
を，初期条件 $t = 0$ で $x(0) = 0.1$ として解きます。
# 解くべき微分方程式を eq として定義
t = Symbol('t')
x = Function('x')(t)
eq = Eq(Derivative(x,t), (1-x)*x)
eq
ans = dsolve(eq, x, ics={x.subs(t,0): Rational(1,10)})
ans
得られた $x(t)$ をグラフにする際，どうやって答えを関数に定義するか悩みます。以下の例でなんとかなりますが，少しトリッキーですかねぇ。
def xans(s):
return ans.rhs.subs(t,s)
p4 = plot(xans(t), (t, 0, 10),
axis_center = (0,0),
xlim = (0, 10), ylim = (0, 1),
title = "ロジスティック方程式の解析解");
# Web ページ用に画像ファイル保存
p4.save('./sk-sympy4.svg');
SymPy による2階常微分方程式の解法
常微分方程式の数値解法をしてくれる関数は SymPy には用意されていないようです。簡単な場合は dsolve() で解析的に解くことができます。
例として，簡単な減衰＋強制振動の方程式
$$ \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$，
として解きます。
まず，$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() で解析的に解きます。
ans0 = dsolve(eq(0, 0), x,
ics = {x.subs(t, 0):3, diff(x, t).subs(t,0):0})
ans0
次に，$a = \frac{1}{2}, b = \frac{2}{10}$ の場合。以下のように Rational() 関数を使うと，実数値に変換しないで（約分された）分数のままで扱います。
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 に u を代入して関数 x1(u) を定義します。
def x1(u):
return ans1.rhs.subs(t, u)
p5 = plot(x1(t), (t, 0, 20),
axis_center = (0,0),
xlim = (0, 20), ylim = (-2, 3),
title = "減衰＋強制振動");
# Web ページ用に画像ファイル保存
p5.save('./sk-sympy5.svg');
