長岡係数¶

長岡係数$C$は次式で与えられる． \begin{gather} C = \frac{4}{3\pi} \frac{1}{k'} \left\{ \frac{k'^2}{k^2} (K-E) + E-k \right\} \end{gather} ただし， \begin{gather} k^2 = \frac{4a^2}{4a^2 + l^2}\\ k' = \sqrt{1-k^2} \end{gather} ここで，$K$は$k$を母数とする第一種完全だ円積分， $E$は$k$を母数とする第二種完全だ円積分を示す． \begin{gather} K = \int_0^{\frac{\pi}{2}} \frac{d\varphi}{\sqrt{1-k^2 \sin ^2 \varphi}}\\ E = \int_0^{\frac{\pi}{2}} \sqrt{1-k^2 \sin ^2 \varphi} \ d\varphi \end{gather}

In [1]:
import numpy as np
import scipy as sp
from scipy import integrate # 数値積分
from scipy.special import ellipk,ellipe # 楕円積分の関数
import matplotlib.pyplot as plt  
import scienceplots # 論文向けプロットスタイル
plt.style.use(['science', 'notebook']) # プロットスタイルの設定
np.set_printoptions(precision=4) # NumPyの出力表示桁数の設定
In [2]:
#plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定
In [3]:
def table(*args): # 表形式で数値を出力する関数
    n1, n2 = np.shape(args)
    for i in range(n2):
        [print(f'{args[j][i]:8.3f}',end=' ') for j in range(n1)]
        print("")
    return
In [4]:
# 長岡係数の計算に用いるパラメータ al = 2a/l（細長さの指標）
al = 0.0001#= 2*a/l
al = np.append(al,np.linspace(0.05, 1, 20)) # 細いコイル領域
al = np.append(al,np.linspace(1.1, 2, 10)) # 中程度のコイル
al = np.append(al,np.linspace(2.5, 5, 6)) # 太いコイル
al = np.append(al,np.linspace(6, 10, 5)) # 非常に太いコイル
#al
In [5]:
# 楕円積分で使用される変数
m = al**2/(al**2+1)
k = np.sqrt(m)
kd = np.sqrt(1-m)
K = ellipk(m)# 第一種完全楕円積分
E = ellipe(m)# 第二種完全楕円積分
In [6]:
C = 4/3/np.pi/kd * (kd**2/k**2*(K-E) + E - k) # 完全楕円積分を用いた長岡係数の計算
In [7]:
fig = plt.figure() # グラフ領域の作成
plt.grid(color = "gray", linestyle="--") # グリッド表示
plt.plot(al,C) # 長岡係数のプロット
plt.xlim(0.0, 10.0)
plt.ylim(0.0, 1.0)
plt.xlabel(r"$2a/l$") # x軸のラベル設定
plt.ylabel("Nagaoka coefficient $C$") # y軸のラベル設定
#plt.legend(ncol=1, loc='lower left', fancybox=False, frameon = True, fontsize=12)
fig.savefig('p3_emt_nagaoka_c.pdf')
No description has been provided for this image

適応積分：ガウスの数値積分 (Gaussian quadrature)¶

In [8]:
def func1(phi,k):
    return 1/np.sqrt(1.0-k**2*np.sin(phi)**2)
def func2(phi,k):
    return np.sqrt(1.0-k**2*np.sin(phi)**2)
In [9]:
# 数値積分による楕円積分の計算（scipy.integrate.quad使用）
n_al = len(al)
K_ = np.empty(n_al)
E_ = np.empty(n_al)
for i in range(n_al):
    K_[i], _ = integrate.quad(func1, 0, np.pi/2, args=(k[i],))
    E_[i], _ = integrate.quad(func2, 0, np.pi/2, args=(k[i],))
In [10]:
integrate.quad(func1, 0, np.pi/2, args=(k[1],)),k[1]
Out[10]:
((1.5717770023396094, 1.7450230175737428e-14), 0.04993761694389224)
In [11]:
print(f"{'2a/l':>8s} {'Nagaoka':>8s} {'k':>8s} {'ellipk':>8s} {'Quad: K':>8s} {'ellipe':>8s} {'Quad: E':>8s}")
m=10 # 表の間引き表示間隔
table(al[::m],C[::m],k[::m],K[::m],K_[::m],E[::m],E_[::m])

    2a/l  Nagaoka        k   ellipk  Quad: K   ellipe  Quad: E
   0.000    1.000    0.000    1.571    1.571    1.571    1.571 
   0.500    0.818    0.447    1.660    1.660    1.489    1.489 
   1.000    0.688    0.707    1.854    1.854    1.351    1.351 
   2.000    0.526    0.894    2.257    2.257    1.178    1.178 
   9.000    0.219    0.994    3.598    3.598    1.019    1.019