計算機基礎B - 10. グラフの可視化
Table of Contents
1 ルジャンドル多項式
1.1 定義
n 次のルジャンドル多項式 Pn(x) は区間 [-1, 1] の実数 x に対して次の漸化式で定義されます。
P0(x) = 1 P1(x) = x Pn(x) = ( (2*n-1) * x * Pn-1(x) - (n-1) * Pn-2(x) ) / n (n >= 2)
Pn(x) を計算する関数 lp(n, x) は次のように書けます。
def lp(n, x): a = 1.0 b = x for i in range(2, n+1): # i = 2, 3, 4, ..., n c = ((2*i-1) * x * b - (i-1) * a) / i # c には P2(x), P3(x), ..., Pn(x) が入る a, b = b, c return c
1.2 配列 x に対して計算
今度は区間 [-1, 1] の数列 (Numpy 配列) x に対してルジャンドル関数の値を計算してみましょう。
import numpy as np x = np.linspace(-1.0, 1.0, 21) # [-1.0, -0.9, -0.8, ..., 1.0] y = np.zeros(21) # x と同じ長さの Numpy 配列を要しておく n = 3 # 3次のルジャンドル関数を計算 for i in range(0, len(x)): y[i] = lp(n, x[i])
これで y[i] には x[i] に対応する P3(x[i]) の値が入りました。
1.3 可視化
次に P3(x) のグラフを描いてみましょう。グラフを書くには matplotlib.pyplot を利用します。
%matplotlib inline import matplotlib.pyplot as plt plt.plot(x, y)
%matplotlib inline はグラフをノートブック内に表示するためのおまじないです。通常は，ノートブックの一番はじめに実行しておきます。
# ノートブックのはじめに書く定番の3行 %matplotlib inline import numpy as np import matplotlib.pyplot as plt
2 二分探索法
区間 [a, b] で定義される単調増加関数 y = f(x) があり， f(a) < 0, f(b) > 0 であることがわかっているとします。つまり，方程式 f(x) = 0 は区間 [a, b] に1つだけ解を持つとします。このとき，次のような方法で f(x) = 0 の数値解を求めることができます。
- step 1: a, b の中点を m とする。
- step 2: a, b の距離が十分近かったら中点 m を解として終了する。そうでなかったら step 3 へ。
- step 3: f(m) > 0 なら b を m で置き換え，f(m) < 0 なら a を m で置き換える。このように区間を半分に縮小し step 1 に戻る。
このような方法を「二分探索法」といいます。
例として，2次関数 f(x) = x*x-5 を区間 [0, 5] で考えてみましょう。 f(x) は x >= 0 で単調増加です。二分探索法で f(x) = 0 の解を求め，解析解 √5 と比較してみましょう。
import numpy as np def f(x): return x * x - 5.0 def find_root(func, a, b): EPS = 1.0e-7 while True: m = (a + b) * 0.5 if np.abs(a - b) < EPS: return m if func(m) >= 0.0: b = m else: a = m print(find_root(f, 0.0, 5.0), np.sqrt(5.0))
3 練習問題
3.1 方程式の解
exp(x) = x*x の解を求めよ。
3.2 方程式の解 (2)
x + exp(x) + sin(x) = 0 の解を求めよ。
3.3 チェビシェフ多項式
次数 n のチェビシェフ (Chebyshev) 多項式 Tn(x) は次の3項間漸化式で定義される。
T0(x) = 1 T1(x) = x Tn(x) = 2 * x * Tn-1(x) - Tn-2(x) (n>=2)
n と x を与えて Tn(x) の値を返す関数 tc(n, x) を定義せよ。この関数を利用して
x = -1, -0.95, -0.90, ..., 0.95, 1
における x と Tn(x) の値の組 (xi, yi) を計算し，n=5 の場合についてグラフを図示せよ。
3.4 エルミート多項式
n 次のエルミート多項式 Hn(x) を返す関数を作成せよ。任意の実数 x に対して Hn(x) は次の漸化式で定義される。
H0(x) = 1 H1(x) = 2 * x Hn(x) = 2 * x * Hn-1(x) - 2 * (n-1) * Hn-2(x) (n >= 2)
n = 1, 2, 3, に対し，関数 Hn(x) * exp(-x*x/2) のグラフを [-10, 10] の範囲でプロットせよ。
3.5 黄金比
黄金比 (1 + √(5)) / 2 は次のような連分数で表される。
1 + √5 1
------- = 1 + ---------------
2 1
1 + ----------
1
1 + -----
:
右辺全体を X とすると，右辺第2項の分母が X に等しいから
1
X = 1 + ---
X
という関係が成立する（実際，この2次方程式を解くと黄金比が得られる）。そこで次のような漸化式を考える。
1
Xn+1 = 1 + ---- (n = 1, 2, ...)
Xn
初項を X1 = 1 として Xn の値が一定値に収束するまで計算を行い，結果が黄金比に一致することを確認せよ。初項の値をいろいろに変えた場合について数値実験を行い，結果がどのように変わるか調べること。
3.6 素因数分解
与えられた正の整数 n を素因数分解するプログラムを作成せよ。例えば n = 60 のときは
2 2 3 5
のように画面に出力すればよい。
ヒント: 与えられた数を 2 からその数の平方根までの整数で順次割って行けばよい。n = 20 のときは次のようになる。
- 20 を 2 で割ると 10 余り 0
- 10 を 2 で割ると 5 余り 0
- 5 を 2 で割ると割り切れないので，次は 3 で割る
- 5 を 3 で割ると割り切れないので，次は 4 で割る
- 5 を 4 で割ると割り切れないので，次は 5 で割る（2 で割り切れなくなるまで割ったので，実 は 4 で割り切れる可能性はない）
- 5 を 5 で割ると 1 余り 0
- 1 はこれ以上割れないので終了
以上より 20 を素因数分解すると 2, 2, 5 となることがわかる。