次の関数 $$ f(x)=-4\,\exp(-x)+2\,\exp(-2\, x) $$ は$x=-1..0$に解$-\ln(2)$を持つ.二分法によって数値解を求めよ.繰り返しは10回程度でいい.また,収束の様子を片対数(logplot)でプロットせよ.
import matplotlib.pyplot as plt
import numpy as np
def func(x):
return -4*np.exp(-x)+2*np.exp(-2*x)
x1, x2 = -1.0, 0.0
x = np.linspace(x1,x2, 101)
y = func(x)
plt.plot(x, y, color = 'k')
plt.plot([x1,x2],[0,0])
plt.grid()
plt.show()
x1, x2 = -1.0, 0.0
f1, f2 = func(x1), func(x2)
print('%+15s %+15s %+15s %+15s' % ('x1','x2','f1','f2'))
print('%+15.10f %+15.10f %+15.10f %+15.10f' % (x1,x2,f1,f2))
x0 = -np.log(2.0)
list_bisec = [[0],[abs(x1-x0)]]
for i in range(0, 20):
x = (x1 + x2)/2
f = func(x)
if (f*f1>=0.0):
x1, f1 = x, f
list_bisec[0].append(i)
list_bisec[1].append(abs(x1-x0))
else:
x2, f2 = x, f
list_bisec[0].append(i)
list_bisec[1].append(abs(x2-x0))
print('%+15.10f %+15.10f %+15.10f %+15.10f' % (x1,x2,f1,f2))
list_bisec
print()
x1 x2 f1 f2 -1.0000000000 +0.0000000000 +3.9049848840 -2.0000000000 -1.0000000000 -0.5000000000 +3.9049848840 -1.1583214259 -0.7500000000 -0.5000000000 +0.4953780742 -1.1583214259 -0.7500000000 -0.6250000000 +0.4953780742 -0.4922979148 -0.7500000000 -0.6875000000 +0.4953780742 -0.0447964325 -0.7187500000 -0.6875000000 +0.2128474183 -0.0447964325 -0.7031250000 -0.6875000000 +0.0810265592 -0.0447964325 -0.6953125000 -0.6875000000 +0.0173789137 -0.0447964325 -0.6953125000 -0.6914062500 +0.0173789137 -0.0138911236 -0.6933593750 -0.6914062500 +0.0016980959 -0.0138911236 -0.6933593750 -0.6923828125 +0.0016980959 -0.0061079375 -0.6933593750 -0.6928710938 +0.0016980959 -0.0022077800 -0.6933593750 -0.6931152344 +0.0016980959 -0.0002555572 -0.6932373047 -0.6931152344 +0.0007210905 -0.0002555572 -0.6931762695 -0.6931152344 +0.0002327219 -0.0002555572 -0.6931762695 -0.6931457520 +0.0002327219 -0.0000114288 -0.6931610107 -0.6931457520 +0.0001106438 -0.0000114288 -0.6931533813 -0.6931457520 +0.0000496068 -0.0000114288 -0.6931495667 -0.6931457520 +0.0000190888 -0.0000114288 -0.6931476593 -0.6931457520 +0.0000038299 -0.0000114288 -0.6931476593 -0.6931467056 +0.0000038299 -0.0000037995
import matplotlib.pyplot as plt
X = list_bisec[0]
Y = list_bisec[1]
plt.plot(X, Y)
plt.yscale("log") # y軸を対数目盛に
plt.show()
大きな数どおしのわずかな差は,丸め誤差にとくに影響を受ける.
(E.クライツィグ著「数値解析」(培風館,2003), p.10, 問題1.1-3改)
from decimal import *
def pretty_p(result,a,b,operator):
print('context.prec:{}'.format(getcontext().prec))
print(' %20.14f' % (a))
print( '%1s%20.14f' % (operator, b))
print('-----------')
print( ' %20.14f' % (result))
getcontext().prec = 5
a=Decimal('0.81321')
b=Decimal('23.173')
c=Decimal('23.094')
pretty_p(b-c,b,c,'-')
print(b-c)
print(a/(b-c))
context.prec:5 23.17300000000000 - 23.09400000000000 ----------- 0.07900000000000 0.079 10.294
TWOPLACES = Decimal(10) ** -2
getcontext().prec = 4
a=Decimal('0.81321').quantize(Decimal(10) ** -4)
b=Decimal('23.173').quantize(Decimal('0.01'))
c=Decimal('23.094').quantize(Decimal('0.01'))
pretty_p(b-c,b,c,'-')
print(b-c)
print(a/(b-c))
context.prec:4 23.17000000000000 - 23.09000000000000 ----------- 0.08000000000000 0.08 10.16
ONEPLACES = Decimal(10) ** -1
getcontext().prec = 3
a=Decimal('0.81321').quantize(Decimal(10) ** -3)
b=Decimal('23.173').quantize(ONEPLACES)
c=Decimal('23.094').quantize(ONEPLACES)
pretty_p(b-c,b,c,'-')
print(b-c)
print(a/(b-c))
context.prec:3 23.20000000000000 - 23.10000000000000 ----------- 0.10000000000000 0.1 8.13
ZEROPLACES = Decimal(10) ** 0
getcontext().prec = 2
a=Decimal('0.81321').quantize(Decimal(10) ** -2)
b=Decimal('23.173').quantize(ZEROPLACES)
c=Decimal('23.094').quantize(ZEROPLACES)
pretty_p(b-c,b,c,'-')
print(b-c)
print(a/(b-c))
context.prec:2 23.00000000000000 - 23.00000000000000 ----------- 0.00000000000000 0
--------------------------------------------------------------------------- DivisionByZero Traceback (most recent call last) <ipython-input-17-055b09461747> in <module> 7 pretty_p(b-c,b,c,'-') 8 print(b-c) ----> 9 print(a/(b-c)) DivisionByZero: [<class 'decimal.DivisionByZero'>]
2桁の時には0割となっている.また,3桁でも相当大きなズレが出ていることが確認できる.似た同士の数値の引き算には注意が必要.
次の4点の内挿式をNewtonの差分商公式を用いて求める.
X:=[-1, 0, 1, 2];
Y:=[4., -2., -1.2, -.52];
最初の3点を用いて求めた2次の補間式は,
$$ - 2.00 - 6.00\,x+ 3.40\, \left( x+1 \right) x $$である.4点目を取り入れて3次の補間式を求めよ.
# 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([-1,0,1,2])
y = np.array([4,-2,-1.2,-0.52])
for i in range(0,4):
plt.plot(x[i],y[i],'o',color='r')
print(_poly_newton_coefficient(x,y))
xx = np.linspace(-1,2, 100)
yy = newton_polynomial(x, y, xx)
plt.plot(xx, yy, color = 'b')
plt.grid()
plt.show()
[ 4. -6. 3.4 -1.15333333]
Newtonの差分商公式の利点であるデータ点の追加に伴う煩雑さの回避は,このコードでは実感できない.pythonは,やっぱり便利なのかな.
次のようなリンクが張られたページ群のページランクを求めよ.
from pprint import pprint
from numpy import array, zeros, diagflat, dot, transpose, set_printoptions
from scipy.linalg import eig
A = array([[0,1,1,1,0],
[1,0,1,0,0],
[0,0,0,1,0],
[0,0,1,0,1],
[1,1,0,0,0]])
n = 5
diag = []
for i in range(0,n):
tmp = 0.0
for j in range(0,n):
tmp += A[i,j]
diag.append(1.0/tmp)
D = diagflat(diag)
tA = dot(transpose(A),D)
set_printoptions(formatter={'float': '{: 0.3f}'.format})
pprint(tA)
array([[ 0.000, 0.500, 0.000, 0.000, 0.500], [ 0.333, 0.000, 0.000, 0.000, 0.500], [ 0.333, 0.500, 0.000, 0.500, 0.000], [ 0.333, 0.000, 1.000, 0.000, 0.000], [ 0.000, 0.000, 0.000, 0.500, 0.000]])
初期ベクトルを等分の値にして,3度ほどホップさせた結果.
x = array([1/5,1/5,1/5,1/5,1/5])
pprint(dot(tA,dot(tA,dot(tA,x))))
array([ 0.125, 0.111, 0.269, 0.328, 0.167])
固有値を求める.固有値がソートされているか自信がないので,表示させてみた.
l, V = eig(tA)
pprint(l)
array([ 1.00000000+0.j , 0.10242347+0.50115386j, 0.10242347-0.50115386j, -0.81317748+0.j , -0.39166946+0.j ])
l[0]に対応する最大固有値のベクトルを取り出す. さらに,初期ベクトルからのホップと比較するために値を揃えている. だいたい一致しているが,ホップ数が少ないので一致はそれほど高くない.
v0 = V[:,0]
pprint(v0[0]/0.294*v0)
array([ 0.29497881-0.j, 0.26220338-0.j, 0.55718219-0.j, 0.65550846-0.j, 0.32775423-0.j])
これより,ページランクは, [4, 3, 5, 1, 2] の順になる.
?eig