フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算。

今回は神経を通る電気信号をモデル化したフィッツヒュー・南雲方程式

dv/dt = v - v³/3 -w +I

dw/dt = ε(v + a - b*w)

をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算してみた。

結果はこんな感じ。スパイクが見える。

パラメータはここと同じにしました。

https://www.comsol.com/blogs/understand-the-dynamics-of-the-fitzhugh-nagumo-model-with-an-app/

ソースコードはこちら：

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

a=0.7
b=0.8
eps=0.08
I=1

def FNmodel(t, u): #odeintのときとt,xの並びが逆
    v=u[0]
    w=u[1] 

    v_dot = (v*(1.0-(v**2)/3)-w+I)
    w_dot = eps*(v+a-b*w)

    return [v_dot, w_dot]

t0=0
tmax = 200
N=5000

u0=[0.0,0.0]

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

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

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

# Plot
fig = plt.figure(figsize=(10,24))

ax1=fig.add_subplot(3, 1, 1)
ax2=fig.add_subplot(3, 1, 2)
ax3=fig.add_subplot(3, 1, 3)

ax1.set_xlabel("Time")
ax1.set_ylabel("v")
ax1.grid(True)
ax2.set_xlabel("Time")
ax2.set_ylabel("w")
ax2.grid(True)
ax3.set_xlabel("v")
ax3.set_ylabel("w")
ax3.grid(True)

ax1.plot(t,sol[:,0],c='Red')
ax2.plot(t,sol[:,1],c='Blue')
ax3.plot(sol[:,0],sol[:,1],c='Green')

plt.show()

カオスを生じる電気回路、Chua’s circuitをLTspiceで回路シミュレーションしてみる。

Chua's circuit

というのを初めて知った。Chuaさんが早稲田大学を訪れていたときに考案した回路だそうだ。

回路的には簡単にできそうなのでLTspiceで回路シミュレーションしてみた。

この辺を参考に：

Chua's Circuit

8. SIMULATING CHUA’S CIRCUIT WITH LTSPICE

では回路はこんな感じ。OPアンプはOP29にしておいた。

そして実行結果は！

おお、本当だ。回路シミュレータでこんな波形を見るのは珍しくて面白い。

高周波（RF・マイクロ波・ミリ波・5G）関連ニュース2021年2月16日 IEEE Microwave Magazineの特集はオールデジタルのRFID、Microwave JournalはＥバンド ミリ波通信に衛星や気球を使う話、アメリカの半導体企業がバイデンに投資を迫る、

Microwave magazineの特集は、オールデジタルのRFID、、、ってもアンテナやアンプ、フィルタはアナログだから、、、ってRFフロントエンドは別扱いか。

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

Microwave Journalの特集は70/80GHzといったミリ波で空から通信する話。

E-Band mmWave Technology for HAPS and LEO Satellite Systems

Holographic Radar Development

も。

CMOSで100Gbpsのオプティカルモジュレータ。

Standard CMOS Yields Integrated Photonic 100-Gb/s Optical Modulator

主要なアメリカの半導体メーカがバイデンに投資を迫る、、、

CEOs Urge President Biden to Fund Chips, Executive Order Expected

「 The CEO signatories on the letter are from AMD, Analog Devices, Cree, GlobalFoundries, Intel, Lattice Semiconductor, Marvell Semiconductor, Maxim, Micron Technology, ON Semiconductor, Qorvo, Qualcomm, Silicon Labs, Skyworks, Texas Instruments, Western Digital, Xilinx and the SIA. In addition, senior executives from Broadcom, IBM and Nvidia also added their signatures. 」だって。

TSMCのLiuさんがISSCCの基調講演で語る。TSMCはほんとうにすごいな。

ISSCC Plenary: A Bright Foundry Future

Labyrinth Chaos(迷宮カオス）を生むThomas-Rössler方程式のパラメータbを色々変えて、Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使って計算してGIFアニメ

先日、こういうのを書いた。

Labyrinth Chaos(迷宮カオス）というのがあるのか！Thomas-Rössler方程式から出る。Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使って計算

こういう方程式（Thomas' cyclically symmetric attractor　or Thomas-Rössler system）

で、bを変えると非常に面白いパターンがでる。ということで今回はGIFアニメにしてみた。

初期値も3種類変えて同時に表示してます。

こちら：

b＝０も面白いけど、途中で星型みたいになるものまた面白い。

参考リンク：

Hyperchaos & Labyrinth chaos: Revisiting Thomas-Rössler systems

Labyrinth chaos: Revisiting the elegant, chaotic and hyperchaotic walks

Thomas' cyclically symmetric attractor

「水晶振動子の等価回路計算」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP！ インピーダンスの大きさと位相がグラフ化できる。

ちょっと頼まれて作ったものがあって、それを転用します。

水晶振動子はこんな風に等価回路化できる。

回路定数は

直列容量C1
直列インダクタンスL1
直列抵抗R1
並列容量C0
負荷容量Cs

を入れるようになっている。

リンクはこちら：

 水晶振動子の等価回路計算

入力画面：

インピーダンスの大きさ|Z|のLog表示：

位相：

水晶はものすごくQが高いのでLでかくとH単位（nHとかμHじゃなくて）になるのが面白いね。

Labyrinth Chaos(迷宮カオス）というのがあるのか！Thomas-Rössler方程式から出る。Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使って計算

Thomas' cyclically symmetric attractor

というのをたまたま見ていた。そこでの脚注で知ったのですが、

こういう方程式（Thomas' cyclically symmetric attractorとかThomas-Rössler systemとかいうらしい）

で、b=0とすると迷路というか迷宮というかランダムウォークみたいになるらしい。

Hyperchaos & Labyrinth chaos: Revisiting Thomas-Rössler systems

とか

Labyrinth chaos: Revisiting the elegant, chaotic and hyperchaotic walks

なんかが参考になりそう。では計算してみる。b=0.3からb=0までやってみよう。

b=0.3

b=0.2

b=0.1

そしてb=0

うわ！本当だ。こんなシンプルな式でこんなパターンがでるのか！

ニューロンのバースト現象のモデル、Hindmarsh–Rose modelをPython＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使って計算

たまたま見てた本でHindmarsh–Rose modelというのがあるのを初めて知った。

https://en.wikipedia.org/wiki/Hindmarsh%E2%80%93Rose_model

こういう3変数の常微分方程式。

Wikipediaにパラメータの大きさも出てるので簡単に計算できそう、ということでやってみる。

Iの大きさでだいぶ振る舞いが変わるよう。

I=2のときは、きれいにバーストが繰り返されるのが見える。

I=1のとき、途中で亡くなった？

I=3のとき、バーストの間隔が長く。

なるほど結構面白いな。

渋滞を表すOVモデル（最適速度モデル）をPython＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使って計算、GIFアニメにしてみた。又吉直樹のヘウレーカが渋滞特集だったので。

先日、

2/3の又吉直樹のヘウレーカ 「どうすれば“密”を避けられる?」（渋滞学で有名な西成活裕さんがゲスト）をTwitterでリアルタイムメモしてました。その記録＆昔、渋滞の微分方程式OVモデルを計算した結果。

というのを見てました。それで昔、OVモデル（最適速度模型、Optimal Velocity model）

dV_i/dt = a [U(X_i+1 - X_i) - V_i]

dX_i/dt = V_i

ここでU(x) = tanh(x-c) + tanh(c)

をExcelのVBAで計算したことを思い出した。今やるならPython＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)かな、ということでやってみた。

クリックで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
import random
from scipy.stats import cauchy

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

Nt=10000
N=40
L=80
c=2.0
a=1.3
dx=L/N

tmin=0.0
tmax=500.0

t0=0
dt = (tmax-tmin)/Nt

t=np.linspace(0,tmax,Nt)
x0=np.empty(N*2)

def Velocity(x):
    return np.tanh(x-2)+np.tanh(2)

p=0
for i in range(N):
    x0[i]=p
    x0[N+i]=Velocity(dx)
    p+=dx
    p+=0.01*(random.uniform(0,1)-0.5)

sol= np.empty((Nt, N*2))
sol[0] = x0

def ov_model(t, x): #odeintのときとt,xの並びが逆
    f=np.zeros(N*2)

    for i in range(N):
        f[i]=x[N+i]
        if i< N-1:
            d=x[i+1]-x[i]
        else:
            d=x[0]-x[i]

        if d>L:
            d -= L

        if d<0:
            d +=L
            f[N+i]=a*(Velocity(d)-x[N+i])

    return f

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

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

    sol[k] = solver.y
    for i in range(N):
        if sol[k,i]>L:
    sol[k,i] -=L
    k+= 1

theta=np.linspace(0,2*np.pi,100)

def update(i):
    plt.cla() #現在描写されているグラフを消去
    plt.xlim(-1.1,1.1)
    plt.ylim(-1.1,1.1) 
    plt.plot(np.cos(theta),np.sin(theta),c='Black',alpha=0.2)
    plt.scatter(np.cos(2*np.pi*sol[i,0:N]/L),np.sin(2*np.pi*sol[i,0:N]/L),c='Red')
    plt.xlabel("x")
    plt.ylabel("y")
plt.title("OV model by dop853")
plt.grid()

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

「ラグランジュ補間(11点まで）」をカシオの高精度計算サイトkeisan.casio.jpに自作式としてUP！

Twitterで久しぶりにラグランジュ補間の話が出ていた。そういや自作式に作ってないな、と思って早速やってみた。

リンクはこちらです。

 ラグランジュ補間(11点まで）

内容は：

x座標が全て異なるn+1点を通るn次の多項式をラグランジュ補間で求めます。
入力する点は、nが11より小さくて不要な場合でも0を入れておいてください。

というもので計算画面はこれ：

せっかくなんで、ルンゲの現象を起こすy=1/(1+25*x^2）を

例題としていれておいた。

「積層セラミックコンデンサ(MLCC)の加速試験による寿命予測」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP!

久しぶりに最先端のMLCCの話を聞く機会があった。それで何か自分でも計算してみようということで、積層セラミックコンデンサ各社が出している寿命予測の経験式をやってみた。

有名どころ三社のリンク：

https://product.tdk.com/info/ja/contact/faq/faq_detail_D/1432655794672.html

https://www.murata.com/ja-jp/support/faqs/products/capacitor/ceramiccapacitor/qlty/0010

https://www.yuden.co.jp/jp/product/support/faq/q020.html

これらを参考に作ったのがこれ：

 積層セラミックコンデンサ(MLCC)の加速試験による寿命予測

積層セラミックコンデンサ（MLCC)の市場での寿命予測に使われる経験式です。
Ln：通常の使用条件での寿命（時間）
La：加速条件での寿命（時間）
Vn：通常の使用電圧（V）
Va：加速条件での電圧（V）
Tn：通常の使用温度（℃）
Ta：加速条件での温度（℃）
n：電圧加速係数
θ: 温度加速係数