Pika’sはきだめ

Twitter→@pikathatremains

たのしい!数値計算

この記事は、奈良高専 Advent Calendar 2022の24日目の記事です。前回の方(Jin)のブログはこちらです。

adventar.org

はじめに

皆さん初めまして。ぴかです。

最近BLEACHの千年血戦篇を見ています。めちゃくちゃおもろいですね。

↑バンビエッタちゃん。かわいい。

Fig.1 アマプラで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法

微分方程式が$ \frac{dy}{dx} = f(x,y) $のとき、解$ y(x) $を、
$$ y(x+\Delta x) = y(x) + f(x,y) \cdot \Delta x $$
により計算する。ここで$ \Delta x $は刻み幅である。

このEuler法によって式(1)で表される自由落下を数値計算します。なおPython数値計算しました。Pythonは個人的に使いやすいなっていうのと、colaboratoryとかweb上で実行できたりするので便利です(VScodeのほうが使いやすいけど)。

colab.research.google.com


コードはこんな感じ(きたないかもだけど許して。あとスマホだとこのコード表示が見にくくなるんやが、どうしたらいいんでしょうか。)。

#簡単な力学系のシミュ(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になる時刻で計算終了としました。

Fig.2 自由落下の初期条件

質点の高度変化をFig3に、速度変化をFig4に示します。

Fig.3 質点の高度変化(横軸は時間、縦軸は高度を表す)

Fig.4 質点の速度変化(横軸は時間、縦軸は速度を表す)

わ~、自由落下運動の数値解を求めれた~......
実はこれ、空気抵抗を考慮していないモデルなので、初期高度が高ければ高いほど落下時における質点の速度は速くなります。
例えば雨を降らす乱層雲は地上から1 km~2 kmぐらいの高さですが、このとき空気抵抗が無いと雨粒は140 ~198 m/sの速度で降ることになります(雨が10 kgの量降るとしたら)。
時速換算すると、504~713 km/h程度です(この速度で雨が降ればヤバイですよね)。

ということで空気抵抗を考慮した自由落下を考えましょう。

閑話休題

空気抵抗を考慮する前にちょっとだけ。

最近Polyphiaっていうバンドの曲をよく聞いてるんですが、まぁえぐいんですね。インストバンドなんですが、ギターがうますぎる。ギターはTim HensonとScott Lepageの2人なんですが、2人ともえげつうまい。

Fig.5 Tim Henson(Ibanezアー写

Fig.6 Scott LePage(Ibanezアー写

Playing Godとかえげついです。クラシックギターでピロピロ弾いてるし。ぜひ聞いてみてね。

www.youtube.com

open.spotify.com

空気抵抗を考慮した自由落下

物理の問題とかで速度に比例する空気抵抗とかよくありますよね*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に示します。

Fig.7 速度に比例する空気抵抗を考慮した自由落下での高度変化

Fig.6 速度に比例する空気抵抗を考慮した自由落下での速度変化

ちゃんと速度が一定のとこ(およそ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)

一昨日、全然寝れてなくて執筆・投稿が少し遅くなりました。ごめんなさい🎅

ちなみに「たのしい!数値計算」ってタイトルですが、アニメーションとか見てるときはたのしいけど、それ以外なんもたのしくなかったです、ねむすぎて。

*1:

*2:SIRモデルとかが有名です。ぜひ調べてみてください。

*3:変数分離形の微分方程式で、ごちゃごちゃ計算するやつ。微分方程式に慣れてくると、めちゃ簡単に爆速で解けるようになる。3年生以上なら物理の授業かなんかでやったことがあるはず。

n岡技科大 H31 電電

n岡技科大の電電を解いた。回路,電磁気関連しか解いていない。(電子回路とかは元気があればやるかも)

 

問題1 電気機器に関する問題

問1

4極の三相誘導電動機がある。この電動機に定格周波数50 Hzを加えて駆動しているときの同期速度$ N_0 $ [r/m]を求めよ。

極数を$ p $、周波数を$ f $とすると同期速度は、

$$ N_0 = \frac{120f}{p} $$

で表される。したがって、

$$ N_0 = \frac{120f}{p} = \frac{120\times 50}{4} = 1500  \rm{[r/m]} $$ 

問2

問1の電動機が回転速度$ N = 1440 $ [r/m]で回転しているとき、すべり$ s $を求めよ。

すべりは、

$$ s = \frac{N_0 - N}{N_0} = \frac{1500-1440}{1500}=0.04 $$

問3

問1の電動機の1相あたりのL型等価回路をFig.1-1のように求めたい。無負荷試験と拘束試験によって以下の結果を得ている。
(無負荷試験) 定格電圧$ 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} $]
を求めよ。

