Python で数値解析のために必要なプログラミングと，いくつかの例を示します。
ここでは数値解析ライブラリである SciPy を使わず，計算機代数システムである SymPy をあえて使った例を示します。

また SymPy 自体でもグラフ作成機能がありますが，ここでは，SymPy Plotting Backends (SPB) を使ってグラフを描いてみます。

ライブラリの import

In [1]:

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

# グラフは SymPy Plotting Backends (SPB) で描く
from spb import *

# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']

# デフォルト設定のため
import matplotlib.pyplot as plt
# mathtext font の設定
plt.rcParams['mathtext.fontset'] = 'cm'

SymPy による数値微分

解析的な微分の定義は（前方差分の極限として）

$\frac{d f}{d x} \equiv lim_{h \to 0} \frac{f (x + h) - f (x)}{h}$

でした。滑らかな関数であれば，（後方差分の極限として）
$\frac{d f}{d x} = lim_{h \to 0} \frac{f (x) - f (x - h)}{h}$
としても良いです。

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

$\begin{array}{rcl} \frac{d f}{d x} & ≃ & \frac{1}{2} {\frac{f (x + h) - f (x)}{h} + \frac{f (x) - f (x - h)}{h}} \\ = & \frac{f (x + h) - f (x - h)}{2 h} \end{array}$

関数 $f (x)$ が初等関数（および初等関数の合成関数）である場合は解析的に微分できますが，そうでない場合，例えば $f (x)$ が数値的にしか与えられていないとか，何か解析的に微分できない場合には，数値微分によって近似的な値を求める必要があります。

また，解析的に微分できそうだが，人力ではけっこう大変そうな場合にも，数値微分によって近似的な値を使うほうが便利な場合があるかもしれません。

以下のように，関数 f(x) を数値微分してくれる関数 ndiff() を定義します。

In [2]:

def ndiff(f, x, h=1e-5):
    # h を省略すると h=1e-5 として計算する
    return (f(x + h) - f(x - h))/(2*h)

例として， $f (x) = \sin x$ の $x = 3$ における数値微分を求めてみます。（答えは $\cos 3$ です。）

In [3]:

def f(x):
    return sin(x)

for h in [1e-4, 1e-5, 1e-6, 1e-7, 1e-8]:
    print('h = %3.1e, ndiff = %.15f' % (h, ndiff(f, 3, h)))

print(' '*10 + 'cos(%3.1f) = %.15f' % (3, cos(3)))

h = 1.0e-04, ndiff = -0.989992494952602
h = 1.0e-05, ndiff = -0.989992496590319
h = 1.0e-06, ndiff = -0.989992496730485
h = 1.0e-07, ndiff = -0.989992494926373
h = 1.0e-08, ndiff = -0.989992490763036
          cos(3.0) = -0.989992496600445

以上の結果からわかるように，h をやみくもに小さくすればするほど，数値微分の精度があがる，というわけではないです。

○練習：関数の数値微分

$f (x) = \frac{x^{3}}{e^{x} - 1}$ の極大値を求めるための準備として，まず， $f (x)$ とその数値微分 ndiff(f, x) を $0 \leq x \leq 5$ の範囲でグラフにする。

（ $f (x)$ の微分は簡単に求められるのだが，ここはせっかく覚えたので数値微分してみます。）

In [4]:

def f(x):
    return x**3/(exp(x) - 1)

graphics(
    line(f(x), (x, 0.001, 5), '$f(x)$'),
    line(ndiff(f, x), (x, 0.001, 5), '$f^{\prime}(x)$'), 
    xlabel = '$x$', ylabel = '　'
);

