四川餐館（京都ポルタ店）で 陳麻婆豆腐（激辛）を食す。

ご飯おかわり自由でお得。陳麻婆豆腐以外にも副菜がいろいろついていて美味しい。

激辛にすると、辛くていいんですが、ちょっと苦味が出てしまうのでおいしく食べるには大辛くらいに抑えたほうがいいかも（って前も思ったが学習しない、、、）。

すき家でニンニクの芽牛丼＋キムチトッピング（特盛）を食す。

ニンニクの芽だけでも美味しいが、キムチでさらに辛さをプラスでさらにおいしくなる。

松のやでわらじかつ丼（ご飯大盛）＋お新香を食す。確かにカツがでかい！のでお味噌汁の蓋に一度退避していただく。

ちょっと笑っちゃうくらいカツがでかかった。ただこのままではごはんまで届かない、、、のでいったんみそ汁の蓋に退避してからいただく。

お新香がいいアクセントになって最後まで美味しくいただきました。

 

やよい軒で「特牛カルビ＆ホルモン焼定食」を食す。ホルモンが美味しい。

ホルモン大好きなのでカルビはあんまり、と思ってましたがカルビも美味しかった。もちろんホルモンも美味しい。

ただ、麺はうーん、ちょっと余計かなと思ってしまった。キャベツかなにかの野菜が良かったなとか。

「万能鑑定士Qの事件簿 0」を読んだ。高校事変をずっと読んでいるのでこれはめっちゃ新鮮に思える久しぶりの人の死なないミステリ。

松岡さんと言えば人の死なないミステリ、だったのですが、最近の作風は人が死にまくる（高校事変とか特に）ものに変わって、それはそれでとても面白いのですがちょっと寂しい感じもしていた。

万能鑑定士Qが10週年記念ということでまた昔の作風が読めてうれしかった。しかも、まだ鑑定士を初めて間もない2009年時点の初々しく悩む莉子さんが出てくる。力士シールもようやく貼られるのが認識された時期。内容もバンクシー、漢委奴国王印、ゴッホが出てきて、またグアムの探偵の登場人物も重要な役割を果たしたりと様々な仕掛けがでてくるのも楽しい。場所も東京、熱海、福岡、グアムと飛び回る。

ただ次があるとしたらもう10年後とか？

「ペンギンが教えてくれた物理のはなし」を読んだ。面白かった！バイオロギングを初めて知った。弊社みたいな超小型モジュールを作れるところが協力すれば、、、と思ったが、金にならないので社長は一切興味ないだろう、、、

これ、めっちゃ面白かったです。マグロの速度が時速80kmというのが全然違っているとか、地球規模の大移動をする生物がいるとか、目から鱗。ただ、その分析を可能にするバイオロギング技術にもものすごく興味を持った。

この本にも書いてありますが、そんなものを大企業は（金にならないので）やらないので個人経営のようなところにお願いしているというのがもったいなくて、、、

うちなんか、たぶんその個人経営のところの1/100のサイズにできる技術があるのだが、絶対にやらない。金にならんからね、、、

なんとか収益を上げるシステムを考えられないだろうか、、、と本気で思った。

紹介されているのは、

・アルゴス

人工衛星とドップラー効果を利用した低消費電力で位置を検出するシステム。

https://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%AB%E3%82%B4%E3%82%B9%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0

https://www.argos-system.org/argos/how-argos-works/

超頭いいシステムだ。これを超低消費電力にする方法ならモジュールメーカは100も考え付く（が、利益がないのでやらない、、、）

・ジオロケータ―

　照度を利用したこれもとんでもなく頭いいシステム。これも予算さえあれば、、、

http://www.birdtracker.co.uk/

・ポップアップタグ

　ああ、これもアイデアが次から次へと浮かぶ、、、が金にならないのでやれない、、、

https://camp-fire.jp/projects/view/36679

　https://www.marinecsi.org/2010/05/21/pop-up-satellite-tags/

大企業に勤めていると金のことしか言われないのでこういう面白いことに関われない、、、

なんとか収益が上がるビジネスモデル考えられないかなあ。。。

フランス絵画の精華 ルネ・ユイグのまなざし -大様式の形成と変容＠大阪市立美術館へ行ってきた。コロナ禍が始まってから4か月くらいぶりに美術館へ行けた。良かった。

月二回くらいのペースで美術館に言っていたのですが、コロナが流行し始めてから全く行く機会がなく寂しい気がしていた。

まだ終息しているわけでもないが、気晴らしに4か月くらいぶりに美術館へ行ってきた。まずは大阪市立美術館だ。

フランス絵画はこういう気がめいっているようなときにぴったりで、絵を観たなあ！という気分にさせてくれる。

もちろん入るときは体温チェックはあるし、学芸員さんたちがカードを持ってしゃべらないように、と距離をあけて、と常に注意を向けている。

それでもやはりストレスが発散された。文化に触れられずに気がめいっている方は気を付けながら行くのがお勧めです。

三宮の生田神社でお参り。

三宮神社の後はこちら。久々に来た。

三豊麺（灘店）で濃厚魚介つけ麺（大盛）を食す。

食欲がないときでもここのつけ麺は食べられる。さすがに特盛は無理で大盛りで。高菜が美味しかった。

”もういちどベートーヴェン”（中山七里）を読んだ。岬洋介がまだ若いころの司法修習生だった時代のお話ですごくさわやか。

岬洋介さんの過去が（どこかでベートーヴェンと共に）わかる作品。

司法試験をトップ合格していた岬洋介さん。司法修習生の中でも一目置かれるが、グループのメンバーの天生と親しくなる。

天生は子供のころピアノの天才と言われていたが挫折し、司法試験を受けて検事を目指している。

天才、岬さんと対比するものとして天生さんが出てきて読者は天生さんに感情移入すると思うが、、、岬さんがあまりにも屈託ない天才なので僻む気にもならないという、、、やっぱりさわやか。

一度はピアノを止めたがある出来事がきっかけで復帰することになるが、それとともにある事件（童話作家の夫婦の妻が夫を殺した事件）の真相にも気が付く岬さんがかっこいい。真相が、、、ああ、確かに、、、と読み返すと思えるようなものでした。

これは読後感がとてもいいのでお勧めです。

「博士を殺した数式」を読んだ。カオス理論の大数学者が殺され、残された数式とは？というものだが、長男が超弦理論の研究者で、著名な物理学者の醜聞（浮気癖）も書かれてたり。そして謝辞がウィッテンやシュワルツなど。

ランダウやファインマンの女癖というか、シュレーディンガーのあの話が出てきて笑ってしまった。（オッペンハイマーもなのか、、、）

パウリの微細構造定数のエピソードや、チューリングの話も出てきます。

あらすじは「シアトルで潰れかけの書店を営むヘイゼルのもとに、養祖父である天才数学者、アイザックが自殺したとの知らせが届く。アイザックからは死ぬ前にヘイゼル宛てに手紙を出しており、”命を狙われている。ある方程式をある人物に届けてほしい”との依頼があった。数学とは縁がないヘイゼルになぜその方程式が託されたのか？アイザックは自分以外にも3人家族が死ぬという。それは本当なのか？」

というもの。ヘイゼルだけでなく、アイザックの長男、フィリップ（超弦理論、M理論の研究者）やヘイゼルの兄、グレゴリー（刑事）のもとにも別の事件が起きていくのが最後にまとまるのが面白い。単に、数学や物理のエピソードをちりばめただけでなくて手紙の暗号や謎の人物の正体、兄の秘密、そして方程式が隠されていた意外な場所（というか、、、、）など飽きさせなくて面白かった。

あとがきでウィッテンやシュワルツなど著名な物理学者に謝辞があったのも面白い。取材をちゃんとして書かれたものだとわかる。ミステリの賞の候補作になったのもうなずける作品でした。

あ、ディラックのあの言葉ってほんと？と調べたらほんとだった、、

https://www.goodreads.com/author/quotes/271154.Paul_A_M_Dirac

“Age is, of course, a fever chill
That every physicist must fear.
He's better dead than living still
When once he's past his thirtieth year.”― Paul Dirac

ところでスナビトリアム教授の高慢な定数ってなに？

Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使う（その16) 時間依存のシュレーディンガー方程式のトンネル効果を計算してGIFアニメに。

時間依存のシュレーディンガー方程式も空間微分を差分化すると普通の常微分方程式になるので、ルンゲクッタ8次のDormand-Prince法で計算できる。
具体的には空間の拡散項は2次の差分
(φi+1 - 2*φi + φi-1) / (dx*dx)
で近似。中央に井戸型ポテンシャルを置いて、左から波束が入ってくるときの計算。

ただし、複素数になるのでライブラリがodeからcomplex_odeに変わる。

ずっと前にExcel VBAでやったもののPythonでの書き直し。これは”計算物理”にでてたもので、もともとはシッフの例題にあるものだそうです。

クリックするとGIFアニメが始まります。

 

ソースはこんな感じで（インデントがないけどココログにコピペすると消えるので、、、）

%matplotlib notebook
import numpy as np
from scipy.integrate import complex_ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
import numba

INTERVAL = 0.05 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

Nt=6000
Nx=2000

t0=0
dx = 0.001
dt = dx * dx / 2
xmax = dx * float(Nx) / 2
xmin = -dx * float(Nx) / 2
tmax=Nt*dt

v_width = 0.064
v0 = 0.6*((70.7 * np.pi) **2)
deltax = 0.035
x0 = -0.3
p0 = np.pi * 50
t=np.linspace(0,tmax,Nt)
x=np.linspace(xmin,xmax,Nx)
v=np.empty(Nx,dtype=np.complex)
u0=np.empty(Nx,dtype=np.complex)

for i in range(Nx):
u0[i]=np.exp(-((x[i]-x0)**2) / (4 * deltax * deltax))*np.exp(p0*x[i]*1j)/((2 * np.pi * deltax * deltax) ** 0.25)
if x[i]>=-v_width/2 and x[i] <=v_width/2:
v[i]=v0
else:
v[i]=0

sol= np.empty((Nt, Nx),dtype=np.complex)
sol[0] = u0

def Schroedinger(t, u): #odeintのときとt,xの並びが逆
f=np.zeros(Nx, dtype=np.complex)

for i in range(Nx):
if i < Nx-1:
p1=i + 1
else:
p1=0
if i > 0:
m1=i - 1
else:
m1=Nx-1

diff2 = (u[m1] + u[p1] -2.0* u[i]) / ( dx ** 2)
f[i] = (-1j)*(-diff2 + v[i]*u[i])

return f

numba_f = numba.jit(Schroedinger,nopython=True)

solver=complex_ode(Schroedinger)
solver.set_integrator('dop853')
solver.set_initial_value(u0,t0) #なぜか関数と並びが逆

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

def update(i):
plt.cla() #現在描写されているグラフを消去
plt.xlim(-1,1)
plt.ylim(0,4)
plt.plot(x,np.abs(sol[i]),c='Red',lw=3)
plt.plot(x,4*abs(v)/v0,c='blue',lw=3)
plt.xlabel("x")
plt.ylabel("Phi")
plt.title("Schroedinger equation by DOP853")
plt.grid()
# Plot
fig = plt.figure(figsize=(10,6))
ani = FuncAnimation(fig, update, frames=range(0,Nt,15), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("sh01.gif", writer = 'pillow',fps=FPS)

松屋で甘唐辛子のトロたまごろチキコンボ牛めしに七味唐辛子をたっぷりかけていただく。

七味をかけておいしさ倍増、というコピーがあるのでかけてみた。確かにそのままだとちょっと甘いのが引き締まる。

チキンは定評あるごろごろで美味しい。

すき家でニンニクの芽牛丼（特盛）をお持ち帰りでいただく。

途中で傾けたので中がばらばらに、、、でも味はいつものニンニクの芽で辛みもあってとても美味しい。単品で増量すればよかった。

Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使う（その15) 蔵本モデルを計算してGIFアニメに。

先日は蔵本・シバシンスキー方程式でしたが、今日は蔵本モデル。

dφi/dt = ωi + K/N Σsin (φj - φi )          （和はj=1 ～ N)

位相がどんどんそろっていくのがわかる。

クリックでGIFアニメが始まります。

吉野家で牛皿麦とろ御膳（特盛）＋ご飯増量をいただく。

牛タンにしようかと思ったけれど、何となく牛皿気分なので特盛で。とろろと麦飯、オクラがとてもいいコンビネーションで美味しい。

松屋で店舗限定のめちゃ大盛麻婆豆腐丼を食す。確かに餡の量がはんぱない。

ちょうど京都で店舗限定の麻婆豆腐丼があったのでいただく。辛さはあまりないが、とにかく量が多い。

七味をかけて辛くするとちょうどいい感じに。お腹いっぱいになります。

Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使う（その14) 蔵本・シバシンスキー方程式（Kuramoto-Sivashinsky equation）を計算してGIFアニメに。

蔵本先生と言えば、最近は千葉先生の研究で蔵本モデルが有名ですが、私にはこっちのほうがなじみがある蔵本・シバシンスキー方程式（Kuramoto-Sivashinsky equation）。1次元の乱流がシンプルに見られる。これはKdV方程式ができればすぐにできる。

∂u/∂t　+u*∂u/∂x = -∂^2u/∂x^2  -∂^4/∂x^4

空間微分をO(h^6)という高次の差分で近似、かつ時間微分をルンゲクッタ8次のスキームであるDormand-Prince法(DOP853)で計算。

結果はこんな感じ。クリックでGIFアニメが始まります。（最初が動きがないので面白くないがだんだん動く）

ソースはKdVと同じなので関数部分だけ。

def kuramoto(t, u): #odeintのときとt,xの並びが逆
    f=np.zeros(Nx)

    for i in range(Nx):
       if i < Nx-1:
          p1 = i + 1
       else:
          p1 = 0

       if i < Nx - 2:
          p2 = i + 2

       elif i == Nx - 2:
          p2 = 0
        else:
          p2 = 1

        if i < Nx - 3:
          p3 = i + 3
        elif i == Nx - 3:
          p3 = 0
        elif i == Nx - 2:
          p3 = 1
        else:
          p3 = 2

        if i > 0:
           m1 = i - 1
        else:
           m1 = Nx-1

        if i > 1:
           m2 = i - 2
        elif i == 1:
           m2 = Nx - 1
        else:
           m2 = Nx - 2

        if i > 2:
           m3 = i - 3
        elif i == 2:
           m3 = Nx - 1
        elif i == 1:
           m3 = Nx - 2
        else:
           m3 = Nx - 3

        d1 = (u[p3] / 9 - u[p2] + 5* u[p1] - 5 * u[m1] + u[m2] - u[m3] / 9) / (20 * dx / 3)
        d2 = (2 * u[p3] - 27 * u[p2] + 270 * u[p1] - 490 * u[i] + 270 * u[m1] - 27 * u[m2] + 2 * u[m3])/(180 * dx * dx)
        d4 = (-u[p3] +12*u[p2]-39*u[p1]+56*u[i]-39*u[m1]+12*u[m2]-u[m3])/(6*dx*dx*dx*dx)

        f[i] = -u[i] * d1 - d2 - d4

    return f

餃子の王将で7月度限定メニュー、酸辣湯麺（Bセット）を食す。そのままでも美味しいが、卓上のお酢をたっぷりかけて辣油もたらすとさらにおいしく。

酸辣湯麺は大好きでどのお店でもメニューにあれば注文する。王将の限定メニューで出たので早速注文。

かなりの具沢山（鶏肉がたっぷり）で熱々ですが、私は酸っぱい＆辛いもの好きなのでもっと、ということで

卓上のお酢をたっぷりかけて、ラー油もたらしていただく。

餃子が3個追加のBセットにした。お酢をかけても本当に熱々で美味しかった。

高周波（RF・マイクロ波・ミリ波・5G）関連ニュース2020年7月15日 Microwave Magazine 8月号はマイクロ波エンジニアのための量子コンピュータ、Microwave Journal8月号はRFツールとAI（機械学習）、SAW vs. BAW、Samsung researchの6Gホワイトペーパー、など。

まずIEEE Microwave magazineの特集は、Quantum computing for microwave engineers。

https://ieeexplore.ieee.org/xpl/mostRecentIssue.jsp?punumber=6668

なるほどこういうのを見るとRFとも深くかかわるのが分かる。

次はMicrowave Journal。

Artificial Intelligence and Machine Learning Add New Capabilities to Traditional RF EDA Tools

https://www.microwavejournal.com/articles/34234-artificial-intelligence-and-machine-learning-add-new-capabilities-to-traditional-rf-eda-tools

LTEアンテナをAIで設計改善する例。

どちらも流行りものとRFのお話でちょっと面白いというかネタが切れているというか。

あと、TF-SAW vs. BAW: RF Companies Taking Different Strategies
についてはノーコメント。

https://www.microwavejournal.com/blogs/9-pat-hindle-mwj-editor/post/34255-tf-saw-vs-baw-rf-companies-taking-different-strategies

Samsung Researchが6Gについてのホワイトペーパーだしてる。まあまだ先なので好き勝手に、、、

Samsung’s 6G White Paper Lays Out the Company’s Vision for the Next Generation of Communications Technology

https://news.samsung.com/global/samsungs-6g-white-paper-lays-out-the-companys-vision-for-the-next-generation-of-communications-technology

あとはこれに驚く、、、

アナログ・デバイセズ、アナログ半導体の主導的地位の強化に向け、Maxim Integratedを買収することを発表

https://www.analog.com/jp/about-adi/news-room/press-releases/2020/7-13-2020-analog-devices-announces-combination-with-maxim-integrated.html

すき家でチーズカルビ牛丼（大盛）を食す。ものすごく甘いタレだが、タバスコと一緒にするとなかなかおいしい。

大盛にすると追いカルビのタレが２つついてきた。これ、かけなくてもかなり濃い甘い味付けのカルビです。

かけなくてもいいかな、、とおもったけどかけて、その代わりにタバスコを大量に投入して食べた。

三種のチーズ牛丼とはやっぱりカルビの分、肉食べたなー、という感じが強かった。結構美味しかったです。

Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使う（その13) KdV方程式の時間発展にDOP853、空間差分はオリジナルのザブスキーとクルスカルのを使ってGIFアニメに。

KdV方程式

∂u/∂t　+u*∂u/∂x+δ^2 ∂^3u/∂x^3 = 0

δ=0.022

ですが、これも空間を差分化すると普通の常微分方程式系と見なせる。そこでルンゲクッタ8次のDOP853を使おう。

空間差分はオリジナルのザブスキーとクルスカルのものと同じとします。

クリックでGIFアニメが始まります。

ソースコードはこんな感じ。

%matplotlib notebook
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
import numba

INTERVAL = 0.05 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

Nt=512
Nx=256

xmax=2.0
xmin=0.0
tmin=0.0
tmax=16/np.pi

t0=0
dt = (tmax-tmin)/Nt
dx = (xmax - xmin) / Nx
d=0.022

t=np.linspace(0,tmax,Nt)
x=np.linspace(0,xmax,Nx)
u0=np.cos(np.pi*x)
sol= np.empty((Nt, Nx))
sol[0] = u0

def kdv(t, u): #odeintのときとt,xの並びが逆
f=np.zeros(Nx)

    for i in range(Nx):
        if i < Nx-1:
            p1 = i + 1
        else:
            p1 = 0

        if i > 0:
            m1 = i - 1
        else:
            m1 = Nx-1

        if i < Nx - 2:
            p2 = i + 2
        elif i == Nx - 2:
            p2 = 0
        else:
            p2 = 1

        if i > 1:
           m2 = i - 2
        elif i == 1:
            m2 = Nx - 1
        else:
            m2 = Nx - 2

        diff3 = (-u[m2] + 2 * u[m1] - 2 * u[p1] + u[p2]) / (2 * (dx ** 3))
        diff1 = (-u[m1] + u[p1]) / (2 * dx)
        f[i] = -(u[m1]+u[i]+u[p1]) * diff1/3 - d*d*diff3

    return f

numba_f = numba.jit(kdv,nopython=True)

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

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

def update(i):
    plt.cla() #現在描写されているグラフを消去
    plt.xlim(0,2)
    plt.ylim(-1,3)
    plt.plot(x,sol[i],c='Red',lw=3)
    plt.xlabel("t")
    plt.ylabel("u")
    plt.title("KdV equation by DOP853")
    plt.grid()

# Plot
fig = plt.figure(figsize=(10,6))
ani = FuncAnimation(fig, update, frames=range(0,Nt,1), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("kdv.gif", writer = 'pillow',fps=FPS)

かつやで全力大人飯（チキンカツ丼）を食す。チキンカツ、焼きそば、から揚げが入った大人様ランチ。

これはなかなかの茶色いビジュアルで、普通は組み合わせないだろう！というものが組み合わさっていて面白い。

量もたっぷりあります。焼きそばが冷たかったのだけが残念。

ココイチでスパイスカレーTHEチキベジ（3辛大盛）を食す。特製スパイスと赤玉ねぎのピクルスがついてくる。

特製スパイスとピクルスが赤玉ねぎのピクルスがついてくるのがいいんですよね。

このスパイスをたっぷりかけるために辛さはちょっと控えめに3辛。

茄子とトマトとオクラがとてもいい食感で、またスパイシーで美味しい。

とはいえもうちょっと辛さが欲しくなって最後の方はとび辛スパイスも投入で最後まで美味しくいただく。

Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使う（その12） 誤差設定を1e-12にしてローレンツ方程式でちょっとだけ初期値が違うものを並べてGIFアニメにする。

DOP853での設定でatol, rtol（絶対誤差、相対誤差）があるがどちらも1E-12くらいにしておくと、ピタゴラスの三体問題なんかが簡単に計算できる。今回はローレンツ方程式で初期値をx0=0.1, 0.1001, 0.100001にしてみた（y0=z0=0.1は共通）。

計算結果はこちら。クリックでGIFアニメが始まります。

新型コロナウイルス、中国、日本、韓国、アメリカでの感染者数を指数関数とロジスティック関数でフィッティングした。（10/18更新）アメリカと日本は再び増加傾向。アメリカはトランプがあれだから、ですが日本も油断している人多すぎでは、、、

(全然終息しなくて記事が長くなったので再分割)

まずは日本のPCR陽性者数は厚生労働省が毎日報告してます。

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

中国、韓国、アメリカのはWikiPediaに詳しい。

https://en.wikipedia.org/wiki/COVID-19_pandemic_in_mainland_China

https://en.wikipedia.org/wiki/COVID-19_pandemic_in_South_Korea

https://en.wikipedia.org/wiki/COVID-19_pandemic_in_the_United_States

さて、感染症の微分方程式として

SIRモデル

SIRモデル（感染症の微分方程式、コロナウイルスなども）の計算をカシオの高精度計算サイトkeisan.casio.jpにUP!　4段4次のルンゲクッタ法を使用。

SEIRモデル

SEIRモデル（コロナウイルス感染の微分方程式その２）をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP!ルンゲクッタ法を使用。

の計算ができる自作式をUPしましたが、データのフィッティングには

もっと簡単なSISモデルから出てくるロジスティック曲線を使います。

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

https://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B8%E3%82%B9%E3%83%86%E3%82%A3%E3%83%83%E3%82%AF%E6%96%B9%E7%A8%8B%E5%BC%8F

7/12更新：

日本が指数関数に近づいている、、、あれだけ連日東京で患者増えればそうなるな、、、

アメリカはどうすんだろう、、、このままではもう、、、

韓国が終息したようなことを言っていたのも今は昔。中国も（発表している分だけは）微増（実際はどうだかわからんが）。

7/19更新：

日本とアメリカは完全に第二波がきたようだ。これでGoToキャンペーンとは、、、

7/26更新

完全に日本が指数関数になってる。8万人ってのもあながちない話ではない。アメリカはずっとひどい、、、

8/2更新

日本は完全に指数関数に乗ってる、、、このままではやばい。

8/9更新

日本の指数関数化が止まらない。こうみると第一波ってなんなんだという感じ。アメリカはさすがに増加スピードは落ちた？

8/16更新

日本の指数関数的増加は止まらず、韓国もとうとう急激に増えた。アメリカは、、、さすがにちょっとだけ増加は減ってるが焼け石に水。

8/23更新：

日本の指数関数的増加は変わらずだが、韓国が1週間前くらいからとんでもないペースで増えている。第三波？　

8/30更新：

日本・アメリカはそれでも少し増加減ったが、韓国がどんどん増えている。

9/6更新：

日本・アメリカ・韓国ともすこし増加速度は落ちているようだ。

9/13更新：

　どこも増加速度が落ちている中、日本は微増している感じ？　

9/20更新：

あれ？アメリカまた増えだした？

（9/27更新）アメリカと日本はまたちょっと増えだした。。。

（10/11更新）トランプまでかかったアメリカは終息の兆しなし、日本もだな、、、

（10/18更新）アメリカと日本は再び増加傾向。アメリカはトランプがあれだから、ですが日本も油断している人多すぎでは、、、

実際最近、マスクしてなかったりしててもちゃんとつけてなかったり、大声で話したりという人が多すぎる。

増山超能力師大戦争を読んだ。前作は割とほのぼの系でもあったが、今回はハード。ダークマターが超能力の原因とは、、、

前作がとても面白かったので今回も楽しみにしていた。前作でのユーモアも満載なのもよかったですが、、、

めっちゃハードというかシリアスになってる！

超能力がダークマター（DM)要因であると解析された世界。企業がDMを検出する装置などを開発している。

そんななか、DMの痕跡を消す装置を開発した責任者が失踪し、奥さんが増山の事務所に捜索を依頼する。

しかし増山超能力師事務所の面々が何者かに襲われる、、、また増山の家族も。一体なぜ技術者は失踪したのか、裏で糸を引いているのは誰か？

というもの。

重くなりそうなところを明美がうまく落としてくれてそこは面白いのですが、味方だと思っていたあの人がそういうことをしていたとは、、、

（ってネタバレなのでやめる）。そして増山の娘のアリスちゃんについても、次回作でどえらいことになりそうな予感。

（メモは事務所じゃないよね？家だよね？）

ということで次回作も楽しみ。

AX(伊坂幸太郎さん)を読んだ。恐妻家の殺し屋の連作短編集だが、途中で驚く、、、そして最後にあの伏線が効いてるとは、、、

伊坂さんの過去の作品に出てきた殺し屋も（本人か噂かはありますが）出てきます。連作短編集の、最初（AX、BEE、Crayon）は殺し屋の話とは言え家族の話も挟まれてユーモアもあったり楽しいのですが、”EXIT”で「え！」と驚く、、、これは予想もしてなかった。

一体どうなるんだ？と思ったが最後の”FINE”でまた「あれがここにくるのか！」と驚く。

これは何言ってもネタバレになりそうなのでここまで、ですが、一番驚いたのは「古山高麗雄」さんが実在の人物で本も実在ってことでした！

Python＋Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その11） 中二病用語、ブルースカイ・カタストロフィを生じるGavrilov Shilnikov modelを計算してGIFアニメに。

千葉逸人さんが命名者が中二病だと言ったので有名な？Blue-Sky catastropheを生じるGavrilov Shilnikov modelをPythonとScipyのDOP853で計算してみる。
スカラーペディアはご本人が書いているので詳しい。

http://www.scholarpedia.org/article/Blue-sky_catastrophe

こんな式です。

計算結果はこちら：クリックでアニメが始まります。

ソースはこんな感じで。

%matplotlib notebook
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter

INTERVAL = 0.001 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

def bluesky(t, x,myu,eps): #odeintのときとt,xの並びが逆

    x_dot = x[0]*(2.0+myu-10.0*(x[0]**2+x[1]**2))+x[2]**2+x[1]**2+2.0*x[1]
    y_dot = -x[2]**3-(1+x[1])*(x[2]**2+x[1]**2+2.0*x[1])-4.0*x[0]+myu*x[1]
    z_dot = (1+x[1])*(x[2]**2)+x[0]**2-eps
    return [x_dot, y_dot, z_dot]

t0=0
tmax = 1000
N=20000

#x0=[0.012277918,-2.356078578,0.018241293]
x0=[0,0,0]
solver=ode(bluesky)
solver.set_integrator('dop853',atol=1.E-12,rtol=1.E-12)
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(0.456,0.0357)

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

def update(i):

    plt.cla() #現在描写されているグラフを消去
    ax.plot(sol[:i,0], sol[:i,1], sol[:i,2],c='Blue',lw=1,alpha=0.8)
    ax.scatter(sol[i,0], sol[i,1], sol[i,2],c='Blue')
    ax.set_xlim(-1.5,1.5)
    ax.set_ylim(-1.5,1.5)
    ax.set_zlim(-1.5,1.5)
    ax.set_xlabel("X Axis")
    ax.set_ylabel("Y Axis")
    ax.set_zlabel("Z Axis")
    ax.set_title("Gavrilov Shilnikov model : Blue-sky catastrophe(DOP853)")
    plt.grid()

# Plot

fig = plt.figure(figsize=(12,12))
ax = fig.gca(projection='3d')
ax.view_init(elev=10, azim=60)
ani = FuncAnimation(fig, update, frames=range(0,len(t),100), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("bluesky.gif", writer = 'pillow',fps=FPS)

Python＋Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その10） 7体問題の昴問題(the Pleadis)

三体問題シリーズが好評だったので、調子に乗って7体問題をやってみよう。

常微分方程式の数値解法I」に載っていた7体問題で「昴問題」(the Pleadis)という名前がついている。

pxi '= ΣGmimj(xj - xi)/rij3

pyi '= ΣGmimj(yj - yi)/rij3

で、

(m1, ... , m7)=(1,2,3,4,5,6,7)

初期値は

(x1, ..., x7) = (3,3,-1,-3,2,-2,2)

(y1, ..., y7) = (3,-3,2,0,0,-4,4)

(x1', ..., x7') = (0,0,0,0,0,1.75,-1.5)

(y1', ..., y7') = (0,0,0, -1.25,1,0,0)

計算結果はこんな感じ。確かに精度よく計算するのは難しい。（クリックでGIFアニメが始まります。）

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その9） 三体問題の周期解いろいろ（3） 1223個見つかった解のうち、I.A68(0.5),I.A2(4),I.B59(0.75)

今回は、さらに見つかった1223個の解について：

1223個の新しい三体問題の周期解が見つかったということで、私もExcel VBAにルンゲクッタ8次のDOP853を移植したもので計算してみた。なかなか面白い動きをする。

初期値は
x0=[-1,0.2138410831, 0,0.0542938396,1,0.2138410831, 0,0.0542938396,0,-2*0.2138410831/0.5, 0,-2*0.0542938396/0.5]

で、m1=m2=1, m3=1/2とする。

計算結果はこちら：クリックするとGIFアニメが始まります。

先日よりさらにすごい動き。

次は、I.A2(4)

v1=0.9911981217
v2=0.7119472124
m1=1
m2=1
m3=4
x0=[-1,v1, 0,v2,1,v1, 0,v2,0,-2*v1/m3, 0,-2*v2/m3]

 

つぎは？I.B59(0.75)

v1=0.4101378717
v2=0.1341894173
m1=1
m2=1
m3=0.75
x0=[-1,v1, 0,v2,1,v1, 0,v2,0,-2*v1/m3, 0,-2*v2/m3]

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その8） 三体問題の周期解いろいろ（２）13個見つかった解のうち、クラスI.A.1のbutterfly I

先日はピタゴラスの三体問題と8の字解を計算した。

シンプレクティック8次公式のExcel VBAのプログラムソース

今回は、ついで見つかった13個の解についてやってみよう。

なんと三体問題の新しい解が13個も見つかっていた！さっそく計算して図示＆GIFアニメ化

初期値は

x0=[-1, 0.306893,0, 0.125507,1,0.306893,0, 0.125507,0,-2*0.306893,0,-2*0.125507]

で、m1=m2=m3=1とする。

クリックするとアニメが始まります。

なかなか面白い動き。

次はこちら：さらに見つかった1223個の解

1223個の新しい三体問題の周期解が見つかったということで、私もExcel VBAにルンゲクッタ8次のDOP853を移植したもので計算してみた。なかなか面白い動きをする。

ホワイトラビット（伊坂幸太郎さん）を読んだ。これは絶対だまされる！帯の全てを、疑え！で疑って読んでいたのだが全く真相が違っていた！

ちょくちょく挟まれるレ・ミゼラブルの話がめっちゃ面白く、いいアクセントになってる（伊坂幸太郎さんも楽しんで書いてるのだろう）。

誘拐をビジネスにしているグループに属している兎田は、自分の新妻を誘拐されて返してもらうにはある男（折尾、オリオオリオ）を見つけ出さなければならなくなった。折尾はオリオン座について語ると止まらない、胡散臭いコンサルタント。GPSの情報を基に忍び込んだ家では折尾が見つからない。そこにいた母子と旦那？（これは伊坂さん作品によくでる黒澤さんです。びっくりした）を拘束し、家に立てこもる。

警察に通報された兎田だが、逆に警察に折尾を探してもらおうとする。折尾が見つかり、新妻の綿子ちゃんを救い出すことができるのか？

というものですが、これは絶対真相わからないよ！テレビで見てたのがそういう場所で見てたのかみたいな。

あの紙切れがそういう伏線だったとか、、、これはだまされてスッキリする作品で、お勧めです。

ごちそうハンバーグトマトカレー（大盛）を松のやでいただく。酸味のあるルーが美味しく、ごちそうハンバーグがすごいボリュームで満足。

最近は売ってないが松屋のトマトカレーが大好きだったので、この松のやのトマトカレーにも期待していった。

酸味がほどよく、これもなかなかおいしい。それよりもごちそうハンバーグの量が半端ないのでもうお腹いっぱいに。

ポテサラもついていて、これはお勧め。

 

松屋でお肉たっぷり回鍋肉定食（小鉢はキムチ）を食す。ちょうどいい辛さで美味しい。

ご飯は大盛（無料）ですが、ご飯の方がなくなるほど味の濃い回鍋肉。キムチとも相性がよく、これは美味しい。

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その7） 三体問題の周期解いろいろ（１） 8の字解

昨日、ピタゴラスの三体問題が計算できるようになったので、いろいろな三体問題も初期値と質量を変えればすぐに計算できる。

有名なのは8の字解

シンプレクティック8次公式のExcel VBAのプログラムソース

ついで見つかった13個の解

なんと三体問題の新しい解が13個も見つかっていた！さっそく計算して図示＆GIFアニメ化

さらに見つかった1223個の解

1223個の新しい三体問題の周期解が見つかったということで、私もExcel VBAにルンゲクッタ8次のDOP853を移植したもので計算してみた。なかなか面白い動きをする。

といろいろある。まずは8の字解かな。

昨日のプログラムで、m1=m2=m3=1として、初期値を

x0=[ 0.97000436,-0.5 * 0.93240737,-0.24308753,-0.5 * 0.86473146, -0.97000436,-0.5 * 0.93240737,0.24308753,-0.5 * 0.86473146,0,0.93240737, 0,0.86473146]

とする。

計算結果はこちら。クリックするとGIFアニメが始まります。

続く。

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その6） ピタゴラスの三体問題を計算してFuncAminationでGIFアニメにしてみる。

昨日は最終的な軌跡を描いただけでしたが、せっかくなんでアニメにしたい。

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その5）　ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。

いろいろ調べて、なんとなくこんな感じでできるんじゃないかということでやってみた（おそらくめっちゃ無駄なことをやっているが、それはおいおい修正するとして、、、）

これ。クリックで動きます。

これならいろいろアニメにして遊べそう。他の三体問題もやってみる。

ソースはこんな感じで。

%matplotlib notebook
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter

INTERVAL = 0.0005 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

def pythagoras(t, x, m1, m2, m3): #odeintのときとt,xの並びが逆
    f=np.zeros(12)

    qx1=x[0]
    vx1=x[1]
    qy1=x[2]
    vy1=x[3]

    qx2=x[4]
    vx2=x[5]
    qy2=x[6]
    vy2=x[7]

    qx3=x[8]
    vx3=x[9]
    qy3=x[10]
    vy3=x[11]

    Ra=np.sqrt((qx2-qx1)**2 + (qy2-qy1)**2)**3
    Rb=np.sqrt((qx3-qx1)**2 + (qy3-qy1)**2)**3
    f[0]=vx1
    f[1]=m2*(qx2-qx1)/Ra + m3*(qx3-qx1)/Rb
    f[2]=vy1
    f[3]=m2*(qy2-qy1)/Ra + m3*(qy3-qy1)/Rb

    Ra=np.sqrt((qx1-qx2)**2 + (qy1-qy2)**2)**3
    Rb=np.sqrt((qx3-qx2)**2 + (qy3-qy2)**2)**3
    f[4]=vx2
    f[5]=m1*(qx1-qx2)/Ra + m3*(qx3-qx2)/Rb
    f[6]=vy2
    f[7]=m1*(qy1-qy2)/Ra + m3*(qy3-qy2)/Rb

    Ra=np.sqrt((qx1-qx3)**2 + (qy1-qy3)**2)**3 
    Rb=np.sqrt((qx2-qx3)**2 + (qy2-qy3)**2)**3
    f[8]=vx3
    f[9]=m1*(qx1-qx3)/Ra + m2*(qx2-qx3)/Rb
    f[10]=vy3
    f[11]=m1*(qy1-qy3)/Ra + m2*(qy2-qy3)/Rb

    return f

t0=0
tmax = 100
N=10000

x0=[1,0,3,0,-2,0,-1,0,1,0,-1,0]
solver=ode(pythagoras)
solver.set_integrator('dop853',atol=1.E-12,rtol=1.E-12)
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(3,4,5)

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

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

def update(i):
    plt.cla() #現在描写されているグラフを消去
    plt.plot(sol[0:i,0],sol[0:i,2],c='Red',lw=1)
    plt.plot(sol[0:i,4],sol[0:i,6],c='Blue',lw=1)
    plt.plot(sol[0:i,8],sol[0:i,10],c='Green',lw=1)
    plt.scatter(sol[i,0],sol[i,2],c='Red')
    plt.scatter(sol[i,4],sol[i,6],c='Blue')
    plt.scatter(sol[i,8],sol[i,10],c='Green')
    plt.xlim(-8,8)
    plt.ylim(-8,8)
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.title("Pythagoras 3 body problem by DOP853")
    plt.grid()

# Plot
fig = plt.figure(figsize=(6,6))
ani = FuncAnimation(fig, update, frames=range(0,len(t),20), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("pythagoras.gif", writer = 'pillow',fps=FPS)

plt.show()

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その5） ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。

非常に計算が難しいので有名なピタゴラスの三体問題。

質量mが3,4,5の三体が辺の長さが3,4,5の直角三角形の頂点におかれている場合にどのような動きをするか、というもの。(Burrauさんが問題提起して、Szebehelyさんが数値的に解決したということです）。

詳細はこちら。

http://www.ucolick.org/~laugh/oxide/projects/burrau.html

でDOP853で相対＆絶対誤差を非常に小さく(rtolとatol)すると割と簡単に計算できる。

こんな感じ。

ソースはこちら。

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

def pythagoras(t, x, m1, m2, m3): #odeintのときとt,xの並びが逆
    f=np.zeros(12)

    qx1=x[0]
    vx1=x[1]
    qy1=x[2]
    vy1=x[3]

    qx2=x[4]
    vx2=x[5] 
    qy2=x[6]
    vy2=x[7]

    qx3=x[8]
    vx3=x[9]
    qy3=x[10]
    vy3=x[11]

    Ra=np.sqrt((qx2-qx1)**2 + (qy2-qy1)**2)**3
    Rb=np.sqrt((qx3-qx1)**2 + (qy3-qy1)**2)**3
    f[0]=vx1
    f[1]=m2*(qx2-qx1)/Ra + m3*(qx3-qx1)/Rb
    f[2]=vy1
    f[3]=m2*(qy2-qy1)/Ra + m3*(qy3-qy1)/Rb

    Ra=np.sqrt((qx1-qx2)**2 + (qy1-qy2)**2)**3
    Rb=np.sqrt((qx3-qx2)**2 + (qy3-qy2)**2)**3
    f[4]=vx2
    f[5]=m1*(qx1-qx2)/Ra + m3*(qx3-qx2)/Rb
    f[6]=vy2
    f[7]=m1*(qy1-qy2)/Ra + m3*(qy3-qy2)/Rb

    Ra=np.sqrt((qx1-qx3)**2 + (qy1-qy3)**2)**3
    Rb=np.sqrt((qx2-qx3)**2 + (qy2-qy3)**2)**3
    f[8]=vx3
    f[9]=m1*(qx1-qx3)/Ra + m2*(qx2-qx3)/Rb
    f[10]=vy3
    f[11]=m1*(qy1-qy3)/Ra + m2*(qy2-qy3)/Rb

    return f

t0=0
tmax = 100
N=10000

x0=[1,0,3,0,-2,0,-1,0,1,0,-1,0]
solver=ode(pythagoras)
solver.set_integrator('dop853',atol=1.E-12,rtol=1.E-12)
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(3,4,5)

t=np.linspace(0, tmax, N)
sol= np.empty((N, 12))
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.rcParams["font.size"] = 18
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Pythagoras 3 body problem by DOP853")
plt.grid()
plt.plot(sol[:,0],sol[:,2])
plt.plot(sol[:,4],sol[:,6])
plt.plot(sol[:,8],sol[:,10])
plt.xlim(-6,6)
plt.ylim(-6,6)

plt.show()