Fig.1-1

(1) 無負荷試験のとき、すべりが無い。したがって$ r_2^{'}\frac{1-s}{s} =\infty $となり、励磁回路のみ残る。等価回路はこんな感じ。

Fig.1-2

簡単な回路にできたので、無負荷試験の条件をそのまま解くと$ 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 $となる。また、電流のほとんどが負荷に流れるため励磁回路は無視できる。等価回路はこんな感じ。

Fig.1-3

拘束試験の条件より簡単に$ 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

Fig.2-1の回路について、角周波数$ \omega $ [rad/s]の交流入力電圧$ e_{in} $ [V]から見たインピーダンス$ Z_{all} $ [$ \rm{\Omega} $]を求めよ。但し、$ R_2 $は十分大きいものとする。

Fig.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

問1において、インピーダンスが最大となる角周波数を求めよ。

角周波数$ \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} $ [V]、出力電圧$ e_{out} $ [V]として、ラプラス変換を用いてブロック線図で表したところ、Fig2-2となった。Fig.2-2の(1)~(4)に入る伝達関数を示せ。

Fig.2-2

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

問3で求めたブロック線図をまとめて、入力電圧$ E_{in}(s) $から出力電圧$ E_{out}(s) $までの伝達関数を求めよ。

伝達関数は、

$$ 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

Fig.2-1の回路に、入力電圧$ e_{in} = 10 $ [V]の直流電圧を加えたときの出力電圧$ e_{out} $ [V]の過渡応答式を求めよ。但し、キャパシタ$ C $の初期電荷およびインダクタ$ L $の初期電流を0とし、$ R_1 =5 \rm{[\Omega]} $、$ R_2 =1 \rm{[\Omega]} $、$ C =0.2 \rm{[F]} $、$ L =1 \rm{[H]} $とする。
回路定数を問4で求めた伝達関数に適用すると、
$$ G(s) = \frac{s}{s^2 +6s +5 } $$
となる。この伝達関数の分母を因数分解すると、
$$ s^2 + 6s + 5 = (s+1 )  (s+5  ) $$
となる。
この伝達関数を逆ラプラス変換すると、
$$ \mathcal{L}^{-1}[G(s)] = g(t) = \frac{1}{4} \Bigl\{ 5e^{-5t} - e^{-t} \Bigr \} $$
となる。
つまり過渡応答式は、
$$ e_{out}(t) = 10g(t) = \frac{5}{2}\Bigl\{ 5e^{-5t} - e^{-t} \Bigr \} $$
となる。
 

問題3 電磁誘導に関する問題

磁場に垂直な面内で円板導体が回転しているため、円板導体中には起電力が生じる。まず、円板導体内の直線OC上でOから距離$ r $の位置にある微小な長さ$ dr $の半径要素に注目する。円板導体が回転する角速度が$ \omega $であるため、半径要素$ dr $が運動する速さは、$ v =(1) $である。これより、半径要素$ dr $の部分で発生する起電力の大きさは、$ de = (2) $となる。したがって、円板導体内の直線OCの部分において発生する起電力の大きさは、$ e = (3) $となる。円板導体中に起電力が発生することにより、抵抗$ R $には電流$ I $がFig.3-1中の$ (4) $に流れる。点Oと点Cにおける接触抵抗や円板導体における抵抗を無視できるとすると、定常状態においては、$ I = (5) $である。ここで求めた電流$ I $は円板導体内の直線OC間にも流れるので、回転する円板導体には力が働く。半径要素$ dr $の部分に働く力の大きさは、$ dF = (6) $となり、$ dF $の回転中心軸まわりのトルクの大きさは、$ dT = (7) $となる。したがって、円板導体全体に対して働く回転中心軸まわりのトルクの大きさは、$ T = (8) $となる。このトルク$ T $は、円板導体の回転を$ (9) $さえる方向に働く。
(1)~(9)に入る数式、語句を埋めよ。

