読者です 読者をやめる 読者になる 読者になる

ゆるふわブログ

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

Lorentz 変換と双曲線

1 年生の夏学期,「基礎方程式とその意味」という講義で時期尚早にもやった Lorentz 変換.
導出の過程は追うことができましたが,直感的な理解ができずしっくりと来ていませんでした.
天下り的に双曲線関数を使うときれいになるということも,どのような背景があるのかさっぱりでした.

今日から,統計と多変量解析の本を読み終わってひと段落ついた (それについての記事はまた後日) ので,途中で挫折した相対論を勉強しようと本を読み始めました.
そして,このサイトを見つけてしまいました.

d.hatena.ne.jp

非常に美しく・わかりやすく Lorentz 変換を説明しています.
また,この資料でさらに理解が深まりました.

https://www.kaijo.ed.jp/wp-content/uploads/2016/02/2013summer-5_2.pdf

内積空間と絡めて説明しています.Lorentz 変換における不変量 (の絶対値) が Minkowski ノルムの 2 乗であることがわかります.

この感動を是非伝えるべく記事にまとめねば (笑) と思い,僕なりにまとめてみました.

github.com

数式混じりの記事は LaTeX にしたほうがいいという結論に至った.

LaTeX 練習

最近教習所に通ってあとはグダグダしててあんまりプログラミングしてないです()

今日あるの知らずにノー勉で突入した仮免学科試験が合格していたことがわかりました.
47/50 だったけど…

春休みはその他学習塾で個別指導のバイトをしていて,受験生が卒業してしまったので今持ってる生徒が高 2 の子 1 人なんですね.暇だし,1 人相手なら余裕あるし,何かできないかなと…

そこで,バイトの生徒のために LaTeX の練習がてら数学のプリントを作りました.TikZ パッケージで図も入れられるようになった.

github.com

数式のフォントは数学ガールのあれです.

何気に GitHub 使ったの初めてです… 使い方よくわかってない…

Haskell で ABC 001 を解く

qiita.com

この記事を読んだからやろうと思った.深い意味はない.

AtCoder Beginner Contest 001 - AtCoder Beginner Contest 001 | AtCoder

A - 積雪深差

main :: IO ()
main = do
  h1 <- readLn
  h2 <- readLn

  print $ h1-h2

B - 視程の通報

import Control.Applicative
import Text.Printf

main :: IO ()
main = do
  m <- (/1000) <$> readLn :: IO Double

  printf "%02.0f\n" $ cond m

cond :: Double -> Double
cond m
  | m < 0.1            = 0
  | m >= 0.1 && m <= 5 = m*10
  | m >= 6 && m <= 30  = m+50
  | m >= 35 && m <= 70 = (m-30)/5+80
  | m > 70             = 89

MultiWayIf 拡張が Compile Error で切れそうだった.

C - 風力観測

import Control.Applicative ((<$>))

main :: IO ()
main = do
  [deg, dis] <- map read . words <$> getLine
  let ansW = w $ roundAt 1 $ dis/60

  putStrLn $ if ansW == 0
    then "C 0"
    else dir (deg/10) ++ " " ++ show ansW

roundAt :: Int -> Double -> Double
roundAt d x
  | dec >= 0.5 = fromIntegral (int+1)/10^d
  | otherwise  = fromIntegral int/10^d
  where
    (int,dec) = properFraction $ x*10^d

dir :: Double -> String
dir deg
  | null cand = "N"
  | otherwise = snd . head $ cand
  where
    cand = filter (\(i,_) -> 11.25+i*22.5 <= deg && deg < 33.75+i*22.5) $ zip [0..] ["NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW"]

w :: Double -> Int
w dis = head $ filter (\i -> fst (lims!!i) <= dis && dis <= snd (lims!!i)) [0..12]
  where
    lims = [(0.0,0.2),(0.3,1.5),(1.6,3.3),(3.4,5.4),(5.5,7.9),(8.0,10.7),(10.8,13.8),(13.9,17.1),(17.2,20.7),(20.8,24.4),(24.5,28.4),(28.5,32.6),(32.7,200.0)]

Haskell の round は HALF_EVEN とかいうそうで偶数に丸めるらしく,挙動が望んだものと違ったので自分で実装せざるを得なかった.このサイトを参考にしました.
oropon.hatenablog.com

D - 感雨時刻の整理

import Control.Applicative ((<$>))
import Control.Monad (replicateM)
import Data.List (sort, splitAt)

type Section = (String, String)

main :: IO ()
main = do
  n <- readLn
  secs <- fmap sort $ replicateM n $ do
    (s,_:e) <- splitAt 4 <$> getLine
    return (roundTF s, roundTS e)

  mapM_ (\(s,e) -> putStrLn $ s ++ "-" ++ e) . reverse $ foldl f [] secs

  where
    f :: [Section] -> Section -> [Section]
    f [] s = [s]
    f (x:xs) s
      | x<&>s     = x<|>s:xs
      | otherwise = s:x:xs

show02d :: Int -> String
show02d i
  | i < 10    = "0" ++ show i
  | otherwise = show i

roundTF :: String -> String
roundTF t
  | m >= 60   = show02d (h+1) ++ show02d (m-60)
  | otherwise = show02d h ++ show02d m
  where
    roundI m = m`div`5*5

    h = read $ take 2 t
    m = roundI . read $ drop 2 t

roundTS :: String -> String
roundTS t
  | m >= 60   = show02d (h+1) ++ show02d (m-60)
  | otherwise = show02d h ++ show02d m
  where
    roundI m
      | m `mod` 5 /= 0 = (m`div`5+1)*5
      | otherwise      = m

    h = read $ take 2 t
    m = roundI . read $ drop 2 t

(<&>) :: Section -> Section -> Bool
(<&>) (a,b) (c,d) = c <= b

(<|>) :: Section -> Section -> Section
(<|>) (a,b) (c,d) = (min a c, max b d)

Data.List.Split が Compiler Error で splitOn "-" ってできなくて切れそうだった.

Haskell で競プロの Straightforward な問題をやるといい練習になる (?)

補足 - Text.Printf.printf の仕組み

printf はどうして可変長引数を取れるんでしょうか.
その秘密は PrintfType 型クラスにありました.
PrintfType 型クラスは,

class PrintfType r where
  printf :: String -> r

と定義されていて,PrintfType のインスタンスには,

instance PrintfType String
instance PrintfType (IO ())
instance (PrintfArg a, PrintfType r) => PrintfType (a -> r)

があります.最初の 2 つが printf の型の終端にあたります.
3 つ目のように,再帰的にインスタンス宣言することで,printf :: String -> a -> a -> … -> a -> (IO () or String) となります.よって,printf は IO アクションとしてだけでなく,C 言語でいう sprintf のように出力先が文字列でも使えます.
以下のサイトを参考にしました.
stackoverflow.com

連成振動と 2 次超曲面

テスト無事終わりました!
今日は振動・波動論のテスト中に気づいたことをまとめてみようと思います.

振動・波動論のテストの大問 1 は連成振動の問題でした.3 連成振子.3 つの質量 M のおもりが互いにバネで接続されているやつです.ここでは一般に N 連成振子を考えます.N 個のおもりの位置を  x_n\ (n=1,2,3,\cdots) とすると,運動方程式は,

 \begin{align*}
  M\ddot{x}_n &= -k(x_n-x_{n-1})+k(x_{n+1}-x_n)\\
  \ \\
  M \ddot{\mathbf{x}} &= -K\mathbf{x} \\
  \ \\
  \mathbf{x} &= \begin{pmatrix}
    x_1 \\
    x_2 \\
    \vdots \\
    x_N
  \end{pmatrix} \\
  \ \\
  K &= \begin{pmatrix}
    2k & -k & 0 & \cdots & \cdots & \cdots & 0 \\
    -k & 2k & -k & 0 & \cdots & \cdots & 0 \\
    0 & -k & 2k & -k & 0 & \cdots & 0 \\
    \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
    0 & \cdots & \cdots & 0 & -k & 2k & -k \\
    0 & \cdots & \cdots & \cdots & 0 & -k & 2k
  \end{pmatrix}
\end{align*}

となります. x_0 = x_{N+1} = 0 という形式的な変数を用意すれば,第 1 式のように綺麗に表されます.それを行列とベクトルでまとめたのが第 2 式です.規則性のおかげで, K は対称行列になります.よって,固有ベクトルを正規直交基底にしたものを束ねた直交行列で対角化できます.この直交行列を  P とすると,

 \begin{align*}
  P^{-1} K P &= \begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix} \\
  \ \\
  M\ddot{\mathbf{x}} &= -P\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix} P^{-1} \mathbf{x} \\
  \ \\
  MP^{-1}\ddot{\mathbf{x}} &= -\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix} P^{-1} \mathbf{x}
\end{align*}
よって,
 \begin{equation*}
  \mathbf{Q} = \begin{pmatrix}
    Q_1 \\
    Q_2 \\
    \vdots \\
    Q_N
  \end{pmatrix} = P^{-1} \begin{pmatrix}
    x_1 \\
    x_2 \\
    \vdots \\
    x_N
  \end{pmatrix} = P^{-1} \mathbf{x}
\end{equation*}
と基準座標を導入すると,
 \begin{align*}
  M\ddot{\mathbf{Q}} &= -\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix} \mathbf{Q} \\
  \ \\
  M\ddot{Q}_n &= -\lambda_n Q_n
\end{align*}
となり,各方程式は単振動型の微分方程式になります.よって, Q_n は解けるので,基準座標の定義を遡って  \mathbf{x} = P \mathbf{Q} と解が求められます.つまり,各  x_n は角振動数  \sqrt{\frac{\lambda_1}{M}},\ \sqrt{\frac{\lambda_2}{M}},\cdots, \sqrt{\frac{\lambda_N}{M}} の単振動の重ね合わせ (線形結合) で表されることがわかります.

つまり,1 つの固有値が 1 つの基準振動数  \omega_n = \sqrt{\frac{\lambda_n}{M}} に対応しているわけですが,固有ベクトルは対応する基準振動の振幅比を表します.
これは, P が直交行列なので  P^{-1} = {}^t\!P であることから,固有ベクトル \mathbf{v_n} とすると,

 \begin{align*}
  P &= \begin{pmatrix}
    \mathbf{v_1} & \mathbf{v_2} & \cdots & \mathbf{v_N}
  \end{pmatrix} \\
  \ \\
  \mathbf{Q} &= P^{-1}\mathbf{x} \\
  &= {}^t\!P\mathbf{x} \\
  &= \begin{pmatrix}
    \mathbf{v_1} \\
    \mathbf{v_2} \\
    \vdots \\
    \mathbf{v_N}
  \end{pmatrix} \begin{pmatrix}
    x_1 \\
    x_2 \\
    \vdots \\
    x_N
  \end{pmatrix}
\end{align*}

となり,基準座標  Q_n \mathbf{v_n} = (v_{n1}, v_{n2}, \cdots, v_{nN}) の各成分で  x_n を線形結合した  \displaystyle\sum_{i=1}^N v_{ni}x_i で表されるからです.

で,テストで,基準座標で全エネルギー  E を表しなさいって問題が出たんですね.まぁ,基準振動のエネルギーの和だろうから,

 \displaystyle \begin{equation*}
  E = \sum_{n=1}^N \frac{1}{2}M \left(\dot{Q}_n^2 + \omega_n^2 Q_n^2\right)
\end{equation*}

で終わりにしようと思ってたんですが,もうちょっとちゃんと書いといたほうがいいのかなと思って.プリントに書いてあったことを思い出してたんですが… プリントには,連成振動の全エネルギーをまず  x_n の方で表してから, Q_n に持って行くみたいな流れだったんです.で, x_n で全エネルギーを表すと,

 \begin{equation*}
  E = \frac{1}{2}M{}^t\!\dot{\mathbf{x}}\dot{\mathbf{x}} + \frac{1}{2} {}^t\!\mathbf{x}K\mathbf{x}
\end{equation*}

って書いてあったんですね! 弾性エネルギーが 2 次超曲面の形で書けちゃうんです!
(タイトル回収)
すごく綺麗だと思いました. \mathbf{x} を一つのデータと捉えた時,2 次曲面はまさにその "2 次" の形なのだなぁと.

なんでこうなるのかちょっと考えてみました.
 x_n運動方程式 \dot{x}_n をかけると,力 × 速度で仕事率になるので,

 \displaystyle \begin{align*}
  M \ddot{x}_n \dot{x}_n &= +kx_{n-1}\dot{x}_n -2kx_n\dot{x}_n +kx_{n+1}\dot{x}_n \\
  \frac{d}{dt}\left(\frac{1}{2}M\dot{x}_n^2\right) &= +kx_{n-1}\dot{x}_n -\frac{d}{dt}\left(\frac{1}{2}2kx_n^2\right) +kx_{n+1}\dot{x}_n
\end{align*}

これを n が 1 から N まで足し合わせると, x_0 = x_{N+1} = 0 は消えて,積の微分から  \displaystyle\dot{x}_nx_{n+1} + x_n\dot{x}_{n+1} = \frac{d}{dt}\left(x_nx_{n+1}\right) とまとめられるので,

 \displaystyle \begin{align*}
  \frac{d}{dt}\left(\sum_{n=1}^N \frac{1}{2}M\dot{x}_n^2\right) &= \frac{d}{dt}\left(\sum_{n=1}^N \left(k x_{n}x_{n+1}-\frac{1}{2}2kx_n^2\right)\right) \\
  &= -\frac{d}{dt}\left(\frac{1}{2}{}^t\!\mathbf{x}K\mathbf{x}\right) \\
  \frac{d}{dt}\left(\frac{1}{2}M{}^t\!\dot{\mathbf{x}}\dot{\mathbf{x}} + \frac{1}{2} {}^t\!\mathbf{x}K\mathbf{x}\right) &= 0
\end{align*}

となります. K は対称行列ですから,i < j を満たす (i,j) 成分 (対角成分の上側) と i > j を満たす (i,j) 成分 (対角成分の下側) が同じです.よって, \dfrac{1}{2} {}^t\!\mathbf{x}K\mathbf{x} \dfrac{1}{2} によって 1 つ分が残ります.また,対角成分は 1 つだけなので, \dfrac{1}{2} によって,半分になります.第 1 式の右辺を見ると, \dfrac{d}{dt} の中で, x_nx_{n+1} には  \dfrac{1}{2} がかかっていませんが, x_n^2 の項にはかかっています.よって,これは 2 次超曲面の形で書けるのです*1

これを  \mathbf{x} = P\mathbf{Q} を代入して基準座標の式に変えます.

 \displaystyle \begin{align*}
  \frac{1}{2}M{}^t\!\dot{\mathbf{x}}\dot{\mathbf{x}} + \frac{1}{2} {}^t\!\mathbf{x}K\mathbf{x} &= \frac{1}{2}M{}^t\!{(P\dot{\mathbf{Q}})}P\dot{\mathbf{Q}} + \frac{1}{2} {}^t\!{(P\mathbf{Q})}KP\mathbf{Q} \\
  &= \frac{1}{2}M{}^t\!\dot{\mathbf{Q}}P^{-1}P\dot{\mathbf{Q}} + \frac{1}{2} {}^t\!\mathbf{Q}P^{-1}KP\mathbf{Q} \\
  \ \\
  &= \frac{1}{2}M{}^t\!\dot{\mathbf{Q}}\dot{\mathbf{Q}} + \frac{1}{2} {}^t\!\mathbf{Q}\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_N
  \end{pmatrix}\mathbf{Q}
\end{align*}

これで,形は少し違いますが,上にあげた基準座標でのエネルギー  E の表現が導けました.
ここで注目したいのが 2 次曲面がどうなったかです.上の方で説明したように,対角化を連立微分方程式を解くのに利用していると捉えることもできますが, 2 次超曲面によって弾性エネルギーを表すと,対角化は 2 次曲面  {}^t\!\mathbf{x}K\mathbf{x} を標準形  \displaystyle \sum_{n=1}^N \lambda_n Q_n^2 に変形していると捉えることができるのです!対角化の応用例というと,この連成振動の話と,微分積分学線形代数学を大学初年次で習っている身としては,Hessian による 2 変数関数の極値判定をあげたくなりますが,あれも, f_x(x,y) = f_y(x,y) = 0 を満たす停留点の周りで 2 次近似をして,2 次曲面の標準形を求め, f(x,y) = x^2+y^2,\ f(x,y) = -x^2-y^2,\ f(x,y) = x^2-y^2 の場合に帰着していると直感的には捉えることができますから*2,この話と繋がっていたのです.一見バラバラだった対角化の応用例が,実は同じ話だったと自分の中で繋がったのでとてもうれしかったです.微分積分学線形代数学,振動・波動論,今期取っていた 3 つの授業で暗に同じことを繰り返し教わっていたのだなぁとしみじみ思いました.

*1:2 次超曲面では, K の (i,j) 成分が  x_ix_j の係数に対応します.

*2:つまり,2 次超曲面の係数行列が  \begin{pmatrix} f_{xx} & f_{xy} \\ f_{yx} & f_{yy} \end{pmatrix} となるということです.

ブロック崩し (by Haskell)

Haskell から JavaScript に変換してくれる Haste というライブラリを使ってとりあえずブロック崩しを作って見ました.
忙しいので詳細は後日書きます.
以下ソースコード

{-# LANGUAGE LambdaCase #-}
import Haste
import Haste.DOM
import Haste.Events
import Haste.Graphics.Canvas
import Lens.Micro
import Control.Applicative
import Control.Arrow ((***))
import Control.Monad
import Data.Array
import Data.IORef
import System.Random

data BlkState = NoDamage | Hit | Destroy
  deriving (Eq, Enum)

data ScrState = Prepare | Game | Clear
  deriving Eq

data State = State {
  scrRef :: IORef ScrState,
  keyFRef :: IORef (Bool, Bool),
  bxRef :: IORef Double,
  bpRef :: IORef Point,
  bvRef :: IORef Point,
  missRef :: IORef Int,
  blocks :: IORef (Array (Int, Int) BlkState)
}

bnw, bnh :: Int
(bnw,bnh) = (12, 12)

main :: IO ()
main = do
  elemById "canvas" >>= \case
    Just cid -> do
      w <- read <$> getProp cid "width"
      h <- read <$> getProp cid "height"
      let (bw, bh) = (floor (w/fi bnw), floor (h/fi (2*bnh)))

      scrRef <- newIORef Prepare
      keyFRef <- newIORef (False, False)
      bxRef <- newIORef (-50)
      bpRef <- newIORef (0, h/2-70)
      theta <- randomRIO (-pi/3, pi/3) :: IO Double
      bvRef <- newIORef (5*cos (theta-pi/2), 5*sin (theta-pi/2))
      missRef <- newIORef 0
      arr <- newIORef $ listArray ((0,0),(bh,bw)) $ repeat NoDamage

      documentBody `onEvent` KeyDown $ \(KeyData code _ _ _ _) -> do
        when (code==37) $ modifyIORef keyFRef (&_1.~True)        -- Left
        when (code==39) $ modifyIORef keyFRef (&_2.~True)        -- Right
        scr <- readIORef scrRef
        when (code==13 && scr /= Clear) $ writeIORef scrRef Game -- Enter
      documentBody `onEvent` KeyUp $ \(KeyData code _ _ _ _) -> do
        when (code==37) $ modifyIORef keyFRef (&_1.~False) -- Left
        when (code==39) $ modifyIORef keyFRef (&_2.~False) -- Right
      
      getCanvas cid >>= \case
        Just cvs -> loop cvs (w,h) (bw,bh) (State scrRef keyFRef bxRef bpRef bvRef missRef arr)
        Nothing  -> error "Canvas could not be found!"
    Nothing  -> error "Canvas ID could not be found!"

loop :: Canvas -> Vector -> (Int, Int) -> State -> IO ()
loop cvs (w,h) (bw,bh) (State scrRef keyFRef bxRef bpRef bvRef missRef arr) = do
  let org = (w/2,h/2)

  scr <- readIORef scrRef
  bx <- readIORef bxRef

  when (scr /= Prepare) $ do
    (leftF, rightF) <- readIORef keyFRef
    when (leftF && bx > -w/2) $ modifyIORef bxRef (subtract 5)
    when (rightF && bx+fi bw-10 < w/2) $ modifyIORef bxRef (+5)

    bv <- readIORef bvRef
    modifyIORef bpRef ((+fst bv)***(+snd bv))

  (bpx, bpy) <- readIORef bpRef
  miss <- readIORef missRef
  blocks <- readIORef arr 

  render cvs $ translate org $ do
    color (RGB 40 40 40) $ fill $ rect (bx,h/2-50) (bx+fi bw-10,h/2-50+fi bh-5)
    color (RGB 40 40 40) $ fill $ circle (bpx,bpy) 10
    color (RGB 40 40 40) $ font "30px ヒラギノ角ゴ" $ text (-w/2+10,-h/2+38) $ "Miss: " ++ show miss
    forM_ [0..bnh-1] $ \i -> do
      forM_ [0..bnw-1] $ \j -> do
        let alpha = case blocks!(i,j) of
                      NoDamage -> 1.0
                      Hit -> 0.2
                      Destroy -> 0

        opacity alpha $ do
          color (hsv ((fi bnh-1-fi i+fi j)*360/fi bnh) 180 240) $ fill $ rect ((w-fi bnw*fi bw)/2+fi bw*fi j+5-w/2, 60+fi i*fi bh-h/2) ((w-fi bnw*fi bw)/2+fi bw*fi j+5-w/2+fi bw-10, 60+fi i*fi bh+fi bh-5-h/2)
          let (rx,ry) = ((w-fi bnw*fi bw)/2+fi bw*fi j+5-w/2+2, 60+fi i*fi bh-h/2+2)
          color (hsv ((fi bnh-1-fi i+fi j)*360/fi bnh) 120 240) $ lineWidth 5 $ stroke $ path [(rx,ry), (rx+fi bw-10-4,ry), (rx+fi bw-10-4,ry+fi bh-5-4), (rx,ry+fi bh-5-4), (rx,ry-2.5)]

    when (scr /= Game) $ do
      opacity 0.8 $ color (RGB 40 40 40) $ fill $ rect (-w/2,-h/2) (w/2,h/2)
      color (RGB 255 255 255) $ font "60px ヒラギノ角ゴ" $ do
        when (scr == Prepare) $ text (-450/2,0) "Enter To Start!"
        when (scr == Clear) $ text (-240/2,0) "Cleared!"

  when (bpx < -w/2+10 || bpx > w/2-10) $ modifyIORef bvRef (&_1%~negate)
  when (bpy < -h/2+10 || bpy > h/2-10) $ modifyIORef bvRef (&_2%~negate) 
  when (bpy > h/2-10) $ modifyIORef missRef (+1)
  reflect (bw,bh) (bx,h/2-50) (bpx,bpy) bpRef bvRef

  forM_ [0..bnh-1] $ \i -> do
    forM_ [0..bnw-1] $ \j -> when (blocks!(i,j) /= Destroy) $ do
      refF <- reflect (bw,bh) ((w-fi bnw*fi bw)/2+fi bw*fi j+5-w/2, 60+fi i*fi bh-h/2) (bpx,bpy) bpRef bvRef

      when refF $ writeIORef arr $ blocks // [((i,j), succ $ blocks!(i,j))]

  when (all (==Destroy) $ elems blocks) $ writeIORef scrRef Clear

  setTimer (Once 10) (loop cvs (w,h) (bw,bh) (State scrRef keyFRef bxRef bpRef bvRef missRef arr))
  return ()

reflect :: (Int, Int) -> Point -> Point -> IORef Point -> IORef Point -> IO Bool
reflect (bw,bh) (bx,by) (x,y) bpRef bvRef = do
  if (y+10 >= by && y-10 <= by+fi bh-5 && x+10 >= bx && x-10 <= bx+fi bw-10)
    then do
      modifyIORef bvRef (&_2%~negate)
      when (y+10 < by+fi bh-5) $ modifyIORef bpRef (&_2%~(subtract 1))
      when (y-10 > by) $ modifyIORef bpRef (&_2%~(+1))
      when (y+10 < by || y-10 > by+fi bh-5) $ modifyIORef bvRef (&_1%~negate)
      return True
    else return False

fi :: (Integral a, Num b) => a -> b
fi = fromIntegral

hsv :: Double -> Double -> Double -> Color
hsv h s v = case hi of
  0 -> fRGB v k m
  1 -> fRGB n v m
  2 -> fRGB m v k
  3 -> fRGB m n v
  4 -> fRGB k m v
  5 -> fRGB v m n
  where
    h' = h-(fromIntegral $ floor $ h/360)*360
    hi = fromIntegral $ floor $ h' / 60
    f = h' / 60 - hi
    m = v * (1-s/255)
    n = v * (1-s/255*f)
    k = v * (1-s/255*(1-f))
    fRGB :: Double -> Double -> Double -> Color
    fRGB r g b = RGB (floor r) (floor g) (floor b)

練習あるのみと思って下手くそなりに書いてみたわけだけど,絶対もっとうまく書けるよね… IORef 使う以外にないのかなぁ…
Haste サンプルが少なくて使い方よくわかってないけど,とりあえず描画とキーの取得さえできればゲームはできる (?)

f:id:Ysmr_Ry:20170119094626g:plain

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]] となる.