python 応用ルンゲクッタ法

Question

import numpy as np
import matplotlib.pyplot as plt
from numpy.lib.type_check import isrealobj

#関数の定義
def f1(fai2):
    return fai2

def f2(fai1,fai2):
    return (i-np.sin(fai1)-fai2)/bc

#初期値
a = 0
b = 100
N = 20000
h = (b-a)/N

bc = 4

#配列
fai1 = np.empty(N)
fai2 = np.empty(N)
fai1_0 = 0
fai2_0 = 0
fai1[0] = fai1_0
fai2[0] = fai2_0

for i in np.linspace(0, 2, 201):

    for j in range (0, N-2, 1):
       k1 = h * f1(fai2[j])
       d1 = h * f2(fai1[j] , fai2[j])
       k2 = h * f1(fai2[j] + d1/2)
       d2 = h * f2(fai1[j] + k1/2, fai2[j] + d1/2)
       k3 = h * f1(fai2[j] + d2/2)
       d3 = h * f2(fai1[j] + k2/2, fai2[j] +d2/2)
       k4 = h * f1(fai2[j] + d3)
       d4 = h * f2(fai1[j] + k3 , fai2[j] + d3)     
       fai1[j+1] = fai1[j] + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) 
       fai2[j+1] = fai2[j] + 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) 

       if j == N-1:
           break

    #6000以降のfai2の合計の平均
    v = np.average(fai2[10000:])
    print(v)
    print(i)
    plt.scatter(v, i, marker='.',color = "green")  

fai2_0 = fai2[19999]
fai1_0 = fai1[19999]


for i in np.linspace(2, 0, 201):

    for j in range (N-2, 0, -1):
       k1 = h * f1(fai2[j])
       d1 = h * f2(fai1[j] , fai2[j])
       k2 = h * f1(fai2[j] + d1/2)
       d2 = h * f2(fai1[j] + k1/2, fai2[j] + d1/2)
       k3 = h * f1(fai2[j] + d2/2)
       d3 = h * f2(fai1[j] + k2/2, fai2[j] +d2/2)
       k4 = h * f1(fai2[j] + d3)
       d4 = h * f2(fai1[j] + k3 , fai2[j] + d3)     
       fai1[j-1] = fai1[j] + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) 
       fai2[j-1] = fai2[j] + 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) 

       if j == 0:
           break

    #6000以降のfai2の合計の平均
    v = np.average(fai2[:10000])
    
    print(v)
    print(i)
    plt.scatter(v, i, marker='.',color = 'blue') 
    
