#!/usr/bin/env python
# coding: utf-8
#
#
#
# 数値計算2016解答例
#
#
#
# file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/fsolve
#
# https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python
#
# 2019/12/06 西谷滋人
#
# Table of Contents
#
# # 数値解の収束性
# In[12]:
get_ipython().run_line_magic('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が解析解であるが,これは修正している.
# In[13]:
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()
# In[14]:
import matplotlib.pyplot as plt
X = list_bisec[0]
Y = list_bisec[1]
plt.plot(X, Y)
plt.yscale("log") # y軸を対数目盛に
plt.show()
# # Gauss-Seidelの収束性
#
# * 以下のcodeはGauss-Seidelに直したコード.
# * np.set_printoptions(precision=6, suppress=True)
# で6桁で表示させて見やすくしている.
# * tt=0.2とした場合は,7回程度のループで収束しているのがわかる.
# * ttを変更した場合は,これで収束せず,iterのrangeを100回程度まで取らないと収束する回数は見えない.
#
# In[22]:
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)
# # FFTの強度表示
# ```
# 1.と2.を同時プロットして表示するよう工夫したバージョン.
# 紅茶花伝の知恵.ありがとう.
# 3.の答えとしては,「ピークの位置が違う」とかいう解答ではなく,
# 「新たに高周波側に40チャンネルあたりにピークが出現している」
# とか
# 「その強度は1/2のはずだが,期待したほど小さくなってない」
# などというコメントがあると点数が上がる.
#
# そこに点数をつける理由は,数値計算で必要なスキルを習得するコツだから.
# 数値計算では,細かい数値だけでなく,プロットの図形の読み取りも大事で,
# それを定性的,定量的に言葉にすることで解の性質を把握する癖がついてくるから.
# ```
# In[40]:
get_ipython().run_line_magic('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()
# In[41]:
get_ipython().run_line_magic('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()
# In[43]:
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を変更した問題です.
# 解答例としてはこっちの方が見やすい(打ち間違いが少ない)でしょうね.
# 桁数制御してないので,ちょっと誤差が出ています.
# In[97]:
import numpy as np
np.set_printoptions(precision=7, suppress=True)
print(np.tan(np.pi/4))
# In[98]:
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)
# In[99]:
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)
# In[100]:
x = 0.7853982
F_0 = y0
F_0 + (x-x_0)*(f1_01)+(x-x_0)*(x-x_1)*f2_012
# In[101]:
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
# In[ ]: