« 2020年5月 | トップページ | 2020年7月 »

2020年6月

2020年6月30日 (火)

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その4) ローレンツ96モデル 36変数の方程式 

ローレンツ方程式(ローレンツモデル)と言えば3変数のものが有名ですが多変数のローレンツ96というモデルもある。

https://en.wikipedia.org/wiki/Lorenz_96_model

F=8でカオスになるとか。

dx_i/dt=(x_i+1 - x_i-2)*x_i-1 -x_i +F 

では計算。Wikipediaではodeintだがdop853で。やっぱりちょっと波形が違うな、、、

Lorenz96

 

プログラムはこちら:

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


N=36
F=8


def lorenz96(t, x): #odeintのときとt,xの並びが逆
"""Lorenz 96 model."""
# Compute state derivatives
    d = np.zeros(N)
# First the 3 edge cases: i=1,2,N
    d[0] = (x[1] - x[N-2]) * x[N-1] - x[0]
    d[1] = (x[2] - x[N-1]) * x[0] - x[1]
    d[N-1] = (x[0] - x[N-3]) * x[N-2] - x[N-1]
# Then the general case
    for i in range(2, N-1):
        d[i] = (x[i+1] - x[i-2]) * x[i-1] - x[i]
# Add the forcing term
    d = d + F

# Return the state derivatives
return d

t0=0.
tmax=30.
dt=0.01
x0 = F * np.ones(N) # Initial state (equilibrium)
x0[19] += 0.01 # Add small perturbation to 20th variable
t = np.arange(t0, tmax, dt)

solver=ode(lorenz96)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆

sol= np.zeros((len(t),N))
sol[0] = x0

k=1
while solver.successful() and solver.t < tmax-dt:
    solver.integrate(t[k])
    sol[k] = solver.y
    k+= 1


# Plot
fig = plt.figure(figsize=(12,12))
ax = fig.gca(projection='3d')

ax.plot(sol[:,0], sol[:,1], sol[:,2])

ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz96 Attractor(DOP853)")

plt.show()

 

2020年6月29日 (月)

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その3) Aizawa Attractor 

今回はみかん袋(あみあみ)を丸めたようなAizawa attractor。

http://analog-ontology.blogspot.com/2015/06/the-aizawa-attractor-iv.html

dx/dt = (z-b)*x - d*y
dy/dt = d*x+(z-b)*y
dz/dt = c+a*z-(z^3)/3 - (x^2 +y^2)*(1+e*z)+f*z*x^3

のような形で、計算結果はこちら。

Aizawa

プログラムはこんな感じで。

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def aizawa(t, x, a, b, c, d, e, f): #odeintのときとt,xの並びが逆
    x_dot = (x[2]-b)*x[0] - d*x[1]
    y_dot = d*x[0]+(x[2]-b)*x[1]
    z_dot = c+a*x[2]-(x[2]**3)/3 - (x[0]**2 +x[1]**2)*(1+e*x[2])+f*x[2]*x[0]**3
    return [x_dot, y_dot, z_dot]

t0=0
tmax = 1000
N=100000

x0=[0.1, 0.0, 0.0]
solver=ode(aizawa)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(0.95,0.7,0.6,3.5,0.25,0.1)

t=np.linspace(0, tmax, N)
sol= np.empty((N, 3))
sol[0] = x0

k=1
while solver.successful() and solver.t < tmax:
    solver.integrate(t[k]) 
    sol[k] = solver.y
    k+= 1


# Plot
fig = plt.figure(figsize=(12,12))
ax = fig.gca(projection='3d')

ax.plot(sol[:,0], sol[:,1], sol[:,2], lw=0.5)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Aizawa Attractor(DOP853)")

plt.show()

2020年6月28日 (日)

Qualcommのリファレンスデザインは富士通が作ってるが、全く売れる見込みもないということでスマホが出始めのころの国内メーカの変化に対応できなくて没落していった事例をちょっとだけ思い出す、、、

まあ

How Qualcomm Snapdragon Modular Platforms can drive 5G rollout and adoption

というのはあるんですけど、富士通はこれでスマートフォンの世界で富士通が返り咲くとでも思ったのか、、、という。私の知る限りまともなところへの採用0です。

Ca_takada_halfinline

いやまあ終わってるんですけど。業界に長くいるので、ガラケー(あれも今となっては技術的には終わっているが、一つだけ、パカパカするところのケーブル、フレキ類だけは今でも結構生きのこっている)からスマートフォンに変わるときにそれまで偉そうにしていたメーカーの技術力のなさが露呈したという、、、

私がちょっと(かなり?)携わった事例です。メーカー名は言わないですが、So,Sh,P,N,Fのどれかです。

・Qcomのリファレンスデザインで電源管理IC(PMIIC)はキャパシタ(積層コンデンサ)を使いすぎだ!俺の設計では9割減らせる!

 →やってみたら動かない、、、急遽ヘルプが来てランドがないのにいっぱいキャパシタを乗せた

  (これは本気でひどかった。日本の携帯終わったな、と思った最初の事例)

・オーディオIC、薄くするためにダイを削ったら音が変わるので許さない(デジタルのCODECですよ?)!

  オカルト部長と私は言っていた。百歩譲ってアナログのアンプなどなら0.0000000000001%の可能性はあるかもしれないが、

  CODEC ICですよ、、、

・Qcomの回路図を理解していないで、勝手に1個キャパシタを抜いて(どういう機能なのか全く理解せずに)不具合が起きたら

 その先にあるモジュールのせいにする。

 これは私が直接かかわった例。夜中の3時に会議したという、、、

 ”御社が勝手にキャパシタ抜いたから不具合が起きてるのですよ”と言ったら怒られた。お前らのせいだ!何とかしろと。

 その部長さんは首になって弊社に再就職を希望したがもちろん蹴られた。工場でもどこでもいいから雇ってくれとか言っていた。

・キャパシタ関連は某最大手でi何とかというスマホを作っているメーカーも同じように、”不具合が起きた。バイパスコンデンサに突入電流が流れているのが原因だ!モジュールの不具合だ!”と言ってきた。

 あ、それが正常動作です。と言って中学生向けの電気回路の本から抜粋して(英訳して)説明した。

 ※こんなやつが私の数倍給料をもらっているという、、、インド人も優秀な人ばかりではない。

 

なんかキャパシタだらけになったけど、本当は高周波(RF)関係でもう国内は終わっている事例がたくさん(これはぼかしても分かってしまうのでかけない、、、)

これだけ。

・キャリアアグリゲーションの回路構成を、”あんなの使われるはずない”とかいっていた国内メーカの方々ですが

 世界標準になってまったく取り残された。

 ※この辺で国内メーカーは完全に終わったと思う。3GPPでもどこでも、国内最大手キャリアが何も影響力がないことに

 日本のメーカーは気づかなかった。しょせんPDCくらいの時代の人たちですよ。。。

 (いまだにPDCの技術力がすごいというじいさんが最近いたので、、、CDMAはスペックが緩いとか、、、5Gの時代にPre2Gの話をされても、、、)

 

 

・放熱が必要なことを全く理解していないでカイロになる。

https://bbs.kakaku.com/bbs/J0000017526/SortID=20107343

/

・アンテナの設計がへぼすぎて全く規格を満たしていないが、認証機関との関係性からごまかす(これは昔のK国のSMSNがほぼそうだった)

 指摘すると切れられる。※アンテナの特性が出てない事例は数限りなくあります。

・一番多い事例は数10k台しかでないくせに、iPhoneでやっているようなカスタム設計を要求してそして偉そう。

 じいさんに多い。ガラケー時代はそうやっていたんだろうな、とか。数億台でるからカスタム設計するのであって

 一体何を偉そうにしているのか。

・最近の事例:RFIを出して、中華圏メーカが適当に答えた回答をうのみにして採用。

 性能は全く満たしていない。中華圏メーカも数が出ないから適当に答えていて、真面目に見積もったメーカーは馬鹿を見る。

 なので出来上がったスマホの性能は終わっている。

・偉そうにしているスマホメーカーが大手の部品メーカーと会議して、企画台数を聞かれた瞬間にその部品メーカーが一言もしゃべらずに帰る。(だから数億台の企画台数もないのに偉そうにするのはやめてくれ、、、これは弊社じゃないです。そんな失礼なことはしないです。)

 

あとはiPhoneなんか日本で使われるはずがないといった偉い人ですかね。

 

2020年6月27日 (土)

Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その2) やはり最初はローレンツ方程式。

昨日で大体使い方が分かったので、まずは一番メジャーなローレンツ方程式か。

dop853で計算してみる。本当はrtolとatolも設定しないといけないが手抜きで。

Lorenz_dop853_20200623205801

おなじみの形が得られた。プログラムはこんな感じで。

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def lorenz(t, x, s, r, b): #odeintのときとt,xの並びが逆
    x_dot = s*(x[1] - x[0])
    y_dot = r*x[0] - x[1] - x[0]*x[2]
    z_dot = x[0]*x[1] - b*x[2]
    return [x_dot, y_dot, z_dot]

t0=0
tmax = 100
N=100000

x0=[0.1, 0.1, 0.1]
solver=ode(lorenz)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(10,28,8/3)

t=np.linspace(0, tmax, N)
sol= np.empty((N, 3))
sol[0] = x0

k=1
while solver.successful() and solver.t < tmax:
    solver.integrate(t[k])
    sol[k] = solver.y
    k+= 1


# Plot
fig = plt.figure(figsize=(12,12))
ax = fig.gca(projection='3d')

ax.plot(sol[:,0], sol[:,1], sol[:,2])
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor(DOP853)")

plt.show()

---

odeintでデフォルトで計算してみた結果も。やっぱり結構違うね。

Lorenz_odeint

 

 

2020年6月26日 (金)

Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その1) まずは準備編

PythonでWikipediaに乗っているカオスの写像リストを描いて遊んでいたが、そろそろネタ切れ、、、

そこで写像から常微分方程式に移る。とはいってもodeintを使うと面白くないので、このサイトでいつもやっている(ExcelのVBAに移植した)

DOP853(Dormand&Prince, ルンゲクッタ8次)を使ってみよう。

(参考リンク)

 Dormand-Prince(ルンゲ・クッタ8次)の有名なルーチンDOP853をExcel VBAに移植

Dormand-Prince(ルンゲ・クッタ8次)DOP853をExcel VBAに移植(説明編)

では準備として単振動する

dx/dt=y

dy/dt=-x

を計算してみよう。これを100周期(2π*100)繰り返して厳密解のcos(t),sin(t)と比較してみる。

Pythonのプログラムはこんな感じ。odeintを使っている方は、そもそも関数のtとxの並びが反転しているとか、1ステップずつ計算しないといけないとか、初期値の与え方がまた並びが反転しているとか謎が多い、、、

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def f1(t, x): #odeintのときとt,xの並びが逆
    x_dot = x[1]
    y_dot = -x[0]
    return [x_dot, y_dot]

t0=0
tmax = 2.*np.pi*100
N=10000

x0=[0., 1.0]
solver=ode(f1)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆


t=np.linspace(0, tmax, N)
sol= np.empty((N, 2))
sol[0] = x0

k=1
while solver.successful() and solver.t < tmax:
    solver.integrate(t[k])
    sol[k] = solver.y
    k+= 1


# Plot
fig = plt.figure(figsize=(12,12))

plt.xlabel("Time")
plt.ylabel("Numerical - Exact")
plt.title("DOP853")
plt.grid()
plt.plot(t,sol[:,0]-np.sin(t))
plt.plot(t,sol[:,1]-np.cos(t))


plt.show()

 

結果はこちら。100周期繰り返しても誤差が1e-14オーダー。

Sin_dop853

ではodeintを使うと?

 

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def f2(x, t):
    x_dot = x[1]
    y_dot = -x[0]
    return [x_dot, y_dot]

t0=0
tmax = 2.*np.pi*100
N=10000

t=np.linspace(0, tmax, N)

x0=[0., 1.0]
x=odeint(f2, x0, t)
# Plot
fig = plt.figure(figsize=(12,12))
plt.xlabel("Time")
plt.ylabel("Numerical - Exact")
plt.title("odeint")
plt.grid()
plt.plot(t,sol[:,0]-np.sin(t))
plt.plot(t,sol[:,1]-np.cos(t))

plt.plot(t,x[:,0]-np.sin(t))
plt.plot(t,x[:,1]-np.cos(t))
plt.show()

 

簡単なのだが誤差は0.00001くらい。

Sin_odeint

やっぱりこの手のStiffじゃない問題ならDOP853がすごく誤差が小さい。

ではこれを使っていろいろ計算してみる(続く)。

2020年6月25日 (木)

Pythonの3Dのscatter+cmap+SciPyの特殊関数で水素原子の波動関数の電子分布を描いてみる。

Pythonのmatplotlib.pyplotの散布図を描くscatter, 色として実数を入れてcmapを設定すると値に応じて色の濃淡がつけられることを知った。

それなら実験的に水素原子の波動関数を描いてみよう。

https://ja.wikipedia.org/wiki/%E6%B0%B4%E7%B4%A0%E5%8E%9F%E5%AD%90%E3%81%AB%E3%81%8A%E3%81%91%E3%82%8B%E3%82%B7%E3%83%A5%E3%83%AC%E3%83%BC%E3%83%87%E3%82%A3%E3%83%B3%E3%82%AC%E3%83%BC%E6%96%B9%E7%A8%8B%E5%BC%8F%E3%81%AE%E8%A7%A3

波動関数は、、、ああめんどくさいなと思ったら関数を作ってくれている方がいたのでこのCalcPsiを使わしてもらった。

https://amorphous.tf.chiba-u.jp/lecture.files/chem_computer/11_scipy%E3%81%AE%E5%9F%BA%E6%9C%AC%E3%81%A8%E5%BF%9C%E7%94%A8/11.html

表示はこんな感じで。nlmによってLmaxの値とpsi2maxを割る値を変えないとうまく行かない。

Nmax=100
Lmax=5
Lmin=-5
nlm=[2,1,0]
L = (Lmin, Lmax, Nmax)
psi2=CalcPsi(nlm,L)
x=np.zeros(Nmax**3)
y=np.zeros(Nmax**3)
z=np.zeros(Nmax**3)
c=np.zeros(Nmax**3)

n=0
psi2max=np.max(psi2)
for i in range(len(xi)):
    for j in range(len(yi)):
        for k in range(len(zi)):
i          f np.abs(psi2[i,j,k])>psi2max/8:
              x[n]=xi[i]
              y[n]=yi[j]
              z[n]=zi[k]
              c[n]=psi2[i,j,k]
              n+=1

fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='3d')
ax.scatter3D(x,y,z,c=c,s=1,alpha=0.5,cmap=cm.Spectral)
plt.show()

 

結果:

さっきのサイトに出ていたnlm=653

Hydrogen_653

nlm=430

Hydrogen_430

nlm=320

Hydrogen_320png

nlm=210

Hydrogen_210png

うーん、もう一声かな、、、黄色いのを何とか消さないとうまくない。

2020年6月24日 (水)

Python+NumbaでWikipediaのList of chaotic mapsを1つずつ描く(その5) Cliford attractor

今日やるのは、、、

https://en.wikipedia.org/wiki/List_of_chaotic_maps#cite_note-14

このCliford atractorです。

http://paulbourke.net/fractals/clifford/

こんな式でa,b,c,dがパラメータ。

x[i+1]=np.sin(a*y[i])+c*np.cos(a*x[i])
y[i+1]=np.sin(b*x[i])+d*np.cos(b*y[i])

リンク先に乗っていたパラメータを描いてみた。

Cliford01 Cliford02 Cliford04 Cliford03

 

コードはこんな感じ。

import matplotlib.pyplot as plt
import numpy as np
from numba import jit

a = -1.4
b = 1.6
c = 1.0
d = 0.7

a = 1.6
b = -0.6
c = -1.2
d = 1.6

a = 1.7
b = 1.7
c = 0.6
d = 1.2

a = -1.8
b = -2.0
c = -0.5
d = -0.9

@jit
def cliford_calc(x,y,N):

    for i in range(N-1):
        x[i+1]=np.sin(a*y[i])+c*np.cos(a*x[i])
        y[i+1]=np.sin(b*x[i])+d*np.cos(b*y[i])

N=500000


fig=plt.figure(figsize=(12,12))
plt.xlim(-2.5,2.5)
plt.ylim(-2.5,2.5)



x=np.zeros(N)
y=np.zeros(N)
x[0],y[0]=0.0, 0.0
cliford_calc(x,y,N)
plt.scatter(-x,y,s=1,c='black',alpha=0.05)


plt.show()

 

2020年6月23日 (火)

Python+NumbaでWikipediaのList of chaotic mapsを1つずつ描く(その4) Ikeda Map

今日は

https://en.wikipedia.org/wiki/List_of_chaotic_maps

の中の、Ikeda Mapです。

https://en.wikipedia.org/wiki/Ikeda_map

t[i]=0.4-6./(1+x[i]**2+y[i]**2)
x[i+1]=1+u*(x[i]*cos(t[i])-y[i]*sin(t[i]))
y[i+1]=u*(x[i]*sin(t[i])+y[i]*cos(t[i]))

の形式を使って計算する。Wikipediaに乗っているu=0.918で初期値を振ってみた。

こっちは狭い領域で有名な形。

Ikeda_u0_918_small

これは広い領域でWikipediaのトップの絵に相当するもの。

 

Ikeda_u0_918_big

 

リストはこちら。

import numpy as np
import matplotlib.pyplot as plt
from numba import jit


up=0.918

@jit
def ikeda_calc(x,y,N):

    for i in range(N-1):
        t=0.4-6./(1+x[i]**2+y[i]**2)
        x[i+1]=1+up*(x[i]*np.cos(t)-y[i]*np.sin(t))
        y[i+1]=up*(x[i]*np.sin(t)+y[i]*np.cos(t))



N=500
xmax=10
xmin=-10
ymax=10
ymin=-10
Nx=100
Ny=100

fig=plt.figure(figsize=(10,10))
plt.xlim(-10,10)
plt.ylim(-10,10)


for i in range(Nx):
    for j in range(Ny):
        x=np.zeros(N)
        y=np.zeros(N)
        x[0],y[0]=xmin+(xmax-xmin)*i/Nx, ymin+(ymax-ymin)*j/Ny
        ikeda_calc(x,y,N)
        plt.scatter(x,y,s=1,c='black',alpha=0.01)


plt.show()

2020年6月22日 (月)

Python+NumbaでWikipediaのList of chaotic mapsを1つずつ描く(その3) Bogdanov map

さて今日は

https://en.wikipedia.org/wiki/List_of_chaotic_maps

の中の、

Bogdanov mapだ。

https://en.wikipedia.org/wiki/Bogdanov_map

x(n+1)=x(n)+y(n+1)

y(n+1)=y(n)+ε*y(n)+k*x(n)*(x(n)-1)+μ *x(n)*y(n)

という写像。っても何やってもWikipediaのような図にならない、、、

ということで

https://arxiv.org/ftp/chao-dyn/papers/9402/9402006.pdf

も参考にして初期値をいろいろ変えてプロットしてみる。

Bogdanov_00_k_12_00

おお、それっぽい図になった。パラメータはε=0,k=1.2,μ=0.0。

別のパラメータはε=0.0001,k=1.44,μ=-0.1。

Bogdanov_00001_k_144_01

 

リストはこちら:

import numpy as np
import matplotlib.pyplot as plt
from numba import jit


ep=0.0
kp=1.2
up=0.0

@jit
def Bogdanov_calc(x,y,N):

    for i in range(N-1):
        y[i+1]=(1+ep)*y[i]+kp*x[i]*(1-x[i])+up*x[i]*y[i]
        x[i+1]=x[i]+y[i+1]



N=300
xmax=2
xmin=-2
ymax=2
ymin=-2
Nx=100
Ny=100

fig=plt.figure(figsize=(10,10))
plt.xlim(-0.5,2)
plt.ylim(-1,1)


for i in range(Nx):
    for j in range(Ny):
        x=np.zeros(N)
        y=np.zeros(N)
        x[0],y[0]=xmin+(xmax-xmin)*i/Nx, ymin+(ymax-ymin)*j/Ny
        Bogdanov_calc(x,y,N)
        plt.scatter(x,y,s=1,c='black',alpha=0.05)


plt.show()

 

 

 

 

2020年6月21日 (日)

特別定額給付金が支給された。振込より後にはがきで連絡がくる。

これが通知書。振込より後に送られてきた。

20200621-152551

 

2020年6月20日 (土)

6月19日に公開された厚生労働省の新型コロナウイルス接触確認アプリ(COCOA) をiPhoneにインストールした。

厚生労働省がさっき公開した新型コロナウイルス接触確認アプリ(COCOA) COVID-19 Contact-Confirming Applicationを早速インストールした。

https://www.mhlw.go.jp/stf/seisakunitsuite/bunya/cocoa_00138.html

Cocoaでアプリを検索しても出てこなかったので、上のリンクのQRコードから行くのが確実。

立ち上げた画面はなぜかアマビエ、、、

20200619-175712

設定はほとんどない(Bluetoothオンとプッシュ通知だけ)で、後はこのアプリの使用を了承するかどうかを

何回も聞かれる。個人情報(メール、名前など)は入力しないようだ。

 

20200619-175520

20200619-175537

しかし陽性になった人が本当にちゃんと連絡してくれるのかと、アプリの使用人数がスマホ持ってる人の80%くらいないと意味ないんじゃ、、、とかいう不安はあるし、位置情報も個人情報も記録しないでただ

「お前は濃厚接触者になった。どこで、とも誰が、ともいうことはできない」

という通知がくるのはちょっと恐怖、、、

(参考)

AppleとGoogleのAPI

https://www.apple.com/covid19/contacttracing

https://www.google.com/covid19/exposurenotifications/

COCOAの仕様

https://cio.go.jp/node/2613

  

2020年6月19日 (金)

円内部(2次元)の一様乱数でアホなミスをする、、、球内部(3次元)のときだけ角度に気を付けると思っていた、、、半径もか、、、Pythonで図示して確認する。

ちょっとしたモンテカルロシミュレーション(ってほどでもないが)で円の内部の一様分布している乱数が必要となった。

ああ、これは3次元はヤコビアンとかで気を付けないと、2次元はそのままでよかったな、と思って

x=r cosθ

y=r sinθ

で[0,1)の一様分布する乱数u,vを使ってr=u, θ=2πvとして計算したら、、、

Circle_uniform_random1

あれ?真ん中が濃いぞ?

・・・ああ、2次元でも半径は変換しないとダメか! r dr dθ = d(r^2/2) dθなので

r=√u としないとダメだった、、、

例えばこんなところ。

https://qiita.com/aa_debdeb/items/e416ae8a018692fc07eb

変換すると、、、