#plot
plt.plot(v,i)   
plt.xlim(0, 1.8)
plt.ylim(0, 2.0)
plt.xlabel("v")
plt.ylabel("i")
plt.show()
```　このコードは、規格化されたRCSJモデル
　　　　　i = sinφ　+ dφ/dτ + βc*d^2φ/dτ^2
という電流iと電圧vのグラフをプロットするために書いたコードです。（ここではφをfai1,dφ/dτをfai2とおいています）
　電流を0から2まで増やしていったあと、そのまま電流を2から0まで減らしていくと以下のグラフ1のように帰ってくる値が０ではなく、与えたβcの値によって変わるヒステリシスという現象が起こります。このコードをルンゲクッタ法を用いて書いていきます。
　本来であれば、電流iは規格化されているので臨界電流としてi＝1.0となるまでは電圧vはグラフのようにほぼ0の値をとり続けるのですが、実際にコードを回してみると、なぜかi＝1.0に行く前に電圧が発生しているグラフがプロットされます。結果の数値を表示すると、i=0.94でvの値が大きくなってしまっていました（グラフ２）
それ以前の問題だと思い、φ2とτのグラフを確認してみたところ、それらのグラフはきちんとプロットされていたので、vとiのグラフをプロットした途端に正しく出ないのはおかしいと思いました。
　どこかコードがおかしいのでしょうか。それとも、小数点以下の誤差なので、どこかループでの誤差の影響なのでしょうか。アドバイスを頂けたらと思います。

グラフ1![イメージ説明](7f5eb9303100c07a76838b9c3867e472.png)

グラフ2![イメージ説明](b7d983fd26d34267deba8efc111ecb7f.png)

Comments

[meg_]

> グラフプロットの際の誤差
計算は合っているのにグラフ化するとおかしくなる、という現象が起きているのでしょうか？

[Ackngawe.-]

申し訳ありません。計算した際の値も表示しているのですが、そこで値がおかしくなっていました。
すぐにタイトルを修正します。

[退会済みユーザー]

あるべき姿（？）を手書きでもよいので画像に加えていただくことはできますか？

[Ackngawe.-]

本来であれば、緑色の点（電流iを増加させたときの点）がi=1.0に到達してから電圧vが発生するので、あるべき姿はグラフ1の赤色の点線のようになるはずです。

[melian]

規格化されたRCSJモデルに関する資料(ジョセフソン接合関連)を２本程斜め読みしたのですが、f2 関数内の計算式にある np.sin(fai1) は np.cos(fai1) じゃないかな、と。

def f2(fai1,fai2):
##return (i-np.sin(fai1)-fai2)/bc
return (i-np.cos(fai1)-fai2)/bc

np.cos に変更して実行すると i=0 に到達してから電圧が発生します。

[Ackngawe.-]

確かにcosにすると正しくプロットできていました、、
ですが、ジョセフソン効果で、I＝Ic*sinθの関係からそれを用いて規格化しているので、そこをcosに変えてしまうとそもそもがおかしくなってしまわないでしょうか、、？

[melian]

v の平均値に影響はなく、位相が π/2 ずれるだけかと思います。ズレがどのくらいかと言うと、bc * i_tick * (cos(0) - sin(0)) です。i_tick は for i in np.linspace(0, 2, 201) としているので 0.01 です。つまり、 4 * 0.01 * 1 = 0.04 ずれるので 0.96 で電流が発生してしまう事になります（実際には 0.94 ですが）。

[Ackngawe.-]

規格化された通りの式でコードを入力しているつもりなのですが、なぜそのようなズレが生じてしまっているのでしょうか？
それとも、コードのどこかでそのズレを生じさせてしまっているということでしょうか？

Answer 1

import numpy as np
import matplotlib.pyplot as plt

def f1(tau,fai2:float):
    return fai2

def f2(tau,fai1,fai2:float):
    return (i-np.sin(fai1)-fai2)/BC

i = 1.0
BC = 4

a = 0
b = 100
N = 20000
h = (b-a)/N
fai1 = 0.0
fai2 = 0.0
s = float()
v = float()

ln, = plt.plot([], [])

for i in np.arange(0, 2, 0.01):

    for j in range(0, N, 1):
        tau = j * h
        k1 = h * f1(tau,fai2)
        d1 = h * f2(tau,fai1 , fai2)
        k2 = h * f1(tau,fai2 + d1/2)
        d2 = h * f2(tau,fai1 + k1/2,  fai2 + d1/2)
        k3 = h * f1(tau,fai2 + d2/2)
        d3 = h * f2(tau,fai1 + k2/2 , fai2 +d2/2)
        k4 = h * f1(tau,fai2 + d3)
        d4 = h * f2(tau,fai1 + k3 , fai2 + d3 )     
        fai1 += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) 
        fai2 += 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) 
        
        if j > 10000:
           s += fai2
    
    v = s / 10000
    print(v)
    print(i)
    ln.set_xdata(np.append(ln.get_xdata(), v))
    ln.set_ydata(np.append(ln.get_ydata(), i))
    plt.scatter(v, i, marker='.',color = 'blue') 
    s = 0

for i in np.arange(2, 0, -0.01):

    for j in range(0, N, 1):
        tau = j * h
        k1 = h * f1(tau,fai2)
        d1 = h * f2(tau,fai1 , fai2)
        k2 = h * f1(tau,fai2 + d1/2)
        d2 = h * f2(tau,fai1 + k1/2,  fai2 + d1/2)
        k3 = h * f1(tau,fai2 + d2/2)
        d3 = h * f2(tau,fai1 + k2/2 , fai2 +d2/2)
        k4 = h * f1(tau,fai2 + d3)
        d4 = h * f2(tau,fai1 + k3 , fai2 + d3 )     
        fai1 += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) 
        fai2 += 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) 
        
        if j > 10000:
           s += fai2
    
    v = s / 10000
    print(v)
    print(i)
    ln.set_xdata(np.append(ln.get_xdata(), v))
    ln.set_ydata(np.append(ln.get_ydata(), i))
    plt.scatter(v, i, marker='.',color = 'green') 
    s = 0

y = ln.get_ydata()     
plt.plot(v,i)
plt.xlim(0, 1.8)
plt.ylim(0, 2.0)
plt.xlabel("v")
plt.ylabel("i")
plt.show()
```配列を使わないコードでなんとかできました。

Comments

[退会済みユーザー]

解決できて何よりです。手元で動かしましたが計算が遅かったのを覚えています。
参考までに高速化方法を書きます。

from numba import jit
として、
@jit()
def f2(fai1,fai2):
とするとかなり早くなります。