ゆるふわブログ

東京大学理科 I 類 2 年の学生です.プログラミング,大学の勉強,日常生活で感じたことをゆるふわに書いていけたらなと思います.技術的に拙いところがあっても温かい目で見守っていただければ幸いです.

sinh 型の置換積分

A セメスター期末試験まであと 1 週間,今日で大学 1 年生の全講義が終了しました.感慨深いものがあります.
テスト勉強で忙しいので手軽に書ける記事を.

電磁気の過去問を同じクラスの強い人が解いて答えを上げてくれたのですが,そこで,あまり見る頻度の少ない置換積分をしていたのでメモっとこうと.
直線電流の周りの磁場を求める積分です.アンペールの法則使えば一瞬とか言いっこなし.

 \displaystyle\begin{equation*}
  \frac{\mu_0 I}{4\pi}\int_{-\infty}^\infty \frac{R}{\sqrt{R^2+z^2}^3}\,dz
\end{equation*}
僕は  \displaystyle z=R\tan\theta,\,dz = \frac{R}{\cos^2\theta}\,d\theta と置換していました.
 \begin{align*}
  \frac{\mu_0 I}{4\pi}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{R}{\sqrt{R^2(1+\tan^2 \theta)}^3}\, \frac{R}{\cos^2 \theta}d\theta &= \frac{\mu_0 I}{4\pi R}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos \theta\,d\theta \\
&= \frac{\mu_0 I}{4\pi R} \left[\sin \theta\right]_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \\
&= \frac{\mu_0 I}{2\pi R}
\end{align*}
しかし,  z=R\sinh t と置換することもできます.
 \begin{align*}
  \frac{\mu_0 I}{4\pi}\int_{-\infty}^\infty \frac{R}{\sqrt{R^2(1+\sinh^2 t)}^3}R\cosh t\,dt &= \frac{\mu_0 I}{4\pi R}\int_{-\infty}^\infty \frac{1}{\cosh^2 t}\,dt \\
&= \frac{\mu_0 I}{4\pi R} \left[\tanh t\right]_{-\infty}^\infty \\
&= \frac{\mu_0 I}{2\pi R}
\end{align*}
これは  \displaystyle \left(\sinh^{-1} \frac{t}{a}\right)' = \frac{1}{\sqrt{a^2+t^2}} に由来するものです.
式変形の途中双曲線関数の公式  \cosh^2 t-\sinh^2 t=1 を利用しました.
また, \displaystyle\lim_{t \to \pm\infty} \tanh t = \lim_{t \to \pm\infty} \frac{e^t-e^{-t}}{e^t+e^{-t}} = \pm1,\, (\tanh t)' = \frac{1}{\cosh^2 t} です.

以下が参考になるかも.

mathtrain.jp
mathtrain.jp

二次元波動のシミュレーション

物理とプログラミングを組み合わせられれば,僕の興味と合致して大変よろしいのではないかと思い至り,物理シミュレーションについてググったところ,以下のページを見つけました.


物理シミュレーションを通して学ぶプログラミング


うちの駒場キャンパスで 2013 年の冬学期でやっていたよう.僕もこんな心躍る授業を取りたかった.(シラバスを隅から隅まで見たつもりですが,今年はあまり面白い講義は開講されてませんでした…)

この講義では,3 次元の描画に VPython が使われていたようです.これほど簡単に 3D アニメーションを扱えるライブラリはあまりないらしいので,使ってみることにしました.

といったものの,どう使うのがうまいのかわからず,また Python もよく知らないので,わけのわからないものを作って遊んでました.


f:id:Ysmr_Ry:20170115011057g:plain


f:id:Ysmr_Ry:20170115011107g:plain


f:id:Ysmr_Ry:20170115030530g:plain


さすがにもうちょっとマシなものを作ろうと思い,何かネタはないかなと考えていたところ,昔とある先輩に教えてもらった物理シミュレーションをやってみようと思いました.それが波動のシミュレーションです.