上のグラフから，数値微分した微分係数 $f^{'} (x)$ の値がゼロになる $x$ は， $2 < x < 3$ のあたりであること，そのあたりで確かに $f (x)$ が極大（最大）となっていることが見てとれます。

念のために，解析的に微分した結果と比較します。目視ではほぼ区別がつかないですね。

In [5]:

df = diff(f(x), x)

graphics(
    line(ndiff(f, x), (x, 0.0001, 5), '数値微分',
         rendering_kw={'lw':2}),
    line(df, (x, 0.0001, 5), '解析的微分',
         rendering_kw={'lw':0.5})
);

SymPy による数値積分：`integrate().evalf()`

$\int_{0}^{π} \sin (\sin x) d x$ の計算をします。解析的に解けないと，SymPy はそのままの形で返します。

In [6]:

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

Out[6]:

\int_{0}^{π} \sin (\sin (x)) d x

.evalf() すると数値積分した値を返します。（SymPy が数値積分してくれるので，わざわざ自分で定義する必要はありません。）

In [7]:

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

Out[7]:

1.78648748195005

.evalf() のかわりに .n() を使ってもよいようです。

In [8]:

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

Out[8]:

1.78648748195005

以下のように N() を使ってもよいです。

In [9]:

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

Out[9]:

1.78648748195005

参考：シンプソン法による数値積分

シンプソン法で $\int_{a}^{b} f (x) d x$ を求める。

積分区間 $[a, b]$ を $N (= 2 n)$ 等分（偶数等分）し，

$h = \frac{b - a}{N}, f_{i} = f (a + i h)$

とする。シンプソン法の公式は

$\begin{array}{rcl} \int_{a}^{b} f (x) d x & ≃ & \frac{h}{3} (f_{0} + 4 f_{1} + f_{2}) + \frac{h}{3} (f_{2} + 4 f_{3} + f_{4}) + \\ \frac{h}{3} (f_{4} + 4 f_{5} + f_{6}) + \frac{h}{3} (f_{6} + 4 f_{7} + f_{8}) + \dots \\ + \frac{h}{3} (f_{N - 2} + 4 f_{N - 1} + f_{N}) \\ = & \frac{h}{3} (f_{0} + 2 \sum_{i = 1}^{n - 1} f_{2 i} + 4 \sum_{i = 1}^{n} f_{2 i - 1} + f_{2 n}) \end{array}$

SymPy では .evalf() を使って数値積分することができるのでわざわざ自作の関数をつくる必要はないのですが，教育的見地からあえて，シンプソン法の公式を使って数値積分の値を返す関数 simpson() を以下のように定義してみます。

In [10]:

def simpson(f, a, b, N):
    h = float((b - a)/N)
    n = Rational(N, 2)
    tmp = h/3 * (f(a) + 
                 2*summation(f(a + 2*i * h), (i, 1, n-1)) +
                 4*summation(f(a + (2*i-1) * h), (i, 1, n)) +
                 f(b)
                 )
    return tmp

N の値を変えて $\int_{0}^{π} \sin (\sin x) d x$ を simpson() で数値積分して，精度を確認してみます。（N が大きい場合は時間がかかります。）

In [11]:

def f(x):
    return sin(sin(x))

In [12]:

for N in [10*2**i for i in range(9)]:
    print('N = %4d, %.15f' %(N, simpson(f, 0, pi, N)))

ninteg = integrate(f(x), (x, 0, pi)).evalf()
print('.evalf()  %.15f' % ninteg)

N =   10, 1.786721213744078
N =   20, 1.786501256275661
N =   40, 1.786488331266679
N =   80, 1.786487534856180
N =  160, 1.786487485253951
N =  320, 1.786487482156504
N =  640, 1.786487481962953
N = 1280, 1.786487481950860
N = 2560, 1.786487481950103
.evalf()  1.786487481950052

以上の結果から，どんなことがわかるかな？

SymPy による方程式の数値解：`nsolve()`

SymPy には数値解を求める関数 nsolve() が用意されていますので，これを使ってみます。（わざわざ自分で定義する必要はありません。）

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

以下では例としてベッセル関数 $J_{0} (x) =$ besselj(0,x) のゼロ点を求めてみます。ベッセル関数 besselj() は

from sympy import *

すれば使えます。 $x$ 軸を横切る点がわかりやすいように， $y = 0$ （ $x$ 軸）も描いてみます。

In [13]:

graphics(
    line(besselj(0, x), (x, 0, 10)), 
    line(0, (x, 0, 10), ''), # x 軸
    xlabel = "$x$", ylabel = "$J_0(x)$", 
    xlim = [0, 10], 
    title = "ベッセル関数 $J_0(x)$"
);

方程式の数値解を求めてくれる関数は nsolve() です。 $x = 2$ の近くで $x$ 軸と交わっていますから，このあたりの解を探します。

In [14]:

# besselj(0,x) = 0 となる解を
# x が 2 のあたりで探す

nsolve(besselj(0,x), x, 2)

Out[14]:

2.40482555769577

[2, 4] のように範囲を指定して解を探すこともできます。

In [15]:

# 2 <= x <= 4 の間で
# besselj(0, x) = 0 となる x を探す

nsolve(besselj(0,x), x, [2, 4])

Out[15]:

2.40482555769577

In [16]:

# ほんとに解になってるか確認。ゼロに近ければよい。
besselj(0,_).evalf()

Out[16]:

- 6.10876525973673 \cdot 10^{- 17}

○練習：関数の極大値

$f (x) = \frac{x^{3}}{e^{x} - 1}$ が極大値をとるときの $x$ の値，およびそのときの極大値を求める。
$g (x) = \frac{x^{5}}{e^{x} - 1}$ が極大値をとるときの $x$ の値，およびそのときの極大値も同様にして求めてみよ。

まず $f (x)$ の数値微分 ndiff(f, x) をグラフにしてみる。

In [17]:

def f(x):
    return x**3/(exp(x) - 1)

graphics(
    line(ndiff(f, x), (x, 0, 5), '$f(x)$ の数値微分'), 
    line(0, (x, 0, 5), ''), 
    xlabel='$x$', ylabel='$y$'
);

ndiff(f, x) は $2 < x < 3$ のあたりでゼロになっているので，そのあたりを nsolve() で探してみる。

In [18]:

nsolve(ndiff(f, x), x, [2, 3])

Out[18]:

2.82143937214806

In [19]:

# そのときの f(x) の極大値は...
f(_)

Out[19]:

1.42143547274774

念のため，数値微分 ndiff(f, x) ではなく，解析的な微分 df $= \frac{d f}{d x}$ がゼロとなる点も nsolve() で探してみる。

In [20]:

def f(x):
    return x**3/(exp(x) - 1)

df = diff(f(x), x)
display(df)

nsolve(df, x, [2, 3])

- \frac{x^{3} e^{x}}{{(e^{x} - 1)}^{2}} + \frac{3 x^{2}}{e^{x} - 1}

Out[20]:

2.82143937212208

In [21]:

# そのときの f(x) の極大値は...
f(_)

Out[21]:

1.42143547274774

ルンゲ・クッタ法による1階常微分方程式の数値解法

常微分方程式の数値解法をしてくれる関数は SymPy には用意されていないようです。備え付けの関数に頼らず，1階常微分方程式

$\frac{d x}{d t} = f (t, x)$

を4次のルンゲ・クッタ法で解いてみます。

初期条件を $t = t_{0}$ で $x (t_{0}) = x_{0}$ とし， $t = t_{1}$ までを $N$ 分割して解きます。

$h = \frac{t_{1} - t_{0}}{N}, t_{i} = t_{0} + i h, x_{i} = x (t_{i})$

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

$\begin{array}{rcl} k_{1} & = & h f (t_{i}, x_{i}) \\ k_{2} & = & h f (t_{i} + h / 2, x_{i} + k_{1} / 2) \\ k_{3} & = & h f (t_{i} + h / 2, x_{i} + k_{2} / 2) \\ k_{4} & = & h f (t_{i} + h, x_{i} + k_{3}) \\ x_{i + 1} & = & x_{i} + \frac{1}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) \end{array}$