Fig.3-1

円運動するときの半径$ 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:共振周波数じゃん!!

*2:ブロック線図をうまくまとめるか、関係式を根気よく解く。私は代数的に解くほうが好き。

*3:こういのを誘導電場という

どきどき物理

こんにちは。あったかくなりましたね。

突然ですがみなさんは最近ドキドキしましたか?
恋人と遊んだり、彼女に包丁を向けられたり、狭い道で対向車のトラックがセンターライン超えてたり、物理を勉強したり...

いろんなドキドキがあります。

さて今日は力学の問題を解いていきます。

問題1

運動方程式の解が$ x=(v_0 \mathrm{cos}\theta )t $、$ y=0 $、$ z=(v_0 \mathrm{sin}\theta)t-\frac{1}{2}gt^2 $で与えられているとする。また運動の軌跡は$ z=(\mathrm{tan}\theta)x-\frac{g}{2{v_0}^2\mathrm{cos}^2\theta}x^2 $である。このときの最高地点の高さ$ z $と到達時間を求めよ。

真っ先に思いつくのが高さ$ 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

運動方程式の解が$ x=(v_0 \mathrm{cos}\theta )t $、$ y=0 $、$ z=(v_0 \mathrm{sin}\theta)t-\frac{1}{2}gt^2 $で与えられているとする。また運動の軌跡は$ z=(\mathrm{tan}\theta)x-\frac{g}{2{v_0}^2\mathrm{cos}^2\theta}x^2 $である。このときの水平到達距離と到達時間を求めよ。また$ v_0 $が一定の場合、水平到達距離の最大値とそのときの角度$ \theta $はいくらか。

これも問題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方程式

真空透磁率を$ \mu_{0} $、外部磁場を$ \boldsymbol{B} $、超伝導体への磁場侵入長を$ \lambda $とすると電流密度$ \boldsymbol{j} $は、

$$ \nabla \times \boldsymbol{j} = - \frac{1}{\mu_{0} \lambda ^2} \boldsymbol{B} $$

という関係式で表される。

このLondon方程式を使った問題を見ていきます。

問題*3

超伝導体ではOhmの法則が成立せず、London方程式が成立すると考える。このとき$  xy $面を表面として$ z \geq 0 $の領域に配置された一様な超伝導体内部の$ \boldsymbol{B} $について、磁場$ \boldsymbol{B}(z) $の分布を求めよ。但しここで$ \boldsymbol{B}(z=0) = \boldsymbol{B}_0 $とし、また変位電流$ \boldsymbol{j}_d $は無視でき、超伝導体内の透磁率は$ \mu_0 $とする。

解法としてまず思いつくのが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) $$

となります。これが超伝導体内部の磁場の分布です。

 

 

*1:厳密にいうと気体分子の平均運動エネルギー

*2:このLondon方程式が成立するのは厳密に言うと電気抵抗が完全に0の超伝導体内部だけだそうです。

*3:早〇田先進理工の院試

変な不等式

以前花粉症の薬をもらいに耳鼻科に行ったときのことです。

耳鼻科の先生「この薬はね、効き目が悪いと思ったら倍量飲んでくださいね」

私「はい」

このとき私はn倍飲んでいいんだ、合法ODじゃんとか思ってました。インターネットで調べてみると2倍の量って書いてました。つまり2錠まで飲めるよってことでした。

youtu.be

さて今日も数学の話題についてです。昨日のJensenの不等式を使うと簡単に解ける問題です。

pika-remains.hatenablog.jp

 

問題

関数$ f(x) = \ln {(1+\exp x)} $は任意の実数$ x $に対して$ \frac{d^2 f(x)}{dx^2} \geq 0 $を満たす。この事実とJensenの不等式を用いて、$\alpha _i \geq 0 , $ $ \sum\limits  _{i=1} ^n \alpha _i =1 $を満たす任意の$ \alpha _1, ......,\alpha _n $と任意の正数$ y_i $、$ z_i $、$ i = 1, ...... , n $に対して以下の不等式が成立することを示せ。
(ヒント:$ 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}} $$
この不等式を示せってだけ書かれてたら難易度爆上がりだなとか感じます。うれしいことにヒントとJensenの不等式を使えるので存分に使っていきます。
 