元ネタはこの動画です.



当時高校生だった僕は,わかったようなわからないような感じでしたが,現在は振動・波動論を今期取っているので理解することは容易でした.解説していこうと思います.

水面波の浅い場合は,

 \displaystyle \begin{equation*} \frac{\partial^2 z(x,y,t)}{\partial t^2} = v^2\,\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)z(x,y,t) \end{equation*}

と近似できることは振動・波動の本に書いてあります.( 2 次元の場合.ただし, v=\sqrt{gh} )

これは,弦の振動や電磁波と同じ,二次元の波動方程式です.このシミュレーションでやることはこの方程式を伝搬のルールとして毎フレーム更新し,あとは好きなポイントを振動させてやれば,その振動がそれらしく伝わるというものです.

右,奥,上方向をそれぞれ x, y, z 軸として,2 次元配列 z[x][y] が点 (x,y) における z 座標を表すものとします.つまり,二次元配列の値をあたかも 2 変数関数 z(x,y) に見立てるということです.ここで気をつけたいのが,2 変数関数 z(x,y) の定義域,すなわち (x,y) の取りうる値は  \mathbb{R}^2 全体であるのに対し,z[x][y] の定義域は  \mathbb{Z}^2 であることです.つまり,関数の方では x, y は全実数という連続的な値を取れたのに対し,配列の方では x, y は整数という離散的な (とびとびな) 値しか取れないということです.実数は,"次の数" といったものがなく,どこまでも緻密に詰まっていますよね.0 の次の数として 1, 0.1, 0.01, 0.001 ... という風にいくらでも細かく出来てしまいます."次の数" の数との差はいわば無限小です.しかし,整数の場合は,"次の数" といったものが存在し,例えば 0 の次は 1 です."次の数" の数との差は必ず 1 になります.

波動方程式によると,z を x と y でそれぞれ偏微分しなければなりませんが,x と y が離散的な値しか取れない場合はどうすればよいでしょうか.ここで,先ほど考えた "次の数" との差といったものが効いてきます.偏微分と言っても所詮は 1 変数の常微分ですから,高校で習う微分の定義を思い出してみましょう.導関数は以下の極限で定義されます.

 \displaystyle \begin{equation*} f'(x) = \frac{df}{dx} = \lim_{h \to 0} \frac{f(x+h)-f(x)}{h} \end{equation*}

この式は,f(x) の変化量 f(x+h)-f(x) を 実数 x の "次の数" との差である無限小 h で割ったものであると見ることができます.これは,Leibniz の記法  \frac{df}{dx} にも表れています.微小変化量 df を微小な dx で割っているというイメージです.つまり,微分とは今いるところと "次の数" との差分進んだところの関数の値の差を "次の数" との差で割ったものと言えます.

このルールに従うと,離散バージョンの微分は次のように定義できるでしょう.( f が離散的な定義域を取ることがわかるようにあえて下付きの添字で  f_n と表記します.)

 \displaystyle \begin{equation*} \Delta f_n = \frac{f_{n+1} - f_n}{1} \end{equation*}

上にすでに書いてしまいましたが, f_n は数列, \frac{f_{n+1}-f_n}{1} はいわゆる階差数列であり,物理等で  \Delta f と表記します.つまり,差分は微分の離散バージョンだったのです.同様に, \sum{  } ("和分" と呼んだ方がしっくりくるかもしれない) は積分の離散バージョンです.形式的には, d \Delta に変えるだけです. \sum{ } が離散的な, \int{ } が連続的な和を表し,いずれも Sum の S を表す文字であることは聞いたことがあるのではないでしょうか. dx \int{ } で足しあわせて  \int{ } dx と表記するというのは,大学で物理を学べばしばしば出てくると思います.