以上のような計算を行う関数 rk1() を Python の関数として以下のようにして定義します。

In [22]:

def rk1(func, x0, t0, t1, N):
    """ dx/dt = func(t, x) を
        初期条件 t = t0 で x = x0 とし，
        [t0, t1] を N 分割して
        4次のルンゲ・クッタ法で解く
    """
    h = (t1 - t0)/N
    T = [t0]
    X = [x0]
    t = t0
    x = x0
    for i in range(N):
        k1 = h*func(t, x)
        k2 = h*func(t + h/2, x + k1/2)
        k3 = h*func(t + h/2, x + k2/2)
        k4 = h*func(t + h, x + k3)
        x = x + (k1 + 2*k2 + 2*k3 + k4)/6
        t = t + h
        T.append(t)
        X.append(x)
    return [T, X]

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

$\frac{d x}{d t} = f (t, x) = (1 - x) x$

を，初期条件 $t_{0} = 0$ で $x_{0} = x (t_{0}) = 0.1$ ， $t = t_{0} = 0$ から $t_{1} = 10$ までを $N$ 分割して解きます。

ルンゲ・クッタ法の精度は分割数 $N$ に依存します。 $N$ の値を変えて $x (t_{1})$ の値を計算して精度を確認します。

In [23]:

# dx/dt = f(t, x) を解く。
def f(t, x):
    return (1-x)*x