Circle_uniform_random2

ちゃんと一様っぽくなった。こんなことまで忘れていてダメダメ、、、

ついでに3次元も。

http://be.nucl.ap.titech.ac.jp/~koba/cgi-bin/moin.cgi/%E7%AD%89%E6%96%B9%E7%9A%84%E3%83%93%E3%83%BC%E3%83%A0

なにも考えず

x=r cosφ sinθ

y=r sinφ sinθ

z=r cosθ

で、[0,1)の一様分布する乱数u,vを使ってr=u, φ=2πv, θ=πwとして計算したら、、

Sphere_uniform_random1

案の定、θ=0,πの方向で濃くなる。これは

sinθ dθ dφ = d(-cosθ) dφなので

φ=2πv

θ=arccos(2*w-1)

これは覚えていたのでやってみたが、

Sphere_uniform_random2

真ん中に集まる、、、のも半径の変換を忘れているから。

実際は r^2 sinθ dr dθ dφ = d(r^3/3) d(-cosθ) dφ

なので

r=∛u

φ=2πv

θ=arccos(2*w-1)

とすればようやく一様な球に。

Sphere_uniform_random3

※表示はPython+Matplotlibを使ってます。

2020年6月18日 (木)

Python+NumbaでWikipediaのList of chaotic mapsを1つずつ描く(その2) Tinkerbell Map。あのティンカーベルに似ている?

ネタがないのでWikipediaのList of Chaotic Mapsをいろいろ見ていくシリーズ。

https://en.wikipedia.org/wiki/List_of_chaotic_maps

 

今日はTinkerbell Mapだ。

https://en.wikipedia.org/wiki/Tinkerbell_map

x[i+1]=x[i]**2-y[i]**2+a*x[i]+b*y[i]
y[i+1]=2*x[i]*y[i]+c*x[i]+d*y[i]

こんな写像。これはとても簡単に

Tinkerbell

 

綺麗に描ける。うーん、ティンカーベルかー。

前にも書いたけど阪神高速のロゴに見える、、、

Hanshinkosoku

 

Pythonのコードはこちら:

import numpy as np
import matplotlib.pyplot as plt
from numba import jit

a=0.9
b=-0.6013
c=2.0
d=0.50

@jit
def tinker_calc(x,y,N):

    for i in range(N-1):
        x[i+1]=x[i]**2-y[i]**2+a*x[i]+b*y[i]
        y[i+1]=2*x[i]*y[i]+c*x[i]+d*y[i]

N=20000


fig=plt.figure(figsize=(10,10))
plt.xlim(-1.6,0.6)
plt.ylim(-1.6,0.6)



x=np.zeros(N)
y=np.zeros(N)
x[0],y[0]=-0.72,-0.64
tinker_calc(x,y,N)
plt.scatter(x,y,s=1)


plt.show()

2020年6月17日 (水)

Python+NumbaでWikipediaのList of chaotic mapsを1つずつ描く(その1) Gingerbreadman map。クッキーよりゼットンかキングジョーに見える。

コロナ対応で在宅勤務だし休みもどこにも行けないし、でとんでもなくこのブログで書くことがなくなっている、、、

こまったときのカオス、ということでこのWikipediaのリストに載っているものをやってみよう。

https://en.wikipedia.org/wiki/List_of_chaotic_maps

まずは、Gingerbreadman Map。

https://en.wikipedia.org/wiki/Gingerbreadman_map

x(n+1)=1 - y(n) + |x(n)|
y(n+1)=x(n)

という簡単だが面白い図形になる。

こんなの。こんな単純なのに対称形が面白いです。

Ginger01

あまり日本ではなじみがないが欧米でこういうクッキーを焼くところから来ている。

Ginger

しかし私には正面にすると、ゼットンかキングジョーにしか見えない、、、

Ginger03

初期値を変えて重ねると面白い形になる。

Ginger02

 

Pythonのコードはこんな感じ。

 

import numpy as np
import matplotlib.pyplot as plt
from numba import jit

 

@jit
def ginger_calc(x,y,N):

    for i in range(N-1):
        x[i+1]=1 - y[i] + np.abs(x[i])
        y[i+1]=x[i]


N=20000
xmax=20
xmin=-20
Nx=100

fig=plt.figure(figsize=(10,10))
plt.xlim(-20,20)
plt.ylim(-20,20)


for i in range(Nx):
    x=np.zeros(N)
    y=np.zeros(N)
    x[0],y[0]=xmin+(xmax-xmin)*i/Nx, -2.1
    ginger_calc(x,y,N)
    plt.scatter(x,y,s=1)

plt.show()

 

 

2020年6月16日 (火)

高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2020年6月16日 Microwave Magazine 6月号で誘電体共振器アンテナ(DRA)が出てるがあと数か月で有名になる、、、Microwave Journal 6月号は5G Sub 6GHzのアンテナシステム、FractusのVirtual Antenna, UWBでキーレスエントリーを安全に、など。

IEEE Microwave Magazine6月号でこの記事が出ている。

Microwavemagazine202006 Bd3c8b43446192d00e2955ecda305fa9
Applications in Filters and Antennas
一件あまりメジャーじゃなさそうだが、、、とある理由でもう数か月で話題になると思う。9月以降かな。
Microwave Journalで面白いのはこれかな。

Semiconductor Trends in Sub-6 GHz 5G Networks

https://www.microwavejournal.com/articles/34079-semiconductor-trends-in-sub-6-ghz-5g-networks?page=2

Junecover F5

この記事も興味深い。 FratctusのVirtual Antenna。

Optimizing Wireless Connectivity with Embedded Antennas

https://www.microwavejournal.com/articles/34113-optimizing-wireless-connectivity-with-embedded-antennas

最近iPhoneに乗ったけど使われていないUWBのキラーアプリになるのでは、というNXPのこの記事。 