まずJensenの不等式がどんなのだったか見てきます。
$$ \sum \limits _{i=1} ^n  {\alpha _i f(x_i)} \geq  f( \sum \limits _{i=1} ^n  {\alpha _i x_i}) $$
こんなんでした。
続いてJensenの不等式の左辺を計算します。
\begin{align} \sum \limits _{i=1} ^n  {\alpha _i f(x_i)} & =   \sum \limits _{i=1} ^n  {\alpha _i \ln \Bigl(1+ \exp[x_i]\Bigr)} \\\ & =  \sum \limits _{i=1} ^n  {\alpha _i \ln \Bigl(1+\exp[\ln \frac{y_i}{z_i}] \Bigr)}  \\\ & = \sum \limits _{i=1} ^n  {\ln \Bigl(1+ \frac{y_i}{z_i}\Bigr)^{\alpha _i}}  \end{align}
となります。
ここで総和と総乗(総積ともいうっぽい)の関係について触れます。
総和は
$$  \sum \limits _{i=1} ^n  {\alpha _i } = \alpha_1 + \alpha_2 + \alpha_3 +......+\alpha_n  $$
という感じで足し算をまとめたものでした。
総乗は
$$  \prod \limits _{i=1} ^n  {\alpha _i } = \alpha_1  \alpha_2 \alpha_3 \cdot \  ......\ \cdot \alpha_n  $$
みたいに掛け算をまとめたものです。
さて$ \log $の性質として足し算を掛け算にできるみたいなやつがありました。
$$  \log(\alpha) + \log(\beta) = \log(\alpha \beta) $$
この性質より左辺は結局のところ
 \begin{align} \text{(左辺)}& = \sum \limits _{i=1} ^n  {\ln \Bigl(1+ \frac{y_i}{z_i}\Bigr)^{\alpha _i}} \\\  & =\ln{\prod \limits _{i=1} ^n  {\Bigl(1+ \frac{y_i}{z_i}\Bigr)^{\alpha _i}}} \end{align}
となります。
次にJensenの不等式の右辺を計算します。
\begin{align}  f( \sum \limits _{i=1} ^n  {\alpha _i x_i}) &= \ln\Bigl(1+\exp[\sum \limits _{i=1} ^n  {\alpha _i x_i}]\Bigr) \\\ &= \ln\Bigl(1+\exp[\sum \limits _{i=1} ^n  {\ln(\frac{y_i}{z_i})^{\alpha_i}}]\Bigr) \\\ &= \ln\Bigl(1+\exp[\ln \prod_{i=1}^n {(\frac{y_i}{z_i})^{\alpha_i}}]\Bigr) \\\ & = \ln\Bigl(1+\prod_{i=1}^n {(\frac{y_i}{z_i})^{\alpha_i}}\Bigr) \end{align}
となる。
まとめると
$$ \ln{\prod \limits _{i=1} ^n  {\Bigl(1+ \frac{y_i}{z_i}\Bigr)^{\alpha _i}}} \geq \ln\Bigl(1+\prod_{i=1}^n {(\frac{y_i}{z_i})^{\alpha_i}}\Bigr) $$
$$ \prod \limits _{i=1} ^n  {\Bigl(1+ \frac{y_i}{z_i}\Bigr)^{\alpha _i}} \geq 1+\prod_{i=1}^n {(\frac{y_i}{z_i})^{\alpha_i}} $$
両辺に$ \prod \limits _{i=1}^n {z_i^{\alpha_i}} $をかけると
$$ \prod \limits _{i=1} ^n {(y_i + z_i) ^{\alpha _i}} \geq  \prod \limits _{i=1} ^n {y_i ^{\alpha _i}} + \prod \limits _{i=1} ^n {z_i ^{\alpha _i}}  $$ 
となる。
以上より不等式が成立することを示した。
 
この不等式ってなんて名前なんですかね。教えてエロい人。