【Python】フーコーの振り子の寿命の質量依存性を調べる

概要

　フーコーの振り子の寿命の質量依存性を調べた。そのために並列処理を実装した。結果は質量が20kg以上の時、質量依存性がほとんどなくなった。

動作環境

• Windows10(64bit)
• Python 3.7.7

やりたいこと

　フーコーの振り子の寿命の質量依存性を調べたい。そのためにフーコーの振り子の振幅が10000回連続して0.1m以下になったときの時刻を調べる。
また、並列処理を実装した。

アルゴリズム

1. ルンゲクッタ法である時刻の振り子の位置を求める。
2. 振り子の振幅を計算する。
3. 振り子の振幅が0.1m以下だったらcountを1増やす。振り子の振幅が0.1m以上だったらcountを初期化する。
4. countが10000に達したときの時刻を取得する

以上のアルゴリズムを並列処理する。

プログラム

import numpy as np
from mpmath import pi
from mpmath import exp
from mpmath import sin
from mpmath import sqrt
import numpy.linalg as LA
import matplotlib.pyplot as plt
from multiprocessing import Pool
import mpmath
mpmath.mp.dps=100 #計算精度100桁

#関数を定義する

#微分方程式その1
def f1(t,x0,x1,y0,y1):
    return x1

#微分方程式その2
def f2(t,x0,x1,y0,y1,i_m):
    g = 9.8
    th = 40.5
    l = 45.0
    o = 7292115.0*(10**(-11))
    k = 0.0037
    m = i_m
    a = -g/l
    b = 2.0*o*sin(th*pi/180)
    c = -k/m
    return a*x0-b*y1+c*(x1)

#微分方程式その3
def f3(t,x0,x1,y0,y1):
    return y1

#微分方程式その4
def f4(t,x0,x1,y0,y1,i_m):
    g = 9.8
    th = 40.5
    l = 45.0
    o = 7292115.0*(10**(-11))
    k = 0.0037
    m = i_m
    a = -g/l
    b = 2.0*o*sin(th*pi/180)
    c = -k/m
    return a*y0+b*x1+c*(y1)

#ルンゲクッタ法
def g(i_m):
    x0 = 0.0　#初期位置
    x1 = 0.0　#初期位置
    y0 = 1.0　#初期位置
    y1 = 0.0　#初期位置
    t_initial = 0.0　#初期時刻
    t_final = 5.0*60.0*60.0　#最終時刻
    N = 1800000
    h = (t_final-t_initial)/N #刻み幅
    array_t = np.arange(t_initial,t_final,h)
    array_amplitude = np.empty((array_t.size),dtype=float)
    count = 0
    for i in range(N):
        k_11 = h*f1(array_t[i],x0,x1,y0,y1)
        k_12 = h*f2(array_t[i],x0,x1,y0,y1,i_m)
        k_13 = h*f3(array_t[i],x0,x1,y0,y1)
        k_14 = h*f4(array_t[i],x0,x1,y0,y1,i_m)
        k_21 = h*f1(array_t[i]+(h/2),x0+((k_11)/2),x1+((k_12)/2),y0+((k_13)/2),y1+((k_14)/2))
        k_22 = h*f2(array_t[i]+(h/2),x0+((k_11)/2),x1+((k_12)/2),y0+((k_13)/2),y1+((k_14)/2),i_m)
        k_23 = h*f3(array_t[i]+(h/2),x0+((k_11)/2),x1+((k_12)/2),y0+((k_13)/2),y1+((k_14)/2))
        k_24 = h*f4(array_t[i]+(h/2),x0+((k_11)/2),x1+((k_12)/2),y0+((k_13)/2),y1+((k_14)/2),i_m)
        k_31 = h*f1(array_t[i]+(h/2),x0+((k_21)/2),x1+((k_22)/2),y0+((k_23)/2),y1+((k_24)/2))
        k_32 = h*f2(array_t[i]+(h/2),x0+((k_21)/2),x1+((k_22)/2),y0+((k_23)/2),y1+((k_24)/2),i_m)
        k_33 = h*f3(array_t[i]+(h/2),x0+((k_21)/2),x1+((k_22)/2),y0+((k_23)/2),y1+((k_24)/2))
        k_34 = h*f4(array_t[i]+(h/2),x0+((k_21)/2),x1+((k_22)/2),y0+((k_23)/2),y1+((k_24)/2),i_m)
        k_41 = h*f1(array_t[i]+h,x0+(k_31),x1+(k_32),y0+(k_33),y1+(k_34))
        k_42 = h*f2(array_t[i]+h,x0+(k_31),x1+(k_32),y0+(k_33),y1+(k_34),i_m)
        k_43 = h*f3(array_t[i]+h,x0+(k_31),x1+(k_32),y0+(k_33),y1+(k_34))
        k_44 = h*f4(array_t[i]+h,x0+(k_31),x1+(k_32),y0+(k_33),y1+(k_34),i_m)
        x0 = x0+(1/6)*((k_11)+2*(k_21)+2*(k_31)+(k_41))
        x1 = x1+(1/6)*((k_12)+2*(k_22)+2*(k_32)+(k_42))
        y0 = y0+(1/6)*((k_13)+2*(k_23)+2*(k_33)+(k_43))
        y1 = y1+(1/6)*((k_14)+2*(k_24)+2*(k_34)+(k_44))
        array_f = np.array([x0,y0])
        array_amplitude[i] = LA.norm(array_f)
    for i_t in range(array_t.size-2):
        if array_amplitude[i_t] > 0.1:
            count = 0
        elif array_amplitude[i_t] <= 0.1:
            count += 1
            if count == 10000:
                break
    print(i_m)
    print(array_t[i_t-10000])
    return array_t[i_t-10000]

#寿命の質量依存性のグラフを書く
def task():
    plt.figure()
    plt.plot(array_m,array_lifespan.T,'blue',label='Calc')
    plt.title('Relation of mass and lifespan')
    plt.xlabel('mass')
    plt.ylabel('lifespan')
    plt.legend(loc=0)
    plt.savefig("Relation of mass and lifespan_1.png")


    
#各mにおける寿命を計算して、結果を配列にする
m_first = 1
m_last = 100
array_m = np.arange(m_first,m_last,1)

#並列処理
if __name__=="__main__":
    pool = Pool(3)
    array_lifespan = np.array(pool.map(g,array_m))
    task()

結果

　グラフからフーコーの振り子の寿命には質量依存性があることがわかった。また質量依存性は20kgからほとんどなくなる。

参考文献