以上から,波動方程式の右辺の 2 階偏微分を離散的な定義域を持つ関数に実行するには,階差数列の階差数列を考えれば良いことがわかります.すなわち,

 \begin{align*} \left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)z_{x,y} &= \Delta_x^2 z_{x,y} + \Delta_y^2 z_{x,y} \\ &= \Delta_x(z_{x+1,y}-z_{x,y}) + \Delta_y(z_{x,y+1}-z_{x,y}) \\ &= ( (z_{x+2,y}-z_{x+1,y})-(z_{x+1,y}-z_{x,y}) ) + ( (z_{x,y+2}-z_{x,y+1})-(z_{x,y+1}-z_{x,y}) ) \end{align*}

と計算できます*1*2.( z は例によって数列表記にしました)

これに伝わる速さの 2 乗を掛けたものが加速度になるので,点 (x,y) について速さ v と変位 z の 2 次元配列を用意しておき,速さに加速度を,変位に速さを足し合わせれば OK です.このとき,適当な減衰定数をかけると良いらしいです.

for x in range(0,sz):
  for y in range(0,sz):
    if x-1 >= 0 and x+1 < sz and y-1 >= 0 and y+1 < sz:
      acc = 0.10*(((z[x+1][y]-z[x][y])-(z[x][y]-z[x-1][y]))+((z[x][y+1]-z[x][y])-(z[x][y]-z[x][y-1])))
      vel[x][y] = (vel[x][y]+acc)*0.999
      #print "(%s,%s): acc=%s, vel=%s, z=%s" % (x,y,acc,vel[x][y],z[x][y])

for x in range(0,sz):
  for y in range(0,sz):
    z[x][y] += vel[x][y]

このような更新を毎フレームしつつ,クリックしたところに対応する点の変位を 10 に引っ張り上げることで,波が伝搬し始めます.マウス位置は VPython の scene.mouse.project で取得できます.

if scene.mouse.events and scene.mouse.getevent().press:
    mp = scene.mouse.project(normal=(0,1,0))
    mx = int(mp.x+sz/2)
    my = int(mp.z+sz/2)
    if mx >= 0 and mx < sz and my >= 0 and my < sz:
      z[mx][my] = 10

波が伝わる媒質は VPython の faces を組み合わせて作りました.
粗くて微妙なんですが,細かくしすぎると重いので… なんかいい手はないものか…
ちなみに,色はサーモグラフィーみたいになればいいなと思って変位の大きさによって変わるようにしました.color.hsv_to_rgb なんて便利なものが.(最初自前で変換してた…)

def col(x,y):
  return color.hsv_to_rgb((max(0.6-(h(x,y)/2 if h(x,y)>0 else 0),0), 0.8, 1))

for x in range(0,sz):
  for y in range(0,sz):
    fs.append(pos=(x*b-w/2,h(x,y),y*b-w/2), normal=(0,1,0), color=col(x,y))
    fs.append(pos=(x*b-w/2+b,h(x+1,y),y*b-w/2), normal=(0,1,0), color=col(x+1,y))
    fs.append(pos=(x*b-w/2+b,h(x+1,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x+1,y+1))
    fs.append(pos=(x*b-w/2+b,h(x+1,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x+1,y+1))
    fs.append(pos=(x*b-w/2,h(x,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x,y+1))
    fs.append(pos=(x*b-w/2,h(x,y),y*b-w/2), normal=(0,1,0), color=col(x,y))

fs.make_normals()
fs.make_twosided()
fs.smooth()

Python 初めて使ったし全く真面目に勉強していないので,見せられたもんじゃないですが,一応ソースの全貌.
2 次元配列最初 [[0.]*sz]*sz とかやっちゃって Ruby の object_id がどうのこうのってあの意味不明な挙動と同じアレに引っかかりました…*3 内包表記使おう!
あと,Python の論理演算子って and / or / not なんですね…

from visual import *
from math import *

b = 1
w = 30
sz = int(w/b)

z = [[0 for i in range(sz)] for j in range(sz)]
vel = [[0 for i in range(sz)] for j in range(sz)]

tim = 0

fs = faces()

for x in range(0,sz):
  for y in range(0,sz):
    fs.append(pos=(x*b-w/2,0,y*b-w/2), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2+b,0,y*b-w/2), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2+b,0,y*b-w/2+b), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2+b,0,y*b-w/2+b), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2,0,y*b-w/2+b), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2,0,y*b-w/2), normal=(0,1,0), color=(0,128,128))