x0 = 0.1
t0 = 0
t1 = 10

for N in [100, 200, 400, 800, 1600]:
    T1, X1 = rk1(f, x0, t0, t1, N)
    print("N =%5d, x(t1) =%15.12f" % (N, X1[-1]))

N =  100, x(t1) = 0.999591565222
N =  200, x(t1) = 0.999591567379
N =  400, x(t1) = 0.999591567509
N =  800, x(t1) = 0.999591567517
N = 1600, x(t1) = 0.999591567517

上の結果から， $N = 800$ くらいで小数点以下 $12$ 桁程度の精度があることがわかります。ではこのくらいの $N$ で計算をすることにします。

In [24]:

# dx/dt = f(t, x) を解く。
def f(t, x):
    return (1-x)*x

x0 = 0.1
t0 = 0
t1 = 10

ndiv = 40
Ndiv = 20
N = ndiv * Ndiv

T1, X1 = rk1(f, x0, t0, t1, N)

T1，X1 には N+1 点の数値解の値が代入されています。滑らかな曲線のグラフを描くには，全部の点を使えばいいです。

In [25]:

graphics(
    list_2d(T1, X1), 
    xlabel = "$t$", ylabel = "$x(t)$",
    title = "ロジスティック方程式の数値解"
);

ここでは少し間引いて Ndiv+1 個の点をグラフにしてみます。（あとで解析解と比べてみるため。）

In [26]:

# データを間引く
t1 = [T1[i] for i in range(0, len(T1), ndiv)]
x1 = [X1[i] for i in range(0, len(X1), ndiv)]

print('元の要素数は',len(X1))
print('間引くと　　',len(x1))

元の要素数は 801
間引くと　　 21

In [27]:

# データを間引いてグラフに

graphics(
    list_2d(t1, x1, scatter = True), 
    xlabel = "$t$", ylabel = "$x(t)$",
    title = "ロジスティック方程式の数値解"
);

参考：SymPy による1階常微分方程式の解析解

簡単な常微分方程式の場合は SymPy の dsolve() で解析的に解くことができます。

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

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

を，初期条件 $t = 0$ で $x (0) = 0.1$ として解きます。

In [28]:

# 解くべき微分方程式を eq として定義
var('t')
x = Function('x')(t)
eq = Eq(Derivative(x, t), (1-x)*x)
eq

Out[28]:

\frac{d}{d t} x (t) = (1 - x (t)) x (t)

In [29]:

# 初期条件を特に設定せずに微分方程式を解く
dsolve(eq, x)

Out[29]:

x (t) = \frac{1}{C_{1} e^{- t} + 1}

In [30]:

# 初期条件を t=0 で x(0) = 1/10 として解く
sol1 = dsolve(eq, x, ics={x.subs(t,0): Rational(1,10)})
sol1

Out[30]:

x (t) = \frac{1}{1 + 9 e^{- t}}

In [31]:

# ans の右辺を表示 
# rhs は rirhgt-hand-side
simplify(sol1.rhs)

Out[31]:

\frac{e^{t}}{e^{t} + 9}

… ということで，解析解は

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

であることがわかります。
得られた解析解 $x (t)$ とルンゲ・クッタ法で得られた数値解をグラフにする例。

In [32]:

graphics(
    line(sol1.rhs, (t, 0, 10), '解析解'), 
    list_2d(t1, x1, '数値解', scatter = True, 
            # "ms" つまり "markersize" を小さめに
            rendering_kw={"ms":3}), 
    xlabel = "$t$", ylabel = "$x(t)$",
    title = "ロジスティック方程式の解"
);

得られた $x (t)$ をどうやって関数に定義するか悩みます。以下の例でなんとかなりますが，少しトリッキーですかねぇ。

In [33]:

def xsol1(s):
    return sol1.rhs.subs(t,s)

xsol1(t)

Out[33]:

\frac{1}{1 + 9 e^{- t}}

ルンゲ・クッタ法による2階常微分方程式の数値解法

常微分方程式の数値解法をしてくれる関数は SymPy には用意されていないようです。備え付けの関数に頼らず，次のような2階常微分方程式をルンゲ・クッタ法で解いてみます。

$\frac{d^{2} x}{d t^{2}} = F (t, x, \frac{d x}{d t})$

この場合には， $v \equiv \frac{d x}{d t}$ とおけば，次のような連立1階常微分方程式の形に帰着できる。

$\begin{array}{rcl} \frac{d x}{d t} & = & F_{1} (t, x, v) = v \\ \frac{d v}{d t} & = & F_{2} (t, x, v) = F (t, x, v) \end{array}$

初期条件を $t = t_{0}$ で $x (t_{0}) = x_{0}, v (t_{0}) = v_{0}$ とし， $t = t_{1}$ までを $N$ 分割して解きます。

$h = \frac{t_{1} - t_{0}}{N}, t_{i} = t_{0} + i h, x_{i} = x (t_{i}), v_{i} = v (t_{i})$

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

$\begin{array}{rcl} k_{1} & = & h F_{1} (t_{i}, x_{i}, v_{i}) = h v_{i} \\ m_{1} & = & h F_{2} (t_{i}, x_{i}, v_{i}) \\ k_{2} & = & h F_{1} (t_{i} + h / 2, x_{i} + k_{1} / 2, v_{i} + m_{1} / 2) = h (v_{i} + m_{1} / 2) \\ m_{2} & = & h F_{2} (t_{i} + h / 2, x_{i} + k_{1} / 2, v_{i} + m_{1} / 2) \\ k_{3} & = & h F_{1} (t_{i} + h / 2, x_{i} + k_{2} / 2, v_{i} + m_{2} / 2) = h (v_{i} + m_{2} / 2) \\ m_{3} & = & h F_{2} (t_{i} + h / 2, x_{i} + k_{2} / 2, v_{i} + m_{2} / 2) \\ k_{4} & = & h F_{1} (t_{i} + h, x_{i} + k_{3}, v_{i} + m_{3}) = h (v_{i} + m_{3}) \\ m_{4} & = & h F_{2} (t_{i} + h, x_{i} + k_{3}, v_{i} + m_{3}) \\ 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{array}$

以上のような計算を行う関数 rk2() を Python の関数として以下のようにして定義します。

In [34]:

def rk2(Func, x0, v0, t0, t1, N):
    """ d^2 x/dt^2 = Func(t, x, v) を
        初期条件 t = t0 で x = x0, v = v0 とし，
        [t0, t1] を N 分割して
        4次のルンゲ・クッタ法で解く
    """
    h = float((t1 - t0)/N)
    T = [t0]
    X = [x0]
    V = [v0]
    t = t0
    x = x0
    v = v0
    for i in range(N):
        k1 = h*v
        m1 = h*Func(t, x, v)
        k2 = h*(v+m1/2)
        m2 = h*Func(t+h/2, x+k1/2, v+m1/2)
        k3 = h*(v+m2/2)
        m3 = h*Func(t+h/2, x+k2/2, v+m2/2)
        k4 = h*(v+m3)
        m4 = h*Func(t+h, x+k3, v+m3)
        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)
    return [T, X, V]

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

