https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python
cc by Shigeto R. Nishitani 2017
概要:補間と近似¶
単純な2次元データについて補間と近似を考える．補間はたんに点を「滑らかに」つなぐことを，近似はある関数にできるだけ近くなるように「フィット」することを言う．補間はIllustratorなどのドロー系ツールで曲線を引くときの，ベジエやスプライン補間の基本となる．本章では補間とそれに密接に関連した積分について述べる．
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
fig, (axL, axR) = plt.subplots(ncols=2, figsize=(10,4))
x = np.array([0,1,2,3])
y = np.array([1.2,3.2,3.8,3.2])
for i in range(0,4):
axL.plot(x[i],y[i],'o',color='r')
axL.hlines(0, -1, 4, linewidth=1)
axL.vlines(0, -1, 4, linewidth=1)
axL.set_title('Interpolation')
axL.grid(True)
x = np.array([0,1,2,3])
y = np.array([0.2,1.2,1.8,3.2])
for i in range(0,4):
axR.plot(x[i],y[i],'o',color='r')
axR.hlines(0, -1, 4, linewidth=1)
axR.vlines(0, -1, 4, linewidth=1)
axR.set_title('Fitting')
axR.grid(True)
# fig.show()
多項式補間(polynomial interpolation)¶
データを単純に多項式で補間する方法を先ず示そう．$N+1$点を$N$次の多項式でつなぐ．この場合の補間関数は，
$$ F \left(x \right)={\sum_{i=0}^{N} } a _{i }x ^{i }=a_{0}+a_{1}x +a_{2}x^{2}+\cdots +a_{N}x^{N} $$である．データの点を$(x_{i},\,y_{i}),i=0..N$とすると
$$ \begin{array}{cl} a _{0}+a _{1}x _{0}+a _{2}x _{0}^{2}+\cdots +a _{N }x _{0}^{N }& =y _{0}\\ a _{0}+a _{1}x _{1}+a _{2}x _{1}^{2}+\cdots +a _{N }x _{1}^{N }& =y _{1}\\ \vdots& \\ a _{0}+a _{1}x _{N}+a _{2}x _{N}^{2}+\cdots +a _{N }x _{N}^{N }& =y _{N} \end{array} $$が，係数 $a_i$を未知数と見なした線形の連立方程式となっている． 連立方程式の係数行列$\mathbf{A}$(ヴァンデルモンド行列と呼ばれる)は $$ A=\left[ \begin{array}{ccccc} 1&x_0&x_0^2&\cdots&x_0^N \\ 1&x_1&x_1^2&\cdots&x_1^N \\ \vdots& & & & \vdots \\ 1&x_N&x_N^2&\cdots&x_N^N \end{array} \right] $$ となる． 係数$a_i$および データ$y_i$をそれぞれベクトルとみなすと
$A$と$y$から$a_i$を導出：
により未知数ベクトル$a_i$が求まる．これは単純に，前に紹介した Gauss の消去法や LU 分解で解ける．
Lagrange(ラグランジュ) の内挿公式¶
多項式補間は手続きが簡単であるため，計算間違いが少なく，プログラムとして組むのに適している．しかし，あまり"みとうし"のよい方法とはいえない．その点，Lagrange(ラグランジュ)の内挿公式は見通しがよい．これは
$$ F(x)= \sum_{k=0}^N \frac{ \prod_{j \ne k} (x-x_j)} { \prod_{j \ne k} (x_k-x_j)} y_k =\sum_{k=0}^N \frac{ \frac{ (x-x_0)(x-x_1)\cdots(x-x_N)}{ (x-x_k)}} {\frac{ (x_k-x_0)(x_k-x_1)\cdots(x_k-x_N)}{ (x_k-x_k)}} y_k $$と表わされる．数学的に 2つ目の表記は間違っているが，先に割り算を実行すると読み取って欲しい．これは一見複雑に見えるが，単純な発想から出発している．求めたい関数$F(x)$を
$$ F(x) = y_0 L_0(x)+y_1 L_1(x)+y_2 L_2(x) $$とすると
$$ \begin{array}{ccc} L_0(x_0) = 1 &L_0(x_1) = 0 &L_0(x_2) = 0 \\ L_1(x_0) = 0 &L_1(x_1) = 1 &L_1(x_2) = 0 \\ L_2(x_0) = 0 &L_2(x_1) = 0 &L_2(x_2) = 1 \end{array} $$となるように関数$L_i(x)$を決めればよい．これを以下のようにとればLagrangeの内挿公式となる．
$L_0(x)$：
$L_1(x)$：
$L_2(x)$：
$F(x)$：
python code¶
pythonではLagrange補間はinterpolate.lagrangeで用意されている．
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
# もとの点
x = np.array([0,1,2,3])
y = np.array([1,2,3,-2])
for i in range(0,4):
plt.plot(x[i],y[i],'o',color='r')
# Lagrange補間
f = interpolate.lagrange(x,y)
print(f)
x = np.linspace(0,3, 100)
y = f(x)
plt.plot(x, y, color = 'b')
plt.grid()
plt.show()
3 2 -1 x + 3 x - 1 x + 1
/Users/bob/opt/anaconda3/lib/python3.8/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.4
warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
Newton(ニュートン) の差分商公式¶
もう一つ有名なNewton(ニュートン) の内挿公式は，
$$ \begin{array}{rc} F (x )&=F (x _{0})+ (x -x _{0})f _{1}\lfloor x_0,x_1\rfloor+ (x -x _{0})(x -x _{1}) f _{2}\lfloor x_0,x_1,x_2\rfloor + \\ & \cdots + \prod_{i=0}^{n-1} (x-x_i) \, f_n \lfloor x_0,x_1,\cdots,x_n \rfloor \end{array} $$となる．ここで$f_i \lfloor\, \rfloor$ は次のような関数を意味していて，
$$ \begin{array}{rl} f _{1}\lfloor x_0,x_1\rfloor &= \frac{y_1-y_0}{x_1-x_0} \\ f _{2}\lfloor x_0,x_1,x_2\rfloor &= \frac{f _{1}\lfloor x_1,x_2\rfloor- f _{1}\lfloor x_0,x_1\rfloor}{x_2-x_0} \\ \vdots & \\ f _{n}\lfloor x_0,x_1,\cdots,x_n\rfloor &= \frac{f_{n-1}\lfloor x_1,x_2\cdots,x_{n}\rfloor- f _{n-1}\lfloor x_0,x_1,\cdots,x_{n-1}\rfloor}{x_n-x_0} \\ \end{array} $$差分商と呼ばれる．得られた多項式は，Lagrange の内挿公式で得られたものと当然一致する．Newtonの内挿公式の利点は，新たなデータ点が増えたときに，新たな項を加えるだけで，内挿式が得られる点である．
# https://stackoverflow.com/questions/14823891/newton-s-interpolating-polynomial-python
# by Khalil Al Hooti (stackoverflow)
def _poly_newton_coefficient(x,y):
"""
x: list or np array contanining x data points
y: list or np array contanining y data points
"""
m = len(x)
x = np.copy(x)
a = np.copy(y)
for k in range(1,m):
a[k:m] = (a[k:m] - a[k-1])/(x[k:m] - x[k-1])
return a
def newton_polynomial(x_data, y_data, x):
"""
x_data: data points at x
y_data: data points at y
x: evaluation point(s)
"""
a = _poly_newton_coefficient(x_data, y_data)
n = len(x_data) - 1 # Degree of polynomial
p = a[n]
for k in range(1,n+1):
p = a[n-k] + (x -x_data[n-k])*p
return p
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
# もとの点
x = np.array([0,1,2,3])
y = np.array([1,2,3,-2])
for i in range(0,4):
plt.plot(x[i],y[i],'o',color='r')
print(_poly_newton_coefficient(x,y))
xx = np.linspace(0,3, 100)
yy = newton_polynomial(x, y, xx)
plt.plot(xx, yy, color = 'b')
plt.grid()
plt.show()
[ 1 1 0 -1]
Newton補間と多項式補間の一致の検証¶
関数$F(x)$を$x$の多項式として展開．その時の，係数の取るべき値と，差分商で得られる値が一致．
maple
> restart: F:=x->f0+(x-x0)*f1p+(x-x0)*(x-x1)*f2p;
maple
> F(x1);
sf1p:=solve(F(x1)=f1,f1p);
f20の取るべき値の導出
maple
> sf2p:=solve(F(x2)=f2,f2p);
fac_f2p:=factor(subs(f1p=sf1p,sf2p));
ニュートンの差分商公式を変形
maple
> ff11:=(f0-f1)/(x0-x1);
ff12:=(f1-f2)/(x1-x2);
ff2:=(ff11-ff12)/(x0-x2);
fac_newton:=factor(ff2);
二式が等しいかどうかをevalbで判定
maple
> evalb(fac_f2p=fac_newton);
数値積分 (Numerical integration)¶
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
x = np.array([1,2,3,4])
y = np.array([2.8,3.4,3.8,3.2])
for i in range(0,4):
plt.plot(x[i],y[i],'o',color='r')
plt.hlines(0, -1, 5, linewidth=1)
plt.vlines(0, -1, 5, linewidth=1)
plt.grid(True)
# 数値積分の模式図
積分，
$$ I = \int_a^b f(x) dx $$を求めよう．1次元の数値積分法では連続した領域を細かい短冊に分けて，それぞれの面積を寄せ集めることに相当する．分点の数を N とすると，
$$ \begin{array}{c} x_i = a+ \frac{b-a}{N} i = a + h \times i \\ h = \frac{b-a}{N} \end{array} $$ととれる．そうすると，もっとも単純には，
$$ I_N = \left\{\sum_{i=0}^{N-1} f(x_i)\right\}h = \left\{\sum_{i=0}^{N-1} f(a+i \times h)\right\}h $$となる．
中点則 (midpoint rule)¶
中点法 (midpoint rule) は，短冊を左端から書くのではなく，真ん中から書くことに対応し，
$$ I_N = \left\{\sum_{i=0}^{N-1}f\left(a+\left(i+\frac{1}{2}\right) \times h\right)\right\}h $$となる．
台形則 (trapezoidal rule)¶
さらに短冊の上側を斜めにして，短冊を台形にすれば精度が上がりそうに思う． その場合は，短冊一枚の面積$S_i$は，
$$ S_i = \frac{f(x_i)+f(x_{i+1})}{2}h $$で求まる．これを端から端まで加えあわせると，
$$ i_N =\sum _{i=0}^{N-1}S_i =h \left\{ \frac{1}{2} f ( x_0 ) +\sum _{i=1}^{N-1}f ( x_i ) +\frac{1}{2} f \left( x_N \right) \right\} $$が得られる．
Simpson(シンプソン)則¶
Simpson(シンプソン) 則では，短冊を2次関数，
$$ f(x) = ax^2+bx+c $$で近似することに対応する．こうすると，
$$ S_i=\int _{x_i}^{x_{i+1}}f ( x )\, {dx}=\int _{x_i}^{x_{i+1}}(ax^2+bx+c)\,{dx} $$
Simpson則の導出（数式変形）：
$$ \frac{h}{6} \left\{f(x_i)+4f\left(x_i+\frac{h}{2}\right)+f(x_i+h)\right\} $$
となる．これより，
$$ I_N=\frac{h}{6} \left\{ f \left( x_0 \right) +4\,\sum _{i=0}^{N-1}f \left( x_i+\frac{h}{2} \right) +2\,\sum_{i=1}^{N-1}f \left( x_i \right) +f \left( x_N \right) \right\} $$として計算できる．ただし，関数値を計算する点の数は台形則などの倍となっている．
教科書によっては，分割数$N$を偶数にして，点を偶数番目 (even) と奇数番目 (odd) に分けて，
$$ I_{{N}}=\frac{h}{3} \left\{ f \left( x_{{0}} \right) +4\,\sum _{i={\it even}}^{N-2}f \left( x_{{i}}+\frac{h}{2} \right) +2\,\sum _{i={\it odd}}^{N-1}f \left( x_{{i}} \right) +f \left( x_{{N}} \right) \right\} $$としている記述があるが，同じ計算になるので誤解せぬよう．
数値積分のコード¶
次の積分を例に，pythonのコードを示す．
$$ \int_0^1 \frac{4}{1+x^2} \, dx $$先ずは問題が与えられたらできるだけお任せで解いてしまう．答えをあらかじめ知っておくと間違いを見つけるのが容易．プロットしてみる．
scipyで積分計算をお任せでしてくれる関数はquadで，これはFortran libraryのQUADPACKを利用している．自動で色々してくれるが，精度は1.49e_08以下．
import matplotlib.pyplot as plt
from scipy import integrate
def func(x):
return 4.0/(1.0+x**2)
x = np.linspace(0,3, 100)
y = func(x)
plt.plot(x, y, color = 'r')
plt.grid()
plt.show()
print(integrate.quad(func, 0, 1))
(3.1415926535897936, 3.4878684980086326e-14)
なんでと思うかもしれないが，
maple
>int(1/(1+x^2),x);
となるので，納得できるでしょう．
Midpoint rule(中点法)¶
def func(x):
return 4.0/(1.0+x**2)
N, x0, xn = 8, 0.0, 1.0
h = (xn-x0)/N
S = 0.0
for i in range(0, N):
xi = x0 + (i+0.5)*h
dS = h * func(xi)
S = S + dS
print(S)
3.142894729591689
Trapezoidal rule(台形公式)¶
def func(x):
return 4.0/(1.0+x**2)
N, x0, xn = 8, 0.0, 1.0
h = (xn-x0)/N
S = func(x0)/2.0
for i in range(1, N):
xi = x0 + i*h
dS = func(xi)
S = S + dS
print("{0}".format(i))
S = S + func(xn)/2.0
print(h*S)
1 2 3 4 5 6 7 3.138988494491089
Simpson's rule(シンプソンの公式)¶
def func(x):
return 4.0/(1.0+x**2)
N, x0, xn = 8, 0.0, 1.0
M = int(N/2)
h = (xn-x0)/N
Seven, Sodd = 0.0, 0.0
for i in range(1, 2*M, 2): #rangeの終わりに注意
xi = x0 + i*h
Sodd += func(xi)
print("{0}".format(i))
for i in range(2, 2*M, 2):
xi = x0 + i*h
Seven += func(xi)
print("{0}".format(i))
print(h*(func(x0)+4*Sodd+2*Seven+func(xn))/3)
1 3 5 7 2 4 6 3.141592502458707
課題¶
対数関数のニュートンの差分商補間(2014期末試験，25点)¶
2を底とする対数関数(Mapleでは$\log[2](x)$)の$F(9.2)=2.219203$をニュートンの差分商補間を用いて求める．ニュートンの内挿公式は，
$$ \begin{array}{rc} F (x )&=F (x _{0})+ (x -x _{0})f _{1}\lfloor x_0,x_1\rfloor+ (x -x _{0})(x -x _{1}) f _{2}\lfloor x_0,x_1,x_2\rfloor + \\ & \cdots + \prod_{i=0}^{n-1} (x-x_i) \, f_n \lfloor x_0,x_1,\cdots,x_n \rfloor \end{array} $$である．ここで$f_i \lfloor\, \rfloor$ は次のような関数を意味していて， $$ \begin{array}{rc} f _{1}\lfloor x_0,x_1\rfloor &=& \frac{y_1-y_0}{x_1-x_0} \\ f _{2}\lfloor x_0,x_1,x_2\rfloor &=& \frac{f _{1}\lfloor x_1,x_2\rfloor- f _{1}\lfloor x_0,x_1\rfloor}{x_2-x_0} \\ \vdots & \\ f _{n}\lfloor x_0,x_1,\cdots,x_n\rfloor &=& \frac{f_{n-1}\lfloor x_1,x_2\cdots,x_{n}\rfloor- f _{n-1}\lfloor x_0,x_1,\cdots,x_{n-1}\rfloor}{x_n-x_0} \end{array} $$ 差分商と呼ばれる．$x_k=8.0,9.0,10.0,11.0$をそれぞれ選ぶと，差分商補間のそれぞれの項は以下の通りとなる．
$$ \begin{array}{ccl|lll} \hline k & x_k & y_k=F_0( x_k) &f_1\lfloor x_k,x_{k+1}\rfloor & f_2\lfloor x_k,x_{k+1},x_{k+2}\rfloor & f_3\lfloor x_k,x_{k+1},x_{k+2},x_{k+3}\rfloor \\ \hline 0 & 8.0 & 2.079442 & & &\\ & & & 0.117783 & &\\ 1 & 9.0 & 2.197225 & & [ XXX ] &\\ & & & 0.105360 & & 0.0003955000 \\ 2 & 10.0 & 2.302585 & & -0.0050250 &\\ & & & 0.095310 & &\\ 3 & 11.0 & 2.397895 & & &\\ \hline \end{array} $$それぞれの項は，例えば，
$$ f_2\lfloor x_1,x_2,x_3 \rfloor =\frac{0.095310-0.105360}{11.0-9.0}=-0.0050250 $$で求められる．ニュートンの差分商の一次多項式の値はx=9.2で
$$ F(x)=F_0(8.0)+(x-x_0)f_1\lfloor x_1,x_0\rfloor =2.079442+0.117783(9.2-8.0)=2.220782 $$となる．
差分商補間の表中の開いている箇所[ XXX ]を埋めよ．¶
ニュートンの二次多項式¶
$$ F (x )=F (x _{0})+(x -x _{0})f _{1}\lfloor x_0,x_1\rfloor+(x -x _{0})(x -x _{1}) f _{2}\lfloor x_0,x_1,x_2\rfloor $$の値を求めよ．
ニュートンの三次多項式の値を求めよ．¶
ただし，ここでは有効数字7桁程度はとるように．(E.クライツィグ著「数値解析」(培風館,2003), p.31, 例4改)
数値積分(I)¶
次の関数
$$ f(x) = \frac{4}{1+x^2} $$を$x = 0..1$で数値積分する．
- $N$を2,4,8,…256ととり，$N$個の等間隔な区間にわけて中点法と台形則で求めよ．(15)
- 小数点以下10桁まで求めた値3.141592654との差をdXとする．dXと分割数Nとを両対数プロット(loglogplot)して比較せよ(10)
(2008年度期末試験）
解答例[7.3]¶
以下には，中点則の結果を示した．課題では，台形則を加えて，両者を比較せよ．予測とどう違うか．
import numpy as np
def func(x):
return 4.0/(1.0+x**2)
def mid(N):
x0, xn = 0.0, 1.0
h = (xn-x0)/N
S = 0.0
for i in range(0, N):
xi = x0 + (i+0.5)*h
dS = h * func(xi)
S = S + dS
return S
x, y = [], []
for i in range(1,8):
x.append(2**i)
y.append(abs(mid(2**i)-np.pi))
plt.plot(x, y, color = 'r')
plt.yscale('log')
plt.grid()
plt.show()
plt.plot(x, y, color = 'r')
plt.yscale('log')
plt.xscale('log')
plt.grid()
plt.show()