while True:
  rate(60)

  if scene.mouse.events and scene.mouse.getevent().press:
    mp = scene.mouse.project(normal=(0,1,0))
    mx = int(mp.x+sz/2)
    my = int(mp.z+sz/2)
    if mx >= 0 and mx < sz and my >= 0 and my < sz:
      z[mx][my] = 10

  tim += 1

  for x in range(0,sz):
    for y in range(0,sz):
      if x-1 >= 0 and x+1 < sz and y-1 >= 0 and y+1 < sz:
        acc = 0.10*(((z[x+1][y]-z[x][y])-(z[x][y]-z[x-1][y]))+((z[x][y+1]-z[x][y])-(z[x][y]-z[x][y-1])))
        vel[x][y] = (vel[x][y]+acc)*0.999
        #print "(%s,%s): acc=%s, vel=%s, z=%s" % (x,y,acc,vel[x][y],z[x][y])

  for x in range(0,sz):
    for y in range(0,sz):
      z[x][y] += vel[x][y]

  fs.visible = False
  del fs
  fs = faces()

  def col(x,y):
    return color.hsv_to_rgb((max(0.6-(h(x,y)/2 if h(x,y)>0 else 0),0), 0.8, 1))

  def h(x,y):
    return z[x][y] if x >= 0 and x < sz and y >= 0 and y < sz else 0

  for x in range(0,sz):
    for y in range(0,sz):
      fs.append(pos=(x*b-w/2,h(x,y),y*b-w/2), normal=(0,1,0), color=col(x,y))
      fs.append(pos=(x*b-w/2+b,h(x+1,y),y*b-w/2), normal=(0,1,0), color=col(x+1,y))
      fs.append(pos=(x*b-w/2+b,h(x+1,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x+1,y+1))
      fs.append(pos=(x*b-w/2+b,h(x+1,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x+1,y+1))
      fs.append(pos=(x*b-w/2,h(x,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x,y+1))
      fs.append(pos=(x*b-w/2,h(x,y),y*b-w/2), normal=(0,1,0), color=col(x,y))

  fs.make_normals()
  fs.make_twosided()
  fs.smooth()
  
  if scene.kb.keys:
    if scene.kb.getkey() == "z":
      for x in range(0,sz):
        for y in range(0,sz):
          z[x][y] = vel[x][y] = 0

これからも VPython 使って物理と絡めていきたいです…

以下動いてる様子.

f:id:Ysmr_Ry:20170115025951g:plain

[追記]
以前はそんなことなかったような気がしますが,今は 波動 シミュレーションで調べると結構出てきて n 番煎じ感が半端ないですねすみません…
その中の一つ,

畳み込み演算を利用して波動方程式をシミュレーション - Qiita

を見てちょっとおもしろいと思いました.この記事では,波動方程式の解くのに行列を使っています.僕もやってみます.
簡単のために,数列に時刻を表す添字 n を加えて  z^n_{x,y} と表します.そして,行列  Z^n

 \begin{equation*}
Z^n = \begin{pmatrix}
  z^n_{0,0} & z^n_{1,0} & \cdots & z^n_{N_x,0} \\
  z^n_{0,1} & z^n_{1,1} & \cdots & z^n_{N_x,1} \\
  \vdots & \vdots & \ddots & \vdots \\
  z^n_{0,N_y} & z^n_{1,N_y} & \cdots & z^n_{N_x,N_y}
\end{pmatrix}
\end{equation*}
と定義します.つまりこの行列はある時刻 n における各座標での変位,つまり媒質の状態を表しています.時刻の 2 階微分も差分に対応させると,波動方程式は,
 \begin{align*}
  \delta^2 Z^n &= v^2 (Z^n * K) \\
  Z^{n+1}-2Z^n+Z^{n-1} &= v^2 (Z^n * K) \\
  Z^{n+1} &= 2Z^n-Z^{n-1}+v^2(Z^n * K) \\
  \, \\
  K &= \begin{pmatrix}
    0 & 1 & 0 \\
    1 & -4 & 1 \\
    0 & 1 & 0
  \end{pmatrix}
\end{align*}
と行列の三項間漸化式になります.(中心差分は n について)
* は畳み込み演算で,行列の各位置に対してフィルターの中心の成分を対応させ,重なった成分をそれぞれ掛け足し合わせる演算です.画像処理のフィルターの演算に使うらしい.
ラプラシアンフィルター出てきた.
登場する文字を行列という 2 次元的な広がりを持つ数にしたことで x, y についての空間に関する微分が畳み込み演算という単なる行列の演算になり,残る t 一文字に関する微分だけが残り,単なる差分方程式 (数列の漸化式) になったということだと思います.(たぶん) ラプラシアンフィルターはまさにラプラシアン (各方向の 2 階微分の和) を行列の四角形としての広がりに落とし込んだ畳み込みカーネルなのです.
このサイト (たたみ込み演算による画像処理) によると,1 次微分のソーベルフィルターは向きによって違うのに対し 2 次のラプラシアンフィルターは濃度勾配の強さだけを求めるのだという.
1 次元のダランベールの解を構成する 2 つの関数が +x と -x 方向を表す波であり,進行方向がどちらでも満たすように 2 階偏微分になっているという話を彷彿とさせますね!
もし, Z^n を各ピクセルの色情報を持った行列だとみなすと,ぱっと見で輪郭をなしているところほど加速度が速いことになります.
濃度勾配の大きく境界となっているところほど,その外側が強い力を受けて励振されるということなのでしょうか… (適当)
「輪郭抽出に使うフィルターをかけるとあら不思議波になる!」ではなく,波の伝搬の仕組みをよくよく理解すると,いかに輪郭っぽいかという度合いが算出できるってことなのかなぁ…
よくわかりません()
行列で描き直してみた.

from visual import *
from math import *
import numpy as np
from scipy import signal

b = 1
w = 30
sz = int(w/b)

z = np.array([[0. for i in range(sz)] for j in range(sz)])
z_p = np.array([[0. for i in range(sz)] for j in range(sz)])
z_pp = np.array([[0. for i in range(sz)] for j in range(sz)])
lapK = np.array([[0.,1.,0.],[1.,-4.,1.],[0.,1.,0.]])

fs = faces()

for x in range(0,sz):
  for y in range(0,sz):
    fs.append(pos=(x*b-w/2,0,y*b-w/2), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2+b,0,y*b-w/2), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2+b,0,y*b-w/2+b), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2+b,0,y*b-w/2+b), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2,0,y*b-w/2+b), normal=(0,1,0), color=(0,128,128))
    fs.append(pos=(x*b-w/2,0,y*b-w/2), normal=(0,1,0), color=(0,128,128))

while True:
  rate(60)

  if scene.mouse.events and scene.mouse.getevent().press:
    mp = scene.mouse.project(normal=(0,1,0))
    mx = int(mp.x+sz/2)
    my = int(mp.z+sz/2)
    if mx >= 0 and mx < sz and my >= 0 and my < sz:
      z[mx, my] = 5

  z = (2.0*z_p-z_pp+0.1*signal.correlate2d(z_p, lapK, 'same'))*0.99
  z_pp = z_p
  z_p = z

  fs.visible = False
  del fs
  fs = faces()

  def col(x,y):
    return color.hsv_to_rgb((max(0.6-(h(x,y)/2 if h(x,y)>0 else 0),0), 0.8, 1))

  def h(x,y):
    return max(-10, min(z[x, y],10)) if x >= 0 and x < sz and y >= 0 and y < sz else 0

  for x in range(0,sz):
    for y in range(0,sz):
      fs.append(pos=(x*b-w/2,h(x,y),y*b-w/2), normal=(0,1,0), color=col(x,y))
      fs.append(pos=(x*b-w/2+b,h(x+1,y),y*b-w/2), normal=(0,1,0), color=col(x+1,y))
      fs.append(pos=(x*b-w/2+b,h(x+1,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x+1,y+1))
      fs.append(pos=(x*b-w/2+b,h(x+1,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x+1,y+1))
      fs.append(pos=(x*b-w/2,h(x,y+1),y*b-w/2+b), normal=(0,1,0), color=col(x,y+1))
      fs.append(pos=(x*b-w/2,h(x,y),y*b-w/2), normal=(0,1,0), color=col(x,y))

  fs.make_normals()
  fs.make_twosided()
  fs.smooth()
  
  if scene.kb.keys:
    if scene.kb.getkey() == "z":
      for x in range(0,sz):
        for y in range(0,sz):
          z[x, y] = z_p[x, y] = z_pp[x, y] = 0

こっちの方がそれっぽいかな.

f:id:Ysmr_Ry:20170117083949g:plain

*1:ここでは, \Delta f_n = f_{n+1} - f_n で定義しました.これを前進差分といいます.

差分にはいくつか種類があり,以下のように定義されます.
 \begin{align*} 
\Delta f_n &= f_{n+1} - f_n \\
\nabla f_n &= f_n - f_{n-1} \\
\delta\,f_n &= f_{n+\frac{1}{2}}-f_{n-\frac{1}{2}}
\end{align*}
これらをそれぞれ 前進 / 後退 / 中心 差分といいます. ここでは最初に挙げた微分の定義が前進型だったために前進差分で話を進めていますが,後退 / 中心型でも問題ありません.
(中心型は  O(h^2) の勢いで収束するらしい)
ソースコードは説明の中での二階差分を
 \begin{equation*}
\Delta_i^2 = \delta^2 = \Delta\nabla = \nabla\Delta \,\, (i=x,y)
\end{equation*}
と捉えています.(中心差分 2 回が最もわかりやすいかも)
また,係数の 1, -2, 1 は二項係数であり,n 階差分の係数は  (x-1)^n を展開した係数に相当します.
このような微分方程式の数値的解法のことを有限差分法というそうです.下記を参考にしました.

有限差分 - Wikipedia
差分法 - Wikipedia

*2:これは画像の輪郭を抽出するラプラシアンフィルターに関連があるらしい (動画参照)

*3:例えば a = [[0]*4]*4 として a[0][0] = 1 とすると [[1,0,0,0],[1,0,0,0],[1,0,0,0],[1,0,0,0]] となる.

このブログについて

高校時代の後輩に勧められて,ブログを始めてみることにしました.

(といっても,どれくらい続くかわかりませんが…)

 

プログラミング系のことや,大学での勉強について,また日常生活で感じたことなどを綴っていこうと思います.ものの考え方や日本語がおかしかったらゴメンなさい…

 

当初は,僕程度が技術的な記事など書いて何になるのかと思いましたが,自分に対する記録,また,文書にまとめ直すことで頭の中を整理し理解を深めるという意味で,ありなのかなと思ったので始めます.(実はすっかり忘れていたのですが先ほど後輩から催促されまして…)

 

技術的に拙いところ等多々あるかと思いますが,何分修行中の身なので,そこは意地悪せずに温かい目で見守っていただき,もしよろしければ,ご指摘をいただければ幸いです.

 

楽しむことを第一にやっていこうと思います.