$\frac{d^{2} x}{d t^{2}} = F (t, x, v) = - x - a v + b \cos (t), v = \frac{d x}{d 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$ の値をいくつか変えてみて，精度を確認します。（結構時間がかかります。）

In [35]:

def F(t, x, v):
    global a, b
    return -x - a*v + b*cos(t)

a = 0.5
b = 0.2
x0 = 3
v0 = 0
t0 = 0
t1 = 20

for N in [200, 400, 800, 1600]:
    T2, X2, V2 = rk2(F, x0, v0, t0, t1, N)
    # X2 の最後の要素 X2[-1] を表示
    print("N =%5d, x(t1) =%15.12f" % (N, X2[-1]))

N =  200, x(t1) = 0.383967288338
N =  400, x(t1) = 0.383966889383
N =  800, x(t1) = 0.383966861347
N = 1600, x(t1) = 0.383966859500

上の結果から， $N = 800$ くらいで小数点以下 $8$ 桁程度の精度があることがわかります。ではこのくらいの $N$ で計算をすることにします。

In [36]:

ndiv = 20
Ndiv = 40
N = ndiv * Ndiv

T2, X2, V2 = rk2(F, x0, v0, t0, t1, N)

T2, X2, V2 には N+1 点の数値解の値が代入されています。滑らかな曲線のグラフを描くには，全部の点を使えばいいです。

In [37]:

graphics(
    list_2d(T2, X2)
);

ここでは少し間引いてグラフにしてみます。（あとで解析解と比べてみるため。）

In [38]:

# データを間引く
t2 = [T2[i] for i in range(0, len(T2), ndiv)]
x2 = [X2[i] for i in range(0, len(X2), ndiv)]

print('元の要素数は    ',len(X2))
print('%2d 飛びに間引くと %2d' % (ndiv, len(x2)))

元の要素数は     801
20 飛びに間引くと 41

In [39]:

graphics(
    list_2d(t2, x2, scatter = True, 
            # "ms" つまり "markersize" を小さめに
            rendering_kw={"ms":3}), 
    xlabel = "$t$", ylabel = "$x(t)$",
    title = "減衰＋強制振動"
);

参考：SymPy による2階常微分方程式の解析解

簡単な常微分方程式の場合は SymPy の dsolve() で解析的に解くことができます。
例として，簡単な減衰＋強制振動の方程式

$\frac{d^{2} x}{d t^{2}} = - x - a \frac{d x}{d t} + 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) を定義します。

In [40]:

var('t a b')
x = Function('x')(t)

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

eq(a, b)

Out[40]:

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

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

In [41]:

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

Out[41]:

x (t) = 3 \cos (t)

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

In [42]:

a = Rational(1,2)
b = Rational(1,5)
sol2 = dsolve(eq(a, b), x, 
              ics = {x.subs(t, 0):3, diff(x, t).subs(t,0):0})
sol2

Out[42]:

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

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

解 sol2 の右辺 rhs の t に s を代入して関数 xsol2(s) を定義します。

In [43]:

def xsol2(s):
    return sol2.rhs.subs(t, s)

In [44]:

key1 = "$a = %3.1f, b = %3.1f$ の解析解" % (a, b)
key2 = "$a = %3.1f, b = %3.1f$ の数値解" % (a, b)

graphics(
    line(xsol2(t), (t, 0, 20), key1), 
    list_2d(t2, x2, key2, scatter = True, 
            rendering_kw={"ms":3}), 
    xlim = (0, 20), 
    xlabel = "$t$", ylabel = "$x(t)$",
    title = "減衰＋強制振動"
);

○練習：減衰＋強制振動

減衰＋強制振動の方程式を以下の場合に解き，グラフにする。

まず $a$ の影響を調べるため，

$a = 0, b = 0$
$a = 0.2, b = 0$
$a = 0.5, b = 0$
$a = 0.8, b = 0$

次に， $b$ の影響を調べるため，

$a = 0, b = 0$
$a = 0, b = 0.2$
$a = 0, b = 0.5$
$a = 0, b = 0.8$

○練習：単振り子

2階常微分方程式

$\frac{d^{2} θ}{d t^{2}} = - 4 π^{2} \sin θ$

を初期条件

$θ = 45^{\circ}, \frac{d θ}{d t} = 0 at t = 0$

のもとで $t = 0$ から $t = 3$ まで数値的に解け。（解析的には解けませんよ。）

$θ \Rightarrow x, \frac{d θ}{d t} \Rightarrow v$ と置き換えて，rk2() を使う。初期条件の角度はラジアンにして $\sin θ$ に渡すように。