UWB enhances security of passive keyless entry

https://www.edn.com/uwb-enhances-security-of-passive-keyless-entry/

2020年6月15日 (月)

福家警部補の追求を読んだ。福家さんの身体能力に驚く、、、

今回は

・未踏峰への夢を息子に託す登山家が、後援を止めようとする会社の相談役を撲殺(未完の頂上)

・動物を愛するペットショップ経営者が、動物を金儲けの道具としか見ていない義理の弟を殺害(幸福の代償)

の2篇。

倒叙ものはやっぱり面白く(コロンボ大好き)、細かいところが気になるところを徹底的に追求し、理詰めで追い詰めまた情にも訴えかける福家さんがいつもかっこいい。誰かに会うたびに(1人を除いて、、、)みんな幸せになるのもいい。

が、今回の登山の話でのあまりの身体能力の高さ(オリンピック狙える?)と動物の話での意外過ぎる弱点が分かってびっくり。

年齢もわかるし、しかもはっきり寝てないことが判明したり!そこも面白かった。

20200611-195549

2020年6月14日 (日)

今日は私の誕生日、ISC(Inverse Symbolic Calculator)でこの日が何の数式か調べてみる。

GoogleとTwitterは誕生日になるとそれぞれロウソクと風船がトップページになる。

202006141 202006142

ほぼ毎年やってますが、ISC(Inverse Symbolic Calculator)で今日が何の数式か見てみよう。

http://wayback.cecm.sfu.ca/projects/ISC/ISCmain.html

まず2.0200614は

Your input of 2.0200614 was probably generated by one
the following functions or found in one of the given tables.
Answers are given from shortest to longest description

Sum(1/(b**n*P(n)),n=1..infinity), P(n) : 2nd and 3rd degree polynomials.
2020061433073306 = sum(1/(5^n*(4*n^3-22*n^2+59*n-40)),n=1..inf)

 

次は0.6142020は

Your input of .6142020 was probably generated by one
the following functions or found in one of the given tables.
Answers are given from shortest to longest description

Sum(1/(b**n+P(n)),n=1..infinity), P(n) : 2nd and 3rd degree polynomials.
6142020691893405 = sum(1/(5^n+(3/2*n^3-n^2+7/2*n+18)),n=1..inf)

 

おお、結構似た数式になった。カシオの高精度計算サイトkeisan.casio.jpで実際に計算してみよう。

1 , 0.2 ,                                0.037037037037037037
2 , 0.201818181818181818 , 0.0542784163473818646
3 , 0.201988394584139265 , 0.05968382175278727
4 , 0.202004394584139265 , 0.0610406738559080299
5 , 0.202005955559749021 , 0.0613416067478731217
6 , 0.202006121362857829 , 0.0614042948120656994
7 , 0.202006140553262627 , 0.0614170121190487726
8 , 0.202006142941322328 , 0.061419567213267872
9 , 0.202006143256399252 , 0.0614200789350222384
10 , 0.20200614329997372 , 0.0614201813197886956
11 , 0.202006143306234802 , 0.0614202017989783882
12 , 0.202006143307163179 , 0.0614202058949363114
13 , 0.202006143307304493 , 0.0614202067141341706
14 , 0.202006143307326485 , 0.0614202068779740636
15 , 0.202006143307329973 , 0.0614202069107420583
16 , 0.202006143307330536 , 0.061420206917295658
17 , 0.202006143307330627 , 0.061420206918606378
18 , 0.202006143307330643 , 0.061420206918868522
19 , 0.202006143307330645 , 0.0614202069189209508
20 , 0.202006143307330646 , 0.0614202069189314366
21 , 0.202006143307330646 , 0.0614202069189335337
22 , 0.202006143307330646 , 0.0614202069189339532
23 , 0.202006143307330646 , 0.061420206918934037
24 , 0.202006143307330646 , 0.0614202069189340538
25 , 0.202006143307330646 , 0.0614202069189340572
26 , 0.202006143307330646 , 0.0614202069189340579
27 , 0.202006143307330646 , 0.061420206918934058
28 , 0.202006143307330646 , 0.061420206918934058
29 , 0.202006143307330646 , 0.061420206918934058
30 , 0.202006143307330646 , 0.061420206918934058

なるほど確かにこの式で計算できた。

2020年6月13日 (土)

雨のときは走ると歩くのどちらが濡れる?(雨に走れば公式)をカシオの高精度計算サイトkeisan.casio.jpに自作式としてUP!

数学セミナー2020年6月号に面白い記事が出ていた。

それは伊藤大雄さんの記事「雨に走れば公式 雨の中を歩くと走るで濡れる量の違い」だ。

