たのしい!数値計算
この記事は、奈良高専 Advent Calendar 2022の24日目の記事です。前回の方(Jin)のブログはこちらです。
はじめに
皆さん初めまして。ぴかです。
最近BLEACHの千年血戦篇を見ています。めちゃくちゃおもろいですね。
さて、今回は数値計算をテーマになんか書いていこうと思います。なんで数値計算かというと、クリスマスだからです*1。クリスマスは数値計算をしろ。
わたくしの分際でクリスマスイブの記事を書くのは、たいへんおこがましいと思いつつも、なんかごちゃごちゃと書いていこうかなといった次第であります。
目次?
たいそうな文章ではないため、目次とは言えない何かです。うーん、おしながきといったところでしょうか。
数値計算ってなに~
数値計算ってなんなんでしょうね。数値を使って計算するなら算数じゃね?って思われるかもしれません。デジタル大辞泉によると、
コンピューターを使って、代数的または解析的に解くことが困難な問題、または手計算では不可能な計算量が膨大な問題を解くこと。
とのことです。
すなわち、手計算するのがしんどかったり、解析的(大雑把に言うと式変形)に解くのが難しいものを計算する方法が数値計算法というわけですね。
数値計算の対象
数値計算だったり数値解析の対象となるものはめっちゃ多いです。
数値補間、行列式、微積分、微分方程式(常微分も偏微分も)とかその他たくさんあります。
その応用例は、流体解析とか電磁界解析、ウイルスの感染*2をシミュレーションしたりとかです。身近なものだと天気予報もあります。
実際に数値計算してみる
実際に数値計算してみます。むずかしいモデルを数値計算するには、その分プログラム書くのも大変だったりします。
そのため、ここでは古典力学上の簡単な計算例を紹介します。
自由落下問題
質量をもった物体(ここでは質点とします)がある高さから落下していくモデルを例に数値計算します。
まず、自由落下について理論式から考えていきます。質量を$ m $、速度を$ v(t) $、重力加速度を$ g $とすれば、Newtonの運動方程式は、
$$ m\frac{dv(t)}{dt} = mg \tag{1} $$
となります。
この微分方程式を解けば、物体が移動する距離や物体の速度が求まるわけです。
ところで微分ってどうやって数値計算するんでしょうか??
微分はTaylor展開により表すことができます。
$ t = \Delta t $まわりのTaylor展開により$ \frac{dv(t)}{dt} $は、
$$ \frac{dv(t)}{dt} \approx \frac{v(t + \Delta t )-v(t)}{\Delta t} \tag{2} $$
となります。(これは、$ \Delta t \approx 0 $とすれば微分の定義$ \frac{dv(t)}{dt} = \lim_{\Delta t \rightarrow 0} \frac{v(t + \Delta t )-v(t)}{\Delta t} $と一致し、微分を近似表現できていると言えます。)
ここで、速度の微分を$ \frac{dv(t)}{dt} = f(t,v) $とすると、式(2)は、
$$ v(t+\Delta t) = v(t) + f(t,v) \cdot \Delta t \tag{3} $$
こんな感じになります。同様に、移動距離$ x(t) $についても、
\begin{align} x(t+\Delta t) &= x(t) + \frac{dx(t)}{dt} \cdot \Delta t \\ &= x(t) + v(t) \cdot \Delta t \tag{4} \end{align}
と表すことができます。
これら式(3)、(4)により微分方程式の数値解を計算する方法をEuler法といいます。
このEuler法によって式(1)で表される自由落下を数値計算します。なおPythonで数値計算しました。Pythonは個人的に使いやすいなっていうのと、colaboratoryとかweb上で実行できたりするので便利です(VScodeのほうが使いやすいけど)。
コードはこんな感じ(きたないかもだけど許して。あとスマホだとこのコード表示が見にくくなるんやが、どうしたらいいんでしょうか。)。
#簡単な力学系のシミュ(1次元質点運動問題) #ここでは自由落下を扱う(空気抵抗を一切含んでいない F = mg の式) import numpy as np import matplotlib.pyplot as plt #const g = 9.8 #計算エンジン------------------------------------------------------------ #変数 t = 0.0 #時刻t d = 0.01 #微分刻み幅d #格納配列(csv出力とかグラフとかにするため) tlist = [] ylist = [] vlist = [] #入力部 y = float(input("初期高度y0を入力してほしいにゃ:")) #初期位置を入力 v = float(input("初速度v0を入力してほしいにゃ:")) #初速を入力 print("{:.5f} {:.5f} {:.5f}" .format(t,y,v)) #入力時の数値出力 #自由落下計算 while y >= 0: #高度が0になるまでの時間 t += d #時刻更新 v += g*d #速度更新 y -= v*d #位置更新 print("{:.5f} {:.5f} {:.5f}" .format(t,y,v)) #計算課程から終了まで数値出力 tlist.append(t) ylist.append(y) vlist.append(v) #グラフ出力します figy = plt.figure(1) plt.plot(tlist,ylist) plt.xlabel('Time [s]') plt.ylabel('Height [m]') plt.show() figv = plt.figure(2) plt.plot(tlist, vlist) plt.xlabel('Time [s]') plt.ylabel('Velocity [m/s]') plt.show()
初期位置を50 m、初速を0 m/sとして実行します。なお計算条件は、時刻の刻み幅を0.01 s、質点の高度が0になる時刻で計算終了としました。
質点の高度変化をFig3に、速度変化をFig4に示します。
わ~、自由落下運動の数値解を求めれた~......
実はこれ、空気抵抗を考慮していないモデルなので、初期高度が高ければ高いほど落下時における質点の速度は速くなります。
例えば雨を降らす乱層雲は地上から1 km~2 kmぐらいの高さですが、このとき空気抵抗が無いと雨粒は140 ~198 m/sの速度で降ることになります(雨が10 kgの量降るとしたら)。
時速換算すると、504~713 km/h程度です(この速度で雨が降ればヤバイですよね)。
ということで空気抵抗を考慮した自由落下を考えましょう。
閑話休題
空気抵抗を考慮する前にちょっとだけ。
最近Polyphiaっていうバンドの曲をよく聞いてるんですが、まぁえぐいんですね。インストバンドなんですが、ギターがうますぎる。ギターはTim HensonとScott Lepageの2人なんですが、2人ともえげつうまい。
Playing Godとかえげついです。クラシックギターでピロピロ弾いてるし。ぜひ聞いてみてね。
空気抵抗を考慮した自由落下
物理の問題とかで速度に比例する空気抵抗とかよくありますよね*3。アレを考えます。
速度$ v $に比例する空気抵抗を考慮した場合、運動方程式は、
$$ m\frac{dv(t)}{dt} = mg - kv(t) \tag{5} $$
と表すことができます。ここで$ k $は空気抵抗の比例係数です。
この比例定数は粘性抵抗と言われるぽいです。
式(5)について初速を0とし、速度$ v(t) $を導くと、
$$ v(t) = \frac{mg}{k}\Bigr ( 1 - \exp\bigr [ -\frac{kt}{m} \bigl ] \Bigl ) \tag{6} $$
となります。
時間$ t $が十分長いとき、速度は、
$$ \lim_{\Delta t \rightarrow \infty} v(t) = \frac{mg}{k} \tag{7} $$
となります。
この速度を終端速度といい、質点は一定の速度に達します。
さて、この空気抵抗を考慮した自由落下を数値計算しましょう。
コードはこんな感じ。
#簡単な力学系のシミュ(1次元質点運動問題) #速度に比例する空気抵抗(粘性抵抗)がかかる自由落下を扱う(F = mg - kv の式) import numpy as np import matplotlib.pyplot as plt #const g = 9.8 #計算エンジン------------------------------------------------------------ #変数 t = 0.0 #時刻t d = 0.01 #微分刻み幅d #格納配列(csv出力とかグラフとかにするため) tlist = [] ylist = [] vlist = [] #入力部 m = float(input("質量mを入力してほしいにゃ:")) #物体の質量を入力 k = float(input("空気抵抗kを入力してほしいにゃ:")) #空気抵抗を入力 y = float(input("初期高度y0を入力してほしいにゃ:")) #初期位置を入力 v = float(input("初速度v0を入力してほしいにゃ:")) #初速を入力 print("{:.5f} {:.5f} {:.5f}" .format(t,y,v)) #入力時の数値出力 #自由落下計算 while y >= 0: #高度が0になるまでの時間 t += d #時刻更新 v += (g-k/m*v)*d #速度更新 y -= v*d #位置更新 print("{:.5f} {:.5f} {:.5f}" .format(t,y,v)) #計算課程から終了まで数値出力 tlist.append(t) ylist.append(y) vlist.append(v) #グラフ出力します figy = plt.figure(1) plt.plot(tlist,ylist) plt.xlabel('Time [s]') plt.ylabel('Height [m]') plt.show() figv = plt.figure(2) plt.plot(tlist, vlist) plt.xlabel('Time [s]') plt.ylabel('Velocity [m/s]') plt.show()
質量を10 kg、比例定数を10 kg/s、初期位置を50 m、初速を0 m/s、として実行します。なお計算条件は、時刻の刻み幅を0.01 s、質点の高度が0になる時刻で計算終了としました。
質点の高度変化をFig7に、速度変化をFig8に示します。
ちゃんと速度が一定のとこ(およそ9.8 m/s)に漸近していますね。
質点の高度変化とか速度変化をグラフだけで見ていても退屈になってきたので、アニメーションで質点の動きを見てみましょう。
動き始めよりも地面に近づくときのほうが速い落下になってますね。
おわりに
自由落下という簡単な物理現象を数値計算してシミュレーションしました。わりとおもろいですよね。
微分方程式で表される現象は古典力学の運動だけでなく、電磁界を表すMaxwell方程式や流体のNavier-Stokes方程式であったりと、たくさんあります。(私は電磁界解析を研究でやってたりやってなかったり)
気になった方はぜひ数値計算してみてください。
おまけ
アニメーションを含むコードです。
# 簡単な力学系のシミュ(1次元質点運動問題) # 速度に比例する空気抵抗(粘性抵抗)がかかる自由落下を扱う(F = mg - kv の式) import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation # const g = 9.8 # 計算エンジン------------------------------------------------------------ # 変数 d = 0.01 # 微分刻み幅d t = 0.0 # 時刻 # 格納配列(csv出力とかグラフとかにするため) tlist = [] ylist = [] vlist = [] td = np.arange(0, 100, 0.01) # アニメ時刻 # 入力部 m = float(input("質量mを入力してほしいにゃ:")) # 物体の質量を入力 k = float(input("空気抵抗kを入力してほしいにゃ:")) # 空気抵抗を入力 y = float(input("初期高度y0を入力してほしいにゃ:")) # 初期位置を入力 v = float(input("初速度v0を入力してほしいにゃ:")) # 初速を入力 print("{:.5f} {:.5f} {:.5f}" .format(t, y, v)) # 入力時の数値出力 y0 = y # 自由落下計算 while y >= 0: # 高度が0になるまでの時間 t += d # 時刻更新 v += (g-k/m*v)*d # 速度更新 y -= v*d # 位置更新 print("{:.5f} {:.5f} {:.5f}" .format(t, y, v)) # 計算課程から終了まで数値出力 tlist.append(t) ylist.append(y) vlist.append(v) t_posi0 = t # 地面に到着する時間 # 計算 ind_0 = np.argmin(np.abs(td-t_posi0)) ylist = ylist[:ind_0] # データをグラフに(ここでは高度変化) data = np.concatenate([ylist]) plt.plot(data) plt.ylabel("Height [m]") plt.xlabel("Time [s]") plt.show() data.shape # アニメーション fig, ax = plt.subplots(figsize=(6, 4)) ball, = ax.plot(0, data[0]+5, 'C2o', ms=15) ax.hlines(0, -0.25, 0.25) ax.set_xticks([]) ax.set_ylim(-10, data.max()+20) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.spines['bottom'].set_visible(False) ax.set_title("$Height$="+str(y0)+"[m]") def update(num): ball.set_ydata(data[num]+5) ani = animation.FuncAnimation(fig, update, 1000, interval=10) ani.save('ball.gif', dpi=100)
一昨日、全然寝れてなくて執筆・投稿が少し遅くなりました。ごめんなさい🎅
ちなみに「たのしい!数値計算」ってタイトルですが、アニメーションとか見てるときはたのしいけど、それ以外なんもたのしくなかったです、ねむすぎて。
n岡技科大 H31 電電
n岡技科大の電電を解いた。回路,電磁気関連しか解いていない。(電子回路とかは元気があればやるかも)
問題1 電気機器に関する問題
問1
極数を$ p $、周波数を$ f $とすると同期速度は、
$$ N_0 = \frac{120f}{p} $$
で表される。したがって、
$$ N_0 = \frac{120f}{p} = \frac{120\times 50}{4} = 1500 \rm{[r/m]} $$
問2
すべりは、
$$ s = \frac{N_0 - N}{N_0} = \frac{1500-1440}{1500}=0.04 $$
問3
(無負荷試験) 定格電圧$ V_0 = 10 \sqrt{3} $ [V]を加えて無負荷試験したところ、入力電流$ I_0 = 0.5 $ [A]、入力電力$ P_0 =12 $ [W]であった。
(拘束試験) 回転子を拘束して入力電圧$ V_s = 3\sqrt{3} $ [V]、50 Hzを加えたところ、入力電流$ I_s = \sqrt{\frac{3}{2}} $ [A]、入力電力$ P_s = 9 $ [W]であった。なお、固定子1相分の巻線抵抗$ r_1 $は0.5 $ \rm{\Omega} $である。
(1)励磁回路の損失分抵抗$ g_0 $ [S]
(2)励磁リアクタンス$ b_0 $ [S]
(3)固定子から見た回転子巻線の抵抗$ r_2^{'} $ [$ \rm{\Omega} $]
(4)固定子巻線と回転子巻線の漏れリアクタンスの和$ X $ [$ \rm{\Omega} $]
を求めよ。
(1) 無負荷試験のとき、すべりが無い。したがって$ r_2^{'}\frac{1-s}{s} =\infty $となり、励磁回路のみ残る。等価回路はこんな感じ。
簡単な回路にできたので、無負荷試験の条件をそのまま解くと$ g_0 $が求まる。
$$ P_0 = V_0I_0 = g_0 {V_0}^{2} = 300g_0 = 12 $$
$$ g_0 = \frac{12}{300} = 0.04 \rm{[S]} $$
(2) (1)の結果と無負荷試験の条件より$ b_0 $が求まる。
アドミタンスを$ Y $とすると、
$$ Y = \frac{I_0}{V_0} = \frac{\sqrt{3}}{60} = \sqrt{g_0^2 + b_0^2} $$
$$ b_0^2 = \Bigl( \frac{\sqrt{3}}{60} \Bigr) ^2- {0.04}^2 = \frac{-23}{30000} $$
$$ b_0 = \sqrt{\frac{23}{30000}} \rm{[S]} $$
(3)拘束試験のとき、すべりが1となる。したがって$ r_2^{'}\frac{1-s}{s} =0 $となる。また、電流のほとんどが負荷に流れるため励磁回路は無視できる。等価回路はこんな感じ。
拘束試験の条件より簡単に$ r_2^{'} $が求まる。
$$ P_s = V_s I_s = (0.5+r_2^{'})I_s^2 = 9$$
$$ r_2^{'} = \frac{9}{I_s^2} - 0.5 = 5.5 \rm{[\Omega]} $$
(4) (3)の結果と拘束試験の条件より$ X $は求まる。
インピーダンスを$ Z $とすると、
$$ Z = \frac{V_s}{I_s} = 3\sqrt{2} = \sqrt{(0.5+r_2^{'})^2 +X^2} $$
$$ X^2 = (3\sqrt{2})^2 - (0.5+5.5)^2 = -18 $$
$$ X = 3\sqrt{2} \rm{[\Omega]} $$
問題2 電気回路、伝達関数に関する問題
問1
$ R_2 $が十分大きい、つまり開放していると見なせる。すなわち$ i_R = 0 $となる。このときインピーダンスは、
$$ Z_{all} = R_1 + \frac{j\omega L (-j\frac{1}{\omega C})}{j\omega L - j\frac{1}{\omega C}} = R_1 + j \frac{\omega L}{1- \omega ^2 LC} \rm{[\Omega]} $$
問2
角周波数$ \omega $で微分して、0になる点を知ればいい。
$$ \frac{d}{d\omega} \{Z_{all}\} = j \frac{L(\omega ^2 LC -1)}{(1-\omega ^2 LC)^2} =0 $$
つまり、
$$ \omega ^2 LC -1 = 0 $$
$$ \omega = \frac{1}{\sqrt{LC}} \rm{[rad/s]} $$*1
問3
Fig.2-1の回路方程式は、
$$ e_{in} = R_1 i_1 + \frac{1}{C} \int{i_C} dt $$
である。これをラプラス変換すると、
$$ \frac{E_{in}}{s} = R_1 I_1 + \frac{1}{Cs} I_c $$
となる。
こんな感じでやっていけば、
$$ (1) \Rightarrow sR_1 $$
$$ (2) \Rightarrow s^2 L $$
$$ (3) \Rightarrow \frac{1}{C} $$
$$ (4) \Rightarrow sR_2 $$
となる。
問4
伝達関数は、
$$ G(s) = \frac{E_{out}(s)}{E_{in}(s)} = \frac{I_C}{sCR_1 I_1 + I_C} = \frac{sLR_{2}}{R_{1} R_{2} +s(LR_{1}+L R_{2}) +s^{2} LCR_{1}R_{2}} $$
となる。*2
問5
問題3 電磁誘導に関する問題
(1)~(9)に入る数式、語句を埋めよ。
円運動するときの半径$ r $と角速度の関係より、
$$ (1) \Rightarrow v= \omega r $$
となる。
電荷量を$ q $とするとLorentz力は、$ \boldsymbol{F}=q(\boldsymbol{v} \times \boldsymbol{B}) $であるから、
$ \boldsymbol{E} =\boldsymbol{v} \times \boldsymbol{B} $となる。*3
この電場の大きさから起電力を求めると、
$$ de = | \boldsymbol{E} | dr = vBdr $$
となる。つまり、
$$ (2) \Rightarrow de= \omega r \mu_{0} H dr $$
である。
(2)を原点Oから半径$ a $まで積分すると、
$$ \int_{0}^{a} {de} = e = \frac{\mu_0 H \omega a^2}{2} $$
となる。つまり、
$$ (3) \Rightarrow e= \frac{\mu_0 H \omega a^2}{2} $$
である。
ここで起電力の向きについて考えよう。Neumann-Lenzの電磁誘導則よりうず電流を妨げる向きに起電力は生じる。この円板導体には磁場によって反時計回りにうず電流が生じるが、その逆向きに誘導電流が生じることになる。つまり起電力は起電力は点Oから点Cの向きに生じる。したがって電流はBの向きである。
$$ (4) \Rightarrow \text{Bの向き} $$
Ohmの法則より電流$ I $は、
$$ I = \frac{e}{R} = \frac{\mu_0 H \omega a^2}{2R} $$
となる。つまり、
$$ (5) \Rightarrow I = \frac{\mu_0 H \omega a^2}{2R} $$
である。
半径要素$ dr $中の電流が受ける力は、$ F = IBl $より、
$$ dF = IBdr = \frac{\omega (\mu_0 H a)^2}{2R} dr $$
となる。つまり、
$$ (6) \Rightarrow dF = \frac{\omega (\mu_0 H a)^2}{2R} dr $$
である。
この回転中心軸回りの力のモーメントは、$ d\boldsymbol{T} = \boldsymbol{r}\times d\boldsymbol{F} $だから、
$$ dT = r dF = \frac{r \omega (\mu_0 H a)^2}{2R} dr $$
となる。つまり、
$$ (7) \Rightarrow dT = \frac{r \omega (\mu_0 H a)^2}{2R} dr $$
である。
この微小半径に加わる力のモーメントを原点Oから半径$ a $まで積分すると、
$$ \int_{0}^{a} {dT} = T = \frac{\omega (\mu_0 H )^2 a^4}{4R} $$
となる。つまり、
$$ (8) \Rightarrow T = \frac{\omega (\mu_0 H )^2 a^4}{4R} $$
である。
なおこの力のモーメントは回転方向とは逆向きに働くので、回転を減速させる方向といえる。
$$ (9) \Rightarrow \text{減速} $$
どきどき物理
こんにちは。あったかくなりましたね。
突然ですがみなさんは最近ドキドキしましたか?
恋人と遊んだり、彼女に包丁を向けられたり、狭い道で対向車のトラックがセンターライン超えてたり、物理を勉強したり...
いろんなドキドキがあります。
さて今日は力学の問題を解いていきます。
問題1
真っ先に思いつくのが高さ$ z $を微分して最大値を求めるって方法です。(厳密に言うと極値を求める作業)
$ z $ですが変数が時間$ t $の時のほうが微分しやすそうですね。
$ z(t) $を微分すると、
$$ \frac{dz}{dt} = \frac{d}{dt} \Bigl\{ (v_0 \mathrm{sin}\theta)t-\frac{1}{2}gt^2 \Bigr\} = v_0\mathrm{sin}\theta - gt $$
となります。
これが0のとき、
$$ v_0\mathrm{sin}\theta - gt = 0 $$
$$ t=\frac{v_0}{g}\mathrm{sin}\theta $$
とわかりました。これは最高地点到達時の経過時間そのものです。
この最大値をとる$ t $を$ z(t) $に適用すると、
$$ z=\frac{{v_0}^2}{2g}\mathrm{sin}^2\theta $$
と最高地点の高さが求まりました。
問題2
これも問題1と同様に解いていけばいいです。水平距離$ x $と高さ$ z $の式、すなわち運動の軌跡を見てみましょう。
$$ z=(\mathrm{tan}\theta)x-\frac{g}{2{v_0}^2\mathrm{cos}^2\theta}x^2 $$
そもそもこの運動は放物運動なので投げ始めと投げ終わりは$ z = 0 $になるはず。
投げ始めは$ t=0 $とすればよい。そこで今知りたい到達距離は時間に依存しない形にしたいわけです。($ z = 0 $のときの時間$ t $を$ x $に代入するのでもいいが、まわりくどい)
このことから$ z=0 $とすると、
$$ 0 =(\mathrm{tan}\theta)x-\frac{g}{2{v_0}^2\mathrm{cos}^2\theta}x^2 $$
$$ x = \frac{2{v_0}^2\mathrm{cos}^2\theta}{g} \mathrm{tan}\theta = \frac{2{v_0}^2\mathrm{cos}\theta\mathrm{sin}\theta}{g} $$
と$ x $に関する式がつくれました。これが水平到達距離です。
この水平到達距離$ x $を$ x=(v_0 \mathrm{cos}\theta )t $に代入すると、到達時の時間$ t $を求めることができます。
$$ \frac{2{v_0}^2\mathrm{cos}\theta\mathrm{sin}\theta}{g} = (v_0 \mathrm{cos}\theta )t $$
$$ t = \frac{2v_0}{g} \mathrm{sin}\theta $$
となります。
さらに、$ x(t) $から$ x(\theta) $になったので 、この水平到達距離が最大になるときの角度も求めることができます。
やはりこれも極値そのものなので微分して0になるのを求めます。
$$ \frac{dx}{d\theta}=\frac{2{v_0}^2}{g} \Bigl\{ \frac{d}{d\theta}\{\mathrm{cos}\theta\}\mathrm{sin}\theta + \mathrm{cos}\theta \frac{d}{d\theta}\{\mathrm{sin}\theta \} \Bigr\} =\frac{2{v_0}^2}{g} \Bigl\{ -\mathrm{sin}^2\theta + \mathrm{cos}^2\theta \Bigr\} $$
これが0のとき、
$$ 0 = -\mathrm{sin}^2\theta + \mathrm{cos}^2\theta $$
$$ \mathrm{sin}^2\theta = \mathrm{cos}^2\theta $$
$$ \mathrm{sin}\theta = \mathrm{cos}\theta $$
となります。このとき角度$ \theta $は
$$ \theta = \frac{\pi}{4} $$
のみです。
したがって角度が$ \frac{\pi}{4} $のとき水平到達距離は最大となります。
最大値は、
$$ x(\frac{\pi}{4})=\frac{2{v_0}^2}{g}\mathrm{cos}(\frac{\pi}{4})\mathrm{sin}(\frac{\pi}{4}) = \frac{{v_0}^2}{g} $$
となります。
余談ですがボールなり物を投げるとき角度45°がいいと言われますよね。それを上記の計算で確かめることができました。
春先は変人たちが踊りだす
最近暖かくなってきたなと感じます。気づけば3月中旬です、そりゃ暖かいわけだ。私の花粉症の症状も顕著に出てきてます。花粉やめてくれ~
ところで、春先は変質者だったり物騒な事件が多いような気がします。気体分子の運動エネルギーは絶対温度$ T $に比例しますよね。↓この式です。*1
$$ \frac{1}{2} m\overline {v^2} = \frac{3}{2} k T $$
冬から春に遷移していく中で気温差を感じます。冬場は寒くてあんまり元気がない人も春になると活発になるわけですね。そういうわけでこの春先は変質者とかが多いのかなと思いました。本州の人よりも高知の人とか九州、沖縄の人達ってなんか活発的だし変な人多いような気がします。これも向こうが暖かいからなんですかね。
今日はLondon方程式について取り上げてみようかなと思います。
London方程式
イギリスの首都の方程式ではありません。London兄弟が導いた方程式です。(イギリス出身じゃなくてドイツ出身ぽい)
これは超伝導体の性質であるMeisnner効果の現象論的方程式です。Meisnner効果ってのはアレです。超伝導体に外部磁場を加えたとき、超伝導体の表面に外部磁場を打ち消しあうような反磁性電流ってのが流れて、超伝導体内部には外部磁場が侵入しないってやつです。
超伝導体は電気抵抗が0 ΩになるのでOhmの法則が使えなくなります。そこでLondon兄弟は超伝導体の電流を電場由来ではなく磁場由来として考えました。へ~。*2
このLondon方程式を使った問題を見ていきます。
問題*3
解法としてまず思いつくのがLondon方程式とAmpere - Maxwellの式を用いて磁場$ \boldsymbol{B} $に関する式をたてることです。
また変位電流(厳密に言うと変位電流密度)$ \boldsymbol{j}_d $を無視できるので、ここではAmpereの式だけで十分です。
Ampereの式は、
$$ \mathrm{rot} \frac{\boldsymbol{B}}{\mu_0} = \boldsymbol{j} $$
です。
両辺に$ \mu_0 $をかけた後、回転$ \mathrm{rot} $を取ると、
$$ \mathrm{rot\ rot} \boldsymbol{B} = \mu_0 \ \mathrm{rot} \boldsymbol{j} $$
となります。
ところでベクトル三重積的な性質により、
$$ \mathrm{rot\ rot} \boldsymbol{B} = - \nabla ^2 \boldsymbol{B} + \mathrm{grad} \ \mathrm{div} \boldsymbol{B} = -\nabla ^2 \boldsymbol{B} $$
となります。
つまりAmpereの式は
$$ \nabla ^2 \boldsymbol{B} = - \mu_0 \ \mathrm{rot} \boldsymbol{j} $$
となります。こういう楕円型の2階偏微分方程式をPoisson方程式といいます。電磁気やってるとスカラーポテンシャルやベクトルポテンシャルのPoisson方程式をよく見かけます。というかこれもベクトルポテンシャルのPoisson方程式と同じです。
この磁場のPoisson方程式の右辺にLondon方程式を代入すると、
$$ \nabla ^2 \boldsymbol{B} = \frac{1}{\lambda ^2} \boldsymbol{B} $$
という磁場$ \boldsymbol{B} $だけの微分方程式ができました。次はこれを解いていきます。
問題の条件より変数$ z $に関する$ \boldsymbol{B}(z) $を考えていきます。つまり1変数の微分方程式に置き換えられるってわけですね。
簡単化した微分方程式は、
$$ \Bigl (\frac{d^2}{dz^2} - \frac{1}{\lambda^2} \Bigr) \boldsymbol{B}(z) =0 $$
となります。
今、この微分方程式のAnsatzを
$$ \boldsymbol{B}(z) = e^{\gamma z} $$
と置くことにします。
これを微分方程式に適用すると、
$$ \gamma^2 - \frac{1}{\lambda^2}=0 $$
という特性方程式が得られます。
これより適当に置いた指数の定数$ \gamma $は、
$$ \gamma = \pm \frac{1}{\lambda} $$
となります。さてこの$ \gamma $はプラスなのかマイナスなのか。
微分方程式の解が$ z \rightarrow \infty $としたとき発散しなければMeisnner効果を満たす(磁場が超伝導体中でちっこくなってほしい)ので、調べていきます。
$ \gamma = +\frac{1}{\lambda} $だと
$$ \lim_{z \to \infty } \boldsymbol{B}(z) = \lim_{z \to \infty } \exp \Bigl( {\frac{z}{\lambda}}\Bigr) = \infty $$
と明らかに発散します。つまり$ \gamma = - \frac{1}{\lambda} $が特性方程式の解になります。
また条件より$ \boldsymbol{B}(z=0) = \boldsymbol{B}_0 $だから結局のところ求める解は、
$$ \boldsymbol{B}(z) = \boldsymbol{B}_0 \exp \Bigl( {-\frac{z}{\lambda}}\Bigr) $$
となります。これが超伝導体内部の磁場の分布です。
変な不等式
以前花粉症の薬をもらいに耳鼻科に行ったときのことです。
耳鼻科の先生「この薬はね、効き目が悪いと思ったら倍量飲んでくださいね」
私「はい」
このとき私はn倍飲んでいいんだ、合法ODじゃんとか思ってました。インターネットで調べてみると2倍の量って書いてました。つまり2錠まで飲めるよってことでした。
さて今日も数学の話題についてです。昨日のJensenの不等式を使うと簡単に解ける問題です。
問題
(ヒント:$ x_i = \ln y_i - \ln z_i $とせよ)
$$ \prod \limits _{i=1} ^n {y_i ^{\alpha _i}} + \prod \limits _{i=1} ^n {z_i ^{\alpha _i}} \leq \prod \limits _{i=1} ^n {(y_i + z_i) ^{\alpha _i}} $$