%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
def func(x):
return -3*np.exp(-x)+np.exp(-3*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()
上では,funcの中身を修正している.plotの範囲(x1..x2)は変える必要はない. 以下では, x0が解析解であるが,これは修正している.
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 = -1/2*np.log(3.0)
list_bisec = [[0],[abs(x1-x0)]]
for i in range(0, 10):
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 +11.9306914378 -2.0000000000 -1.0000000000 -0.5000000000 +11.9306914378 -0.4644747418 -0.7500000000 -0.5000000000 +3.1367357865 -0.4644747418 -0.6250000000 -0.5000000000 +0.9160812480 -0.4644747418 -0.5625000000 -0.5000000000 +0.1407849543 -0.4644747418 -0.5625000000 -0.5312500000 +0.1407849543 -0.1809993961 -0.5625000000 -0.5468750000 +0.1407849543 -0.0251426693 -0.5546875000 -0.5468750000 +0.0565301134 -0.0251426693 -0.5507812500 -0.5468750000 +0.0153750461 -0.0251426693 -0.5507812500 -0.5488281250 +0.0153750461 -0.0049629758 -0.5498046875 -0.5488281250 +0.0051861813 -0.0049629758
import matplotlib.pyplot as plt
X = list_bisec[0]
Y = list_bisec[1]
plt.plot(X, Y)
plt.yscale("log") # y軸を対数目盛に
plt.show()
で6桁で表示させて見やすくしている.
import numpy as np
np.set_printoptions(precision=6, suppress=True)
tt=0.2
A=np.array([[1,tt,tt],[tt,1,tt],[tt,tt,1]])
b=np.array([2,2,2])
n=3
x0=np.zeros(n)
x1=np.zeros(n)
for iter in range(0, 10):
for i in range(0, n):
x1[i]=b[i]
for j in range(0, n):
x1[i]=x1[i]-A[i][j]*x0[j]
x1[i]=x1[i]+A[i][i]*x0[i]
x1[i]=x1[i]/A[i][i]
x0[i]=x1[i]
print(iter,x0)
0 [ 2. 1.6 1.28] 1 [ 1.424 1.4592 1.42336] 2 [ 1.423488 1.43063 1.429176] 3 [ 1.428039 1.428557 1.428681] 4 [ 1.428552 1.428553 1.428579] 5 [ 1.428574 1.42857 1.428571] 6 [ 1.428572 1.428571 1.428571] 7 [ 1.428571 1.428571 1.428571] 8 [ 1.428571 1.428571 1.428571] 9 [ 1.428571 1.428571 1.428571]
1.と2.を同時プロットして表示するよう工夫したバージョン.
紅茶花伝の知恵.ありがとう.
3.の答えとしては,「ピークの位置が違う」とかいう解答ではなく,
「新たに高周波側に40チャンネルあたりにピークが出現している」
とか
「その強度は1/2のはずだが,期待したほど小さくなってない」
などというコメントがあると点数が上がる.
そこに点数をつける理由は,数値計算で必要なスキルを習得するコツだから.
数値計算では,細かい数値だけでなく,プロットの図形の読み取りも大事で,
それを定性的,定量的に言葉にすることで解の性質を把握する癖がついてくるから.
%matplotlib inline
from scipy import fft
import matplotlib.pyplot as plt
import numpy as np
def func(x):
return np.cos(x/15)+np.cos(x/2)
x = np.linspace(0, 256, 256) #0から2πまでの範囲を100分割したnumpy配列
plt.plot(x, func(x), color = 'b')
plt.plot(x, np.cos(x/15), color = 'r', linewidth=0.8)
plt.plot(x, np.cos(x/2), color = 'r', linewidth=0.8)
plt.grid()
plt.show()
yy = func(x)
out = fft(yy)
def spectrum_power(x):
re, im = x.real, x.imag
return np.sqrt(re**2+im**2)
plt.plot(x,spectrum_power(out))
plt.xlim(0,128)
plt.show()
plt.plot(x,spectrum_power(out))
plt.xlim(0,128)
plt.yscale('log')
plt.show()
%matplotlib inline
def func(x):
return np.cos(x/15)+np.cos(x/2)+1/2.0*np.cos(x)
x = np.linspace(0, 256, 256) #0から2πまでの範囲を100分割したnumpy配列
plt.plot(x, func(x), color = 'b')
plt.plot(x, np.cos(x/15), color = 'r', linewidth=0.8)
plt.plot(x, np.cos(x/2), color = 'r', linewidth=0.8)
plt.grid()
plt.show()
yy = func(x)
out2 = fft(yy)
plt.plot(x,spectrum_power(out2))
plt.xlim(0,128)
plt.show()
plt.plot(x,spectrum_power(out2))
plt.xlim(0,128)
plt.yscale('log')
plt.show()
plt.plot(x,spectrum_power(out2), color = 'b')
plt.plot(x,spectrum_power(out), color = 'r')
plt.xlim(0,128)
# plt.yscale('log')
plt.show()
章末課題に挙げていたlogを変更した問題です. 解答例としてはこっちの方が見やすい(打ち間違いが少ない)でしょうね. 桁数制御してないので,ちょっと誤差が出ています.
import numpy as np
np.set_printoptions(precision=7, suppress=True)
print(np.tan(np.pi/4))
1.0
x0 = 0.6
x1 = 0.7
x2 = 0.8
x3 = 0.9
y0 =np.tan(x0)
y1 =np.tan(x1)
y2 =np.tan(x2)
y3 =np.tan(x3)
print(y0,y1,y2,y3)
0.684136808342 0.842288380463 1.02963855705 1.26015821755
f1_01 = (y1-y0)/(x1-x0) #1.581516
f1_12 = (y2-y1)/(x2-x1) #1.873506
f2_012 = (f1_12-f1_01)/(x_2-x_0)
print(f2_012)
f1_23 = (y3-y2)/(x3-x2) # 2.305190
f2_123 = (f1_23-f1_12)/(x_3-x_1)
print(f2_123)
1.45993022329 2.15847419563
x = 0.7853982
F_0 = y0
F_0 + (x-x_0)*(f1_01)+(x-x_0)*(x-x_1)*f2_012
1.0004615735170512
f3_0123 = (f2_123-f2_012)/(x3-x0)#2.328233
F_0 + (x-x_0)*(f1_01)+(x-x_0)*(x-x_1)*f2_012+(x-x_0)*(x-x_1)*(x-x_2)*f3_0123
0.99992326268286946