そこに出てきた雨に走れば公式(Runnin' in the Rain Formula, RiR)が、今まで見た類似の解析に比べて

とても簡単な形をしていて面白いなと。

こんな形です。

Rir

これをカシオの高精度計算サイトkeisan.casio.jpに自作式としてUPした。

これです。

 雨のときは走ると歩くのどちらが濡れる?(雨に走れば公式)

説明文:

強い横風がないとし、霧雨でない雨の中を普通の体形の人間が、普通の速さで歩くときの濡れる量を1としたとき、そのx倍の速さで歩くまたは走る場合の濡れる量がRiR(x)です。

Rir2

なるほど少しでも走ると濡れなくなるが、よっぽど速く走ってもあまり効果はない、、、ということで直感と合ってる。

 

2020年6月12日 (金)

マンデルブロ集合z→z^2+cで、z^n+cと拡張してnを整数以外も選んで(n=1.5~10)GIFアニメにしてみる(Python+numba)

さっきTwitter見てたらマンデルブロ集合z→z^2+cで、z^n+cと拡張してnを整数以外も選んでアニメーションにしているのを見た。

これは!と思って自分でもやってみた。

プログラムは以前作った

Python+numbaでマンデルブロ集合を描く(今更)、、、colormapsの色を確かめるために!

の2のところを変えているだけ。

ではこちら。ちょっと動きを面白くしてみた。ちょっと癒される、、、

(クリックすると始まります)

Mandel

 

2020年6月11日 (木)

サザエさんじゃんけんが29年で初めて5連続パーというのが話題だが、ではそんなことが起きる条件付き確率をモンテカルロシミュレーションする(Excel VBA+メルセンヌツイスタで)。

サザエさん史上初、じゃんけんで「5回連続」同じ手を出す大事件 「4連続」すら29年間で3度目 

と言うのを見た。

まあ全くランダム(1/3の確率)で出るわけではなくて、人間がやるとできるだけ手を連続にしないようにしてしまうから。

こういうお話もあるし。

 

 

ここではずっと簡単に、前の手を見たときに同じ手を選ぶ条件付き確率があるとしてみましょう。

Aがある目になって、Bが同じ目になる確率を

P(B|A)

として、別の目になる確率は残りの半分ずつ、にします。

計算は

メルセンヌツイスタの高速版(SFMT)がExcel VBAで簡単に使える!

を使ってモンテカルロシミュレーションしてみる。

29年間毎週あるとして大体1500回の試行回数。本来は初期値を変えて何回もやるところですが、

今日は時間がないので1回のみ。それでも結構面白い結果に。

これです。※プログラムにバグがあったので修正しました。

Sazae_janken_20200611211501

だいたいP(B|A)<1/7くらいで29年間で5回連続が1回だけくらいに。

通常の2倍以上ぞろ目は出にくくしている感じ?

 

 

2020年6月10日 (水)

高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2020年6月10日 LGの2画面スマホの画面間通信はミリ波のKeyssa、AkoustisのWi-Fi 6Eに向けたBAWフィルタ、超小型のルビジウム原子時計、QualcommのマルチSIMの勧め、など。

これはなぜか知らなった。V60からじゃなくて前のV50からKeyssaの60GHzの超近距離通信がディスプレイをつなぐのに有線を使わなくていいように使われている。

※Samsungはものすごく変なフレキ使っている。

LG’s Dual Screen wireless technology revealed: It’s mmWave

https://www.xda-developers.com/lg-dual-screen-wireless-technology-revealed-mmwave/

Keyssaのはこんなの。

http://www.keyssa.com/technology/

https://www.anandtech.com/show/8858/a-quick-look-at-keyssa-contactless-usb-30

Dsc_2120_678x452

 

※Samsungのはこれを参照。

https://www.ifixit.com/Teardown/Samsung+Galaxy+Z+Flip+Teardown/131002

次はWiFi6に向けたBAWフィルタ。

Akoustis Announces First 5.5 GHz XBAW Filter and Sample Shipments for WiFi 6E

https://ir.akoustis.com/press-releases/detail/179/akoustis-announces-first-5-5-ghz-xbaw-filter-and-sample

50.8 x 50.8 x 19.5 mmサイズのルビジウム原子時計。

Orolia Delivers Its First Low SWaP-C Miniaturized Rubidium Oscillator

https://www.microwavejournal.com/articles/34063-orolia-delivers-its-first-low-swap-c-miniaturized-rubidium-oscillator

ちなみに防衛産業で重視されるSWaP-Cとは、

SWaP-C (Size, Weight and Power -  Cost)

です。

Mro50400

GSAのレポートが出ました。

Fixed Wireless Access Snapshot June 2020

https://gsacom.com/

QualcommのマルチSimの勧め。

5G global multi-SIM for the ultimate user experience

https://www.qualcomm.com/news/onq/2020/06/08/5g-global-multi-sim-ultimate-user-experience

 

 

2020年6月 9日 (火)

露 天神社(お初天神)でお参り。

御霊神社の後はこちらへ。しかし来るたびに思うのだが心中の舞台が恋人の聖地でいいのだろうか。

20200317-151416 20200317-151442 20200317-151451 20200317-151456 20200317-151500 20200317-151539

2020年6月 8日 (月)

高校事変VIIを読んだ。兵庫県出身の私には身につまされる話、、、今回の舞台は甲子園!

これが書かれたときにはコロナウイルスの影響がここまで大きくなるちょっとだけ前だったらしくあれ?と思うところはありますがそれでも最新の話題をいつも取り入れる松岡さん。

最初は過去の結衣がセンバツの応援に行ったときに起きた事件から始まり、、、

前話が米軍基地での戦闘だったので、ちょっとスケールダウン?と思ったら全然違った。

それは最初だけの話。それから半年以上たち、また甲子園で起きる事件とは、、、

前作よりさらにスケールアップしていた!

それはそうと、甲子園が舞台なこともあり、兵庫県のあの組織の腐敗が出てきて、兵庫県出身の私はトホホな気持ちに、、、

次作でとうとうあの親子と決着が?という期待させる終わり方でした。次も読もう。

20200605-175819

 

2020年6月 7日 (日)

なか卯で7種チーズの親子丼(大盛)を食す。以前に食べた4種チーズよりも親子丼にはマッチしているように思った。

ずっと前に4種チーズの親子丼を食べたことがあったけど、その時はチーズが強すぎる?という感想。

今回、新しくなってかつ7種に増えたもの、どんな変化が?と思って食べましたが、これはかなり親子丼と

マッチしていて美味しくなってる。改良されてるみたいですね。

ブラックペッパーがついてくるのですが、欲を言えばこれ5袋くらいほしい感じ、、、たっぷりペッパーかけられたらもっと美味しくなりそう。

20200607-130243 20200607-130433 20200607-130432 20200607-130434

2020年6月 6日 (土)

松屋でうまトマチーズハンバーグ定食(ポテサラ追加、ご飯大盛)を食す。チーズが異常に伸びて切れないが美味しい!

毎回復活するたびに食べているうまトマハンバーグ定食。今日はチーズも追加で。

このチーズがすごく伸びて全く切れないくらい。スプーンでハンバーグを運びながら箸で切る、、、が箸にまとわりつき手で切る、、、

がそれがまた美味しさを倍増させる。これは美味しいです。

20200606-142949 20200606-142947 20200606-142944

2020年6月 5日 (金)

高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2020年6月4日  Qualcommの5GはSub6GHzよりミリ波が速いという話&ドコモのLTE-Mに向けた9205、高速 DDR4でViaがスタブになる、 NXPのMIFARE DESFire EV3、アンリツの220GHzネットアナ、など。

まあこれはそうじゃないと使う意味がないからという。

mmWave 5G devices powered by Qualcomm Snapdragon can deliver 4x faster 5G

https://www.qualcomm.com/news/onq/2020/06/01/mmwave-5g-devices-powered-qualcomm-snapdragon-can-deliver-4x-faster-5g

この発表も出てた。

Qualcomm Completes Validation of Qualcomm 9205 LTE Modem on NTT DOCOMO’s Low-Power Wide-Area Networks

https://www.qualcomm.com/news/releases/2020/06/03/qualcomm-completes-validation-qualcomm-9205-lte-modem-ntt-docomos-low-power

これは高周波の人には常識だけど、高速デジタルの人も高周波の設計を学ばないといけない時代。

Signal Integrity Characterization of Via Stubs on High-Speed DDR4 Channels

https://www.signalintegrityjournal.com/articles/1731-signal-integrity-characterization-of-via-stubs-on-high-speed-ddr4-channels

Thumb_bd

NXPの新しい非接触通信(NFC)規格の発表。FeliCaとは違うけど、全部同じにしてどこの国でも

NXP Introduces MIFARE DESFire EV3 IC, Ushers In New Era of Security and Connectivity for Contactless Smart City Services

https://media.nxp.com/news-releases/news-release-details/nxp-introduces-mifare-desfire-ev3-ic-ushers-new-era-security-and

こういうのは本当にコネクタ部分が壊れやすいんですが大丈夫なんだろうか不安(昔いっぱい壊した、、、)

VNA Spans 70 kHz to 220 GHz in a Single Sweep

https://www.mwrf.com/technologies/test-measurement/article/21132942/vna-spans-70-khz-to-220-ghz-in-a-single-sweep

Ms7838g_m_rf_figure_15ed6ae0059e90

2020年6月 4日 (木)

六甲八幡神社でお参り。

阪急の駅を降りてすぐです。

20200228-153559 20200228-153804 20200228-153717 20200228-153825 20200228-153941 20200228-153905 20200228-153947 20200228-154001

2020年6月 3日 (水)

すき家でハニマスレタス牛丼(大盛)をいただく。甘めのソースが結構牛丼と合うと思ったり。

私、パンダエクスプレスの甘いオレンジチキンがめちゃくちゃ好きなので、辛いものすきだけど甘いソース+肉もありと思ってます。

これは結構甘い(ハニーだし)ソースとアクセントでクルミがあるのでかなり美味しかった。

でも万人向けではないかも?

20200531-131203 20200531-131201

2020年6月 2日 (火)

高校事変VI(松岡圭祐さん)を読んだ。沖縄の基地でアパッチ2機に襲われてどう切り抜けるか、、、

コロナウイルスのために外出を自粛していて、本屋さんにもほとんど行ってなかったのでもうVIだけでなくてVIIまで出てるのを知った。

まずはVIから読んでみた。

今回はコロナウイルスがまだここまで深刻じゃなかった(触れられてはいる)段階なので、沖縄に修学旅行に行ったりしている。

松岡さんでもここまでの状況は予測できなかったんだな、、、

で、今回の敵は貧困家庭を食い物にする反社会的勢力と、民間軍事会社。特に米軍基地の中の戦いがメインでアパッチ2機をどうやって倒すか、がポイント(いやそんなうまく行くはずが、、、なんて突っ込むのは野暮で面白いからいいのだ)。

さらにその敵を倒してもページ数かなりあるな?と思ったらもう一つ山場が。これはてんこ盛り。

でもVまでの敵が、、、ということなのでここから新展開になるのかな?VIIも読もう。

20200530-183015

 

2020年6月 1日 (月)

「ブルーローズは眠らない」を読んだ。いろんな要素てんこ盛りのミステリだ!

ジェリーフィッシュは凍らない、が面白かったので期待して読んだがこれも面白かった!

あらすじは

「科学者と牧師が同時期に独立に青いバラを作ったというニュースが流れる。マリアと蓮はジェリーフィッシュ事件で知り合った刑事に、科学者・テニエル博士とクリーブランド牧師を調べてほしいと理由も語らず頼まれる。ところがテニエル博士が施錠されてバラの蔓が壁と窓を覆い、出入りが不可能な密室状態の温室で切断された首だけの状態で殺されてしまう。

密室はどうやって作られたのか?また青いバラはこの事件と関係するのか?」

というもの。科学ネタ(DNA)もちりばめられているのですが、そこはマリアが生徒役となっていろいろ聞いて説明はわかりやすい。

2つの関連しそうな事件が交互に同時並行に起こっているのですが、もっとシンプルなものかと最初思っていて、

(単に名前が違うので世代が違うだけ)

そんなシンプルなものではなかった。最後の最後まで誰と誰が対応しているのかがわからない○○トリックも加わった○○○○トリック。

でもあの呼び方はいいのかとおもったら、

https://family.disney.com/baby-names/unisex-names/meaning-of-frankie-girl/

なるほど、そういう場合にも使うのか。私には某国民的コミックに出てくるあの人しか、、、

20200529-173924

 

« 2020年5月 | トップページ | 2020年7月 »

最近の記事

2021年4月
        1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30  
フォト
無料ブログはココログ