In [45]:

# 変数 t, x, v の初期化
var('t x v')

def F(t, x, v):
    # pi のままだと厳密に計算しようとして
    # ものすごく時間がかかってしまうので
    return -4*float(pi)**2*sin(x)

x0 = float(rad(45))
v0 = 0
t0 = 0
t1 = 3

for N in [100, 200, 400, 800, 1600, 3200]:
    T2, X2, V2 = rk2(F, x0, v0, t0, t1, N)
    # X2 の最後の要素 X2[-1] を表示
    print("N =%5d, x(t1) =%15.12f" % (N, X2[-1]))

N =  100, x(t1) = 0.591465894161
N =  200, x(t1) = 0.591544091180
N =  400, x(t1) = 0.591548830281
N =  800, x(t1) = 0.591549121191
N = 1600, x(t1) = 0.591549139198
N = 3200, x(t1) = 0.591549140318

上の結果から，N = 1600 で小数点以下 7~8 桁程度の精度があることがわかる。

以下のように結構時間がかかるが，これで行こう。

セルの最初に %%time と書くことで，このセルを実行するのにかかった時間がわかる。

In [46]:

%%time

T2, X2, V2 = rk2(F, x0, v0, t0, t1, 1600)

CPU times: user 6.91 s, sys: 0 ns, total: 6.91 s
Wall time: 6.91 s

一方， $θ$ があまり大きくないとすれば，テイラー展開の1次の近似で $\sin θ ≃ θ$ であるから

$\frac{d^{2} θ}{d t^{2}} ≃ - 4 π^{2} θ$

$≃$ を $=$ に変えて方程式にすると，これは単振動の方程式であるから，同じ初期条件で以下のように解析的に解ける。

$θ_{s} = θ_{0} \cos (2 π t) = \frac{π}{4} \cos (2 π t)$

この単振動の解も重ねてグラフにしてみよ。周期にはどのような違いがあるか。

In [47]:

# 縦軸のリスト X2 の各要素をラジアンから度へ変換
XX2 = [deg(1)*X2[i] for i in range(len(X2))]

graphics(
    list_2d(T2, XX2, '数値解'), 
    line(deg(pi/4) * cos(2*pi*t), (t, 0, 3), '単振動解'), 
    xlabel = '$t$', ylabel = r'$\theta$ (度)'
);

以上のように，SymPy では解析的に計算（微分・積分）できるところは解析的にやって，できないところは数値的に解く，ということができて便利なのではあるが，数値解析用の関数は一部しか準備されておらず，自前で数値解析のプログラミングをしてみると，実行にかなり時間がかかるという残念なところもある。

もう少し Python で効率的な数値解析をしたい場合には，別途説明するライブラリ SciPy などを使うことになる。

SymPy で（あえて）数値解析

ライブラリの import

SymPy による数値微分

○練習：関数の数値微分

SymPy による数値積分：`integrate().evalf()`

参考：シンプソン法による数値積分

SymPy による方程式の数値解：`nsolve()`

○練習：関数の極大値

ルンゲ・クッタ法による1階常微分方程式の数値解法

参考：SymPy による1階常微分方程式の解析解

ルンゲ・クッタ法による2階常微分方程式の数値解法

参考：SymPy による2階常微分方程式の解析解

○練習：減衰＋強制振動

○練習：単振り子

このページの目次

このセクションの構成

SymPy で（あえて）数値解析

ライブラリの import

SymPy による数値微分

○練習：関数の数値微分

SymPy による数値積分：integrate().evalf()

参考：シンプソン法による数値積分

SymPy による方程式の数値解：nsolve()

○練習：関数の極大値

ルンゲ・クッタ法による1階常微分方程式の数値解法

参考：SymPy による1階常微分方程式の解析解

ルンゲ・クッタ法による2階常微分方程式の数値解法

参考：SymPy による2階常微分方程式の解析解

○練習：減衰＋強制振動

○練習：単振り子

このページの目次

このセクションの構成

SymPy による数値積分：`integrate().evalf()`

SymPy による方程式の数値解：`nsolve()`