#!/usr/bin/env python
# coding: utf-8
#
#
#
# 数値計算試験予行演習問題
#
#
#
# 2020/12/11 実施
#
# cc by Shigeto R. Nishitani 2020
#
#
#
#
# # 簡単な行列計算:25点
#
# 次の行列
# $
# A = \left(\begin{array}{ccc}
# 5 & -2 & -2 \\
# 2 & 1 & -2 \\
# 6 & -2 & -3
# \end{array}
# \right)
# $
# の固有値と固有ベクトルを求めよ.
# また,対角化行列$P$ を用いて,ドット演算 により$P^{-1}.A.P$が対角化されることを確かめよ.
#
# In[1]:
import numpy as np
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})
A = np.array([[5, -2, -2], [2,1,-2], [6,-2,-3]])
l, P = np.linalg.eig( A )
print(l)
print(P)
# In[2]:
from pprint import pprint
dA = np.dot(np.linalg.inv(P),np.dot(A,P))
np.set_printoptions(precision=3, suppress=True)
print(format(dA))
# $P^{-1}.A.P$で固有値(-1,3,1)に対角化されていることが確認できる.
# # FFTの強度表示:25点
#
# FFTによって非整合波の重ね合わせを周波数分解したときの様子を観察する.
# 1. $\displaystyle \cos\left(\frac{x}{17}\right)$と
# $\displaystyle \cos\left(\frac{x}{3}\right)$を重ね合わせた関数に
# FFTをかけて(スペクトル)強度を周波数で表示せよ.
# 1. 上記2関数にさらに$\displaystyle \frac{1}{4}\cos(2x)$を
# 重ね合わせた関数にFFTをかけて(スペクトル)強度を周波数で表示せよ.
# 1. 上で求めた1. 2.のスペクトル強度の違いを述べよ.
#
# In[3]:
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.sin(x/13)+np.sin(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/17), color = 'r', linewidth=0.8)
plt.plot(x, np.cos(x/3), 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[4]:
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.sin(x/13)+np.sin(x/2)+1/4*np.cos(2*x)
x = np.linspace(0, 256, 256) #0から2πまでの範囲を100分割したnumpy配列
plt.plot(x, func(x), color = 'b')
plt.plot(x, np.cos(x/17), color = 'r', linewidth=0.8)
plt.plot(x, np.cos(x/3), color = 'r', linewidth=0.8)
plt.plot(x, 1/4*np.cos(2*x), 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()
# 2のスペクトルは1に高周波のノイズが乗った感じ.logでプロットすると周波数90あたりに変な落ち込みがあるね.なんだろう.誤差かな..?
# # Gauss-Seidelの収束性:25点
# 初期値を$[0,0,0]^{t}$として,
# $A(tt)x=b$に
# ガウス・ザイデルによる連立一次方程式の反復解法プログラムを適用する.
# ただし,
# \begin{equation}
# A(tt)=
# \left(
# \begin{array}{ccc}
# 5&tt&tt \\
# tt&5&tt \\
# tt&tt&5
# \end{array}
# \right)
# ,
# b=
# \left(
# \begin{array}{c}
# 2 \\
# 2 \\
# 2 \\
# \end{array}
# \right)
# \end{equation}
# である.$tt=1,2.5,4.5$に対して有効数字6桁の解を得るための反復回数を求めよ.
# (E.クライツィグ著「数値解析」(培風館,2003), p.89, 問題2.3-9改)
# In[5]:
import numpy as np
np.set_printoptions(precision=10, suppress=True)
for tt in [1.0,2.5,4.5]:
print("tt= ", tt)
A=np.array([[5,tt,tt],[tt,5,tt],[tt,tt,5]])
b=np.array([2,2,2])
n=3
x0=np.zeros(n)
for iter in range(0, 10):
for i in range(0, n):
x1 = b[i]
for j in range(0, n):
x1 -= A[i][j]*x0[j]
x1 += A[i][i]*x0[i]
x1 /= A[i][i]
x0[i] = x1
print(iter,x0)
# tt=1では6回程度で,tt=2.5では10回程度で収束しているが,tt=4.5では10回でもまだ収束していない.
# 100回まで繰り返すと収束しており,大体80回程度で収束している.
# お互いの差をとってやれば収束判定が可能であろう.
#
# In[6]:
for tt in [4.5]:
print("tt= ", tt)
A=np.array([[5,tt,tt],[tt,5,tt],[tt,tt,5]])
b=np.array([2,2,2])
n=3
x0=np.zeros(n)
for iter in range(0, 100):
for i in range(0, n):
x1 = b[i]
for j in range(0, n):
x1 -= A[i][j]*x0[j]
x1 += A[i][i]*x0[i]
x1 /= A[i][i]
x0[i] = x1
print(iter,x0)
if abs(x0[0]-x0[1])<10**(-6):
break
# 第1項と第2項の差で収束判定をすると,78回目で収束している.
#
# これらの結果から,Gauss-Seidel法によって解を求める時には,
# 対角成分と非対角成分が近い値では,収束が遅くなることが確認できた.
# # 常微分方程式:25点
#
# スーパーカーの加速性能をEuler法で求める.
# 1. 空気抵抗のない場合に,自動車を静止状態から加速度$2 \times 17.35$ [m/sec$^2$]で
# 加速した場合の移動距離と速度の変化をtt=0..5でプロットせよ.
# 1. 空気抵抗としてcc=0.5とした場合の移動距離と速度の変化をtt=0..10でプロットせよ.
# 1. 上問1. 2.のプロットの違いを,「ゼロよん」と呼ばれる400mを通過するまでの時間や,
# その時の時速を用いて,説明せよ.
#
# 空気抵抗は,テキストにある通り,速度に比例するとして,その係数は規格化された係数とする.
# 時間の刻み幅dt=0.1[sec]程度をとれ.
# In[7]:
import numpy as np
def euler(x0, v0):
v1 = v0 + g * dt
x1 = x0 + v0 * dt
return x1, v1
import matplotlib.pyplot as plt
def my_plot(xx, vv, tt):
plt.plot(tt, xx, color = 'b', linestyle='--',label="distance")
plt.plot(tt, vv, color = 'r', label="velocity")
plt.legend()
plt.xlabel('time')
plt.ylabel('distance and velocity')
plt.grid()
plt.show()
g, dt =2*17.35, 0.1
tt,xx,vv=[0.0],[0.0],[0.0]
t = 0.0
for i in range(0,50):
t += dt
x, v = euler(xx[-1],vv[-1])
tt.append(t)
xx.append(x)
vv.append(v)
# print(xx)
# print(vv)
my_plot(xx, vv, tt)
# In[8]:
print(xx[-2])
print((50-2)*0.1)
print(vv[-2]*60*60/1000)
# In[2]:
import numpy as np
import matplotlib.pyplot as plt
def my_plot(xx, vv, tt):
plt.plot(tt, xx, color = 'b', linestyle='--',label="distance")
plt.plot(tt, vv, color = 'r', label="velocity")
plt.legend()
plt.xlabel('time')
plt.ylabel('distance and velocity')
plt.grid()
plt.show()
def euler2(x0, v0):
v1 = v0 + (-cc * v0+ g) * dt
x1 = x0 + v0 * dt
return [x1, v1]
g, dt, cc =2*17.35, 0.1, 0.5
tt,xx,vv=[0.0],[0.0],[0.0]
t = 0.0
for i in range(0,100):
t += dt
x, v = euler2(xx[-1],vv[-1])
tt.append(t)
xx.append(x)
vv.append(v)
my_plot(xx, vv, tt)
# In[5]:
print(xx[-24])
print((100-24)*0.1)
print(vv[-24]*60*60/1000)
# xx,vvの値をプリントすることで400mを通過する時間が読み取れる.
#
# 空気抵抗がない場合は(50-2)/10秒で4.8秒,
# 空気抵抗がある場合は(100-24)/10秒で7.6秒である.
#
# またそれぞれのその時点での時速は,
# 空気抵抗がない場合は,612km/h, 空気抵抗がある場合は,245km/hである.現実問題として空気抵抗がない状態で,600km/hでは新幹線の最速を抜くが,ジェット機はまだその先にある.自動車にジェットエンジンを積むのは現実的ではないが,空気抵抗は抑えることができそうである.空気抵抗をいかに抑えることができるかが,車の性能に大きく利いていることが実感できる.F1のマシンの空気抵抗はどの程度なのだろう.
# In[ ]: