年: 2016年

  • 素数夜曲とプログラミング B.2-B.4 Scheme, Python, Haskell, JavaScript 第 2 回

    素数夜曲とこれまでの内容

    素数夜曲

    これまでの内容

    B.2 名前と手続き

    リストの一番左は特別: そこには手続きを置く.
    命名と抽象化について割といいことが書いてある気がする.

    太字で強調してあるところがあって,
    何となく大事そうだからこちらにも書かせてもらおう.

    プログラミングは命令的知識を扱い, 数学は宣言的知識を扱う.

    命令的というのは「如何にして為すか」という意味で,
    宣言的というのは「何であるか」という意味.

    B.2.1 関数の定義

    素数夜曲では「函数」と書いているがめんどいので関数と書くことにする.
    いわゆるラムダの話から入る.
    次の関数は $f(x) = x^2$ だ.
    ラムダで書くと $\lambda x.x^2$.

    ラムダ以外の書き方もあるがとりあえず本の順に沿って書いていこう.

    
    (lambda (x) (* x x))
    (write ((lambda (x) (* x x)) 5))
    

    RESULT

    25
    

    一次関数 $x \mapsto ax + b$ は Scheme だと次のように書く.

    
    (lambda (a) (lambda (b) (lambda (x) (+ (* a x) b))))
    (write ((((lambda (a) (lambda (b) (lambda (x)
                                        (+ (* a x) b))))
              2) 3) 5))
    

    RESULT

    13
    

    こんな鬱陶しいのは嫌だ.
    LISP が嫌いになる理由として括弧が挙げられる理由を強く感じる.
    慣れるとむしろ便利なくらいと聞くが道は遠い.

    何はともあれもう少し簡単な書き方がラムダレベルでもあるし Scheme にもある.
    $\lambda abx.ax+b$ で, Scheme は次の通り.

    
    ((lambda (a b x)
       (+ (* a x) b))
     2 3 5)
    
    (write ((lambda (a b x)
              (+ (* a x) b))
            2 3 5))
    

    RESULT

    13
    

    素数夜曲でいつ出てくるのかわからないが,
    大分長いこと lambda をひっぱるようだ.
    鬱陶しいので適当に検索してきて lambda 抜きの関数定義を書いておく.

    
    (define (hello name)
      (string-append "Hello " name "!"))
    (write (hello "World"))
    

    RESULT

    "Hello World!"
    

    Python の lambda

    $f(x) = x^2$ を書いてみよう.

    
    myfunc = lambda x: x ** 2
    print(myfunc(2))
    

    RESULT

    4
    

    普通の関数で書いてみよう.

    
    def f(x):
        return x**2
    print(f(2))
    

    RESULT

    4
    

    Haskell の lambda

    Haskell での lambda は次の通り.

    
    mysquareInt = \x -> x ^ 2
    mysquareDouble = \x -> x ** 2
    
    main = do
      print $ mysquareInt 2
      print $ mysquareDouble 2.0
    

    RESULT

    4
    4.0
    

    JavaScript の lambda

    一般的な lambda と無名関数の違いがよくわかっていない.

    
    var f = function(a) { return Math.pow(a, 2) };
    console.log(f(3));
    var a = (function(a){ return Math.pow(a, 2) })(3);
    console.log(a);
    console.log((function(a){ return Math.pow(a, 2) })(3));
    

    RESULT

    9
    9
    9
    

    B.2.2 アルファ変換

    いわゆる束縛変数の文字を変えるというやつ.

    B.3 ラムダ算法

    Lambda calculus だし,
    たぶんいわゆるラムダ計算なのだろう.
    プログラムじたいはないのでばっさり省略.

    B.3.3 イータ変換

    省略.

    B.3.4 ラムダ項の定義

    略.

    B.3.5 コンビネータ

    略.

    B.3.6 簡約の戦略

    略.

    B.4 特殊形式 (special form)

    対象に名前をつけてプログラム上のどこからでも使えるようにすることを,
    トップレベル定義 (top-level definition) という.

    1 次関数を定義する.

    
    (define linear
      (lambda (a b x)
        (+ (* a x) b)))
    (write linear)
    (newline)
    (write (linear 2 3 5))
    

    RESULT

    #<procedure linear (a b x)>
    13
    

    最後の (linear 2 3 5) は $2 \times 5 + 3$ を計算している.

    Lambda を使わない省略記法もあって MIT 記法 という.

    
    (define (linear a b x)
      (+ (* a x) b))
    (write linear)
    (newline)
    (write (linear 2 3 5))
    

    RESULT

    #<procedure linear (a b x)>
    13
    

    この手の省略記法は糖衣構文 (syntax sugar) と呼ばれる.

    トップレベルは対話的入力時のプロンプトに象徴される領域で,
    大域環境 (global environment) と呼び,
    この環境下で定義された変数を大域変数 (global variable) と呼ぶ.
    これと対になるのが局所環境 (local environment),
    局所変数 (local variable) だ.

    アルファ変換した linear# も書いておこう.

    
    (define (linear# u v w)
      (+ (* u w) v))
    (write linear#)
    (newline)
    (write (linear# 2 3 5))
    

    RESULT

    #<procedure #{linear#}# (u v w)>
    13
    

    アルファ変換は局所変数を別の文字 (列) に置き換えてもいいとか,
    とりあえずそういう雑な理解をしている.
    計算機科学的にもっと面倒なところまでカバーしている可能性があり,
    そこまで調べていないからいまの理解の雑さは割と真剣に警戒している.

    B.4.1 define

    用語の定義があるから適当に書いておこう.
    対象を挿入する場所 (スロット) は本文と同じくを角括弧 (<>) で書きたいが,
    いろいろな事情から隅括弧 ([]) で書くことにする1.

    Lambda 抜きの糖衣構文を使った define の一般形は次の通り.

    
    定義方法: (define ([name] [fps]) [body])
    利用方法: ([name] [aps])
    

    本の記述が死ぬほどわかりづらいが,
    とりあえず引用しよう.

    手続の作用の対象となるのが引数 (parameter, argument) だ.
    [fps] は関数内部で使われる仮引数 (formal parameters) の略記で,
    [aps] は実際の処理対象である実引数 (actual parameters) の略記.
    [body] は任意個数の (expression) で構成される手続の本体だ.
    Scheme での式はアトムやリテラルを含めて評価値が戻ってくるもの全てを指す.

    とりあえずここまでくり返し書いてきた関数定義を.

    
    (define (linear a b x)
      (+ (* a x) b))
    

    linear が [name], [fps] は a, b, x で,
    [body] が (+ (* a x) b) だ.
    [name] は要は関数の名前,
    [fps] は定義した関数の引数,
    [body] は関数で実際にやる処理のこと.

    よくプログラミングのマニュアルでは上のタイプのよくわからない説明書きがある.
    こういうやつ.

    
    array explode ( string $delimiter , string $string [, int $limit = PHP_INT_MAX ] )
    

    こういうのを見たらとりあえず実際にコード片を書いて実行して確認してみた方がいい.
    めちゃくちゃ引数がたくさんあってそんなことやっていられないこともよくあるけれども.

    Scheme の一般評価規則にしたがわない対象,
    つまり特定の対象の評価値を求めないものを特殊形式 (special form),
    あるは構文 (syntax) と呼ぶ.

    スペシャルフォームと片仮名書きされているのをよく見かける.
    ここからは代表的な特殊形式を紹介していくようだ.

    一応改めて Python や Haskell での関数定義を書いておこう.

    Python での変数・関数定義

    まずは変数定義.

    
    a = 1
    print(a)
    

    RESULT

    1
    

    数を 2 乗する関数を定義する.

    
    def f(x):
        return x ** 2
    
    print(f(2))
    print(f(3))
    

    RESULT

    4
    9
    

    Haskell での変数・関数定義

    私が理解している限り, いわゆる変数はない.
    定数関数があってそれが定数の役割をしてくれるという理解.

    
    two :: Int
    two = 2
    
    main = do
      print two
    

    RESULT

    2
    

    こちらは整数を 2 乗する関数を定義する.

    
    mySquare :: Int -> Int
    mySquare x = x ^ 2
    
    main = do
      print $ mySquare 2
      print $ mySquare 3
    

    RESULT

    4
    9
    

    JavaScript の変数・関数定義

    use strict がないバージョン.
    グローバル変数になるのでよくないと言われるやつ.

    
    a = 2;
    console.log(a);
    

    RESULT

    2
    

    use strict はどこまで使えるようになっているのだろうか.
    use strict をつけると上のキーワードなしの宣言でエラーになる.
    次のキーワードはスコープとかの問題で詳しくはこことかで適当に調べる.

    
    "use strict";
    var a0 = 1;
    console.log(a0);
    var a = 2;
    console.log(a);
    let b = 3;
    console.log(b);
    const c = 4;
    console.log(c);
    

    RESULT

    1
    2
    3
    4
    

    略号まとめ

    • proc: procedure
    • exp: expression
    • num: number
    • predi: predicate
    • func: function
    • var: variable
    • int: integer
    • consq: consequent
    • args: argument
    • val: value
    • str: string
    • altna: alternative
    • char: character
    • obj: object
    • lst: list
    • stm: stream

    引数の有効範囲

    変数定義的な意味での define の使い方.
    正確には次のように言うべきのようだ.

    記憶領域のある場所に a という名前を与え, その場所と値 3.14 を結びつける.

    値を参照するラベルが a とかいうやつだろう.

    
    (define a 3.14)
    (write a)
    

    RESULT

    3.14
    

    define の次の用法の例でもある模様.

    
    定義方法: (define [name] [exp])
    

    次のように読めばいいようだ.

    [name] とは [exp] の名前である.

    束縛と環境

    a を定義したあとなら a を入力してもエラーが出ない.
    REPL 的なアレで挙動確認しているわけではないので,
    次のコード片では write を使っている.

    
    (define a 3.14)
    (write a)
    

    RESULT

    3.14
    

    このように変数に束縛された値を確認する方法を変数参照 (variable reference) という.
    次の実行結果を見ると a と定数 3.14 が関連づけられていることがわかる.

    
    (define a 3.14)
    (write a)
    (newline)
    (write (* 2 a))
    

    RESULT

    3.14
    6.28
    

    この a に値を格納するのを束縛 (bind) と呼ぶ.
    この束縛情報の全体を環境 (environment) と呼ぶ.

    B.4.2 lambda

    
    定義方法: (lambda ([fps]) [body])
    利用方法: ((lambda ([fps]) [body]) [aps])
    

    改めてコード例を書いておこう.
    数を 2 乗する関数を定義する.

    
    (lambda (x) (* x x))
    (write ((lambda (x) (* x x)) 2))
    

    RESULT

    4
    

    関数定義の実験として前置記法 (prefix notation) を中置記法後置記法で書いている.
    前置記法について prefix を作っている.

    
    (define prefix
      (lambda (proc a b)
        (proc a b)))
    (write (prefix + 2 3))
    (newline)
    (write (prefix * 2 3))
    

    RESULT

    5
    6
    

    これを中置記法 (infix notation), 後置記法 (postfix notation) で書く.

    
    (define infix
      (lambda (a proc b)
        (proc a b)))
    (write (infix 2 + 3))
    (newline)
    

    RESULT

    5
    
    
    (define postfix
      (lambda (a b proc)
        (proc a b)))
    (write (postfix 2 3 +))
    

    RESULT

    5
    

    Python, Haskell, JavaScript のコード例はいらないだろう.

    ブロック構造

    手続の評価値は内部の最後の式で決まる.

    
    (define (linear x)
      (define a 2)
      (define b 3)
      (+ (* a x) b))
    (write (linear 5))
    (newline)
    

    RESULT

    13
    

    次の式では最後の (* a b) の結果が返る.

    
    (define (linear1 x)
      (define a 2)
      (define b 3)
      (+ (* a x) b)
      (* a b))
    (write (linear1 3))
    

    RESULT

    6
    

    関数内部で定義した変数は外で参照できない.
    a は関数の中で定義されているだけでトップレベルで定義されていないから,
    関数の外で呼ぼうとするとエラーになる.

    次のコードを org-babel で実行するとエラーになる.
    よくわからないが guile -s で実行した結果を貼っておこう.

    
    (define a 3.14)
    (write a)
    (newline)
    
    (define (linear2 x)
      (define a 2)
      (define b 3)
      (+ (* a x) b)
      (* a b))
    (write (linear2 3))
    (newline)
    
    (write a)
    

    RESULT

    An error occurred.
    

    guile -s での実行結果

    3.14
    6
    3.14
    
    その他の言語

    同じようなスコープの問題がある.
    JavaScript だとグローバル変数に関する有名な問題がある.
    昔は必ず var を使えというのがあった.
    最近は var よりも use strictlet, const を使うのがいいのだろうか.
    識者のご意見求む.

    begin

    lambda の持つ列挙機能を含む特殊形式が begin.
    式を対象にする限りは同じらしい.
    begin だとトップレベル定義になる.

    
    (begin
      (define a 2)
      (define b 3))
    (write a)
    (newline)
    (write b)
    

    上のコード, org-babel で実行できない.
    guile -s での実行結果を書いておこう.

    guile -s での実行結果

    2
    3
    
    その他の言語

    Python, Haskell, JavaScript で似たようなのがあるだろうか?

    今回はこの辺で終わろう.


    1. 実体参照にすればいいのだが,
      今の私の腕ではスロットだけを自動で変換するスクリプトを作れないので諦めた.
      Markdown の引用の「>」や TeX 中の不等式を無視して,
      スロットだけを狙い打って実体参照に書き換えるコードを書く腕がない. 
  • 素数夜曲とプログラミング B-B.1.2 Scheme, Python, Haskell, JavaScript 第 1 回

    はじめに: なぜ素数夜曲か

    数学・物理学習とプログラミングと絡められないかと試行錯誤している.

    数学からの取っつきやすさ,
    物理からの取っつきやすさ,
    プログラミングからの取っつきやすさをそれぞれ考えないといけない.
    ある程度の体系性もほしい.

    数学や物理としても面白い内容にしたいし,
    プログラミングから見ても意味がある内容にしたい.

    どうしたものかとずっと思っていたのだが,
    素数という数学のキラーコンテンツとプログラミングを絡めた素数夜曲があった.

    微妙なところも多いとは思いつつ,
    メインエディタとして Emacs を使っているし,
    LISP 系の言語はきちんとやりたいとも思っている.

    そんなところでいい感じの落とし所という気がしたので,
    少しずつ読み進めつつコードの記録をしていこうと決めた.

    数学・物理の数値計算との相性の良さ,
    私の勉強してみたさとも合わせて Python や Haskell のコードも書いていけたらいいな,
    と思っている.
    あと何となく JavaScript もやってみよう.
    JavaScript は動きが速すぎて何かやってもすぐに動かなくなりそうで嫌なのだが,
    ここでやるくらいのことならそう簡単に陳腐化しないだろうと勝手に信じてやっていこう.

    基本的には付録のプログラミングパート,
    特に付録 B から記録をつけていく予定だ.

    Emacs の org-babel で書いているので,
    素数夜曲本編とは見た目ちょっと違うコードを書いていくかもしれない.
    Scheme, org-babel ともにあまりよくわかっていないがとりあえず進める.

    org-babel の Scheme が標準で Guile を使っていて,
    その変更の仕方がわからなかったのでとりあえず Guile を使う.
    これまで入れていた Gauche を使いたかったが org-babel が対応していないようだ.
    そうでないなら MIT/GNU Scheme が良かった気もするがこれも設定がわからない.

    でははじめよう.

    出力用関数

    
    (write "write test")
    (newline)
    (display "display test")
    

    RESULT

    "write test"
    display test
    

    B.1 評価値を得る

    とりあえず Hello, World! で.

    
    (write "Hello, Scheme!")
    

    RESULT

    "Hello, Scheme!"
    

    B.1.2 四則計算

    基本的な書き方

    加減乗除は次の通り.
    前置記法なのがポイント.
    割り算で分数にしてくれるのは割とポイント高い気がする.

    
    (write (+ 5 3))
    (newline)
    (write (- 5 3))
    (newline)
    (write (* 5 3))
    (newline)
    (write (/ 5 3))
    

    RESULT

    8
    2
    15
    5/3
    

    Python

    
    print(5 + 3)
    print(-5 + 3)
    print(5 * 3)
    print(5 / 3)
    

    RESULT

    8
    -2
    15
    1.6666666666666667
    

    Haskell

    まだ (?) org-babel との相性がよくないっぽい.
    ここ を参考にちょっと改造.

    Haskell でのあまりの計算には mod を使う.

    
    main = do
      print $ 5 + 3
      print $ -5 + 3
      print $ 5 + 3
      print $ 5 / 3
      print $ mod 5 3
      print $ 5 `mod` 3
    

    RESULT

    8
    -2
    8
    1.6666666666666667
    2
    2
    

    JavaScript

    
    console.log(5 + 3);
    console.log((-5 + 3));
    console.log(5 * 3);
    console.log(5 / 3);
    

    RESULT

    8
    -2
    15
    1.6666666666666667
    

    長めの計算をどう書くか

    例えば次の式.

    \begin{align}
    2 \div 3 + 5 \times 7 – 11
    \end{align}

    Scheme (LISP 系言語) だと括弧で細かく区切っていくから,
    演算の優先度みたいな面倒なことをあまり考えなくても済む.

    
    (define a (- (+ (/ 2 3) (* 5 7)) 11))
    (write a)
    

    RESULT

    74/3
    

    Python

    Python だと普通に書く.

    
    a = 2 / 3 + 5 * 7 - 11
    print(a)
    

    RESULT

    24.666666666666664
    

    Haskell

    Python っぽくも Scheme っぽくも書ける.
    演算子を括弧でくくると Scheme のように演算子を前に置けるから.

    
    main = do
      print $ 2 / 3 + 5 * 7 - 11
      print $ (-) ((+) ((/) 2 3) ((*) 5 7)) 11
    

    RESULT

    24.666666666666664
    24.666666666666664
    

    JavaScript

    
    console.log(2 / 3 + 5 * 7 - 11);
    

    RESULT

    24.666666666666664
    

    平方根を取る関数 (演算子?)

    素数夜曲 P.389 によると Scheme で sqrt は単項演算子らしい.

    
    (write (sqrt 2))
    

    RESULT

    1.4142135623730951
    

    入れ子にすれば複数回適用できる.
    Scheme は関数合成あるのだろうか.
    ここを見ると Gauche ならあるようだが.

    
    (write (sqrt (sqrt 4)))
    

    RESULT

    1.4142135623730951
    

    Python

    Python 版は次の通り.
    Python は math をインポートしないといけないのがめんどい.
    Python 3.5.1 を使っているが関数合成はないのだろうか.
    ちょっと調べたところでは標準ではなさそうだった.

    
    import math
    print(math.sqrt(2))
    print(math.sqrt(math.sqrt(4)))
    

    RESULT

    1.4142135623730951
    1.4142135623730951
    

    Haskell

    Haskell は関数合成 (.) がある.

    
    main = do
      print $ sqrt 2
      print $ sqrt $ sqrt 4
      print $ (sqrt . sqrt) 4
    

    RESULT

    1.4142135623730951
    1.4142135623730951
    1.4142135623730951
    

    JavaScript

    
    console.log(Math.sqrt(2));
    console.log(Math.sqrt(Math.sqrt(4)));
    

    RESULT

    1.4142135623730951
    1.4142135623730951
    

    自然対数の底, ネイピア数

    Scheme では exp を使えば出せる.
    $e$ 自体は定義されていない?

    
    (write (exp 1))
    

    RESULT

    2.718281828459045
    

    Python

    Python の場合はやはり math モジュールにある.

    
    import math
    print(math.e)
    

    RESULT

    2.718281828459045
    

    NumPy にもある.

    
    import numpy as np
    print(np.e)
    

    RESULT

    2.718281828459045
    

    Haskell 版

    $e$ 自体は定義されていないので Scheme と同じく指数関数から近似値を出す.

    
    main = do
      print $ exp 1
    

    RESULT

    2.718281828459045
    

    JavaScript 版

    Math.exp() が自然対数の底による指数関数で,
    Math.pow() は一般の底と指数の関数.

    
    console.log(Math.exp(1));
    console.log(Math.pow(2, 3));
    

    RESULT

    2.718281828459045
    8
    

    2 項演算子をいくつか

    まずは等号 = を.

    
    (write (= 1 1))
    (newline)
    (write (= 1 2))
    

    RESULT

    #t
    #f
    

    ここで #t は true, #f は false を意味している.
    == ではないのかとちょっと驚いた.

    ちなみに引数がたくさんある場合にも適用できる.

    
    (write (= 1 2 3 4 5 6))
    (newline)
    (write (= 1 1 1))
    (newline)
    (write (= 1 1 1 2))
    

    RESULT

    #f
    #t
    #f
    

    不等号 > >= などもある.
    不等号はそれぞれ隣の項に対して適用させていった結果を出力するようだ.

    
    (write (< 1 2 2 3))
    (newline)
    (write (< 1 2 3 4))
    

    RESULT

    #f
    #t
    

    Python

    等しくないことの判定は != だ.

    
    print(1 == 2)
    print(1 == 1)
    print(1 != 1)
    print(1 < 2)
    

    RESULT

    False
    True
    False
    True
    

    こちらは = だと代入になってしまうので == と重ねる.

    Haskell

    Scheme の (= 1 2 3 4) みたいなのはどう書けばいいのだろう.
    foldr で書けると思ったので foldr (==) 1 [1, 2, 3, 4] と書いてみたら怒られた.
    相談したら「それ型エラー起こしてるから. ちゃんとして」と言われた.
    Haskell, ちょっとしたところですぐこけるのでとてもつらい.

    結論からいうと適当に 2 つリストを作り,
    その要素を比較していく形で対処するのが普通らしい.
    他のところでもそういう処理をするのがいいというコメントを頂いたことがあるので,
    道具箱におさめておこう.

    
    main = do
      print $ 1 == 1
      print $ 1 /= 1
      print $ 1 < 2
      let xs = [1, 2, 3, 4]
      print $ and $ zipWith (==) xs (tail xs)
      let ys = [1, 2, 3, 4]
      print $ and $ zipWith (<) ys (tail ys)
    --  print $ foldr (==) 1 [1, 2, 3, 4] -- エラーになった
    

    RESULT

    True
    False
    True
    False
    True
    

    let を使うのが癪なので次回やる予定の lambda を使って書いてみる.

    
    main = do
      print $ (\xs -> and $ zipWith (==) xs (tail xs)) [1, 2, 3, 4]
      print $ (\xs -> and $ zipWith (<) xs (tail xs)) [1, 2, 3, 4]
    

    RESULT

    False
    True
    

    見づらい気もするので無理せず関数定義した方がいい気もする.

    
    mycompare :: [Int] -> Bool
    mycompare xs = and $ zipWith (<) xs (tail xs)
    
    main = do
      print $ mycompare [1, 2, 3, 4]
    

    RESULT

    True
    

    Haskell 追記

    またコメントを頂いてしまった.
    これ.

    
    and $ ap (zipWith (<)) tail [1..4]
    True
    
    and $ ap (zipWith (<)) tail [1,2,2,3]
    False
    

    実行してみよう.
    しばらくはまったが何となく ap を使うためにインポートが必要っぽい.

    
    import Control.Monad
    
    main = do
      print $ and $ ap (zipWith (<)) tail [1..4]
      print $ and $ ap (zipWith (<)) tail [1,2,2,3]
    

    RESULT

    True
    False
    

    最高か.

    そして「素人はこんなこともわからない」というのを思い出せさてくれるし,
    とてもいい体験だった.

    ちょっと調べはしたが ap が何なのかあまりよくわかっていない.
    とりあえずこれ.

    
    myap f g = \x -> f x (g x)
    main = do
      print $ and $ myap (zipWith (<)) tail [1..4]
      print $ and $ myap (zipWith (<)) tail [1,2,2,3]
    

    RESULT

    True
    False
    

    stackoverflow もきちんと全部読んだわけではないし,
    この雑な理解でいいとも思わないがとりあえずこれで掴まえられるところもある.
    「適当にしか理解していない (そして半端な理解だと割とどうでもいいところで無駄にはまる)」ことさえ理解しておけば,
    あとでいくらでも調節は効くと数学でいつも経験することをいつも通り信じて追記を終える.

    JavaScript

    比較には ===== がある.
    基本的に === でやるのが安全.

    
    console.log(1 === 2);
    console.log(1 === 1);
    console.log(1 !== 1);
    console.log(1 < 2);
    

    RESULT

    false
    true
    false
    true
    

    JavaScript で複数の比較はどう書くといいだろうか.
    reduce を使おうと思ったが,
    それは Haskell の最初の foldr でエラーになるのと同じ理由で駄目だ.
    ここを見ると zip がないようだ.
    ふつうにループを回すのでは面白くない.

    大文字小文字の区別

    ~~素数夜曲 P.390 の記述によると,~~
    ~~Scheme は原則として大文字小文字を区別しないらしいが実装により異なるとのこと.~~
    いわゆる処理系で変わるとかいうやつか?

    安全性を考えるなら小文字に統一するのがいいらしい.
    区別しないことを前提にキャメルケースを使うのも一手とか何とか.

    今回はこのくらいにしよう.

    Scheme 追記

    下にあるようにコメントを頂いた.

    Scheme の仕様は改定を重ねられているので、どの版であるかによって挙動が異なる部分があります。 Scheme の仕様はその正式名称を略して RnRS と呼ばれていて、 n の部分に何回目の改定であるかの番号が入ります。 最新は R7RS です。
    「大文字小文字を区別しない」というのは R5RS までの仕様です。 (ただ、 R5RS 処理系を名乗る処理系でも区別する処理系や、切り替えられる処理系はあります。) R6RS 以降では区別するのがデフォルトです。
    有理数や無限長整数などは R6RS 以外ではオプショナルな仕様です。 (R6RS で必須になった後、 R7RS では再びオプショナルになりました。) 処理系によっては桁あふれに配慮が必要な場合もあります。 有理数や無限長整数は Guile や Gauche では問題なく使えますので気にする必要はありませんが、その他に色々なところで仕様では未規定の部分があるので処理系を乗り換える機会にはつまらないことで躓くかもしれません。

    思っていたよりも面倒そうな話だった.
    ただ「ハマる可能性があるかもしれない」
    「ハマるところはこの辺」というのが頭の片隅にあるかないかだけでも
    全く違うことはこれまでの経験でよくわかっている.
    こんなありがたいことはない.
    こういうコメントを頂けるのが記事にしておく醍醐味だ.

    地道に続けよう.

  • 数列とは何か? 微分方程式のシミュレーションの観点から/中高数学駆け込み寺 第 7 回

    こちらに PDF があります.
    サイトでは見づらい方, コンテンツを手元に置いておきたい方は
    ダウンロードしてご覧ください.

    数列も関数

    ベクトル, 関数と来て今回は数列です.
    簡単に復習すると,
    関数は数と数の割り当てルールのことで,
    数列も関数なのでした.

    等比数列

    いろんな数列があります.
    高校でもいろいろやります.
    今回の話で 1 番大事な等比数列です.
    長々とやっても仕方ないので軽く復習しときましょう.

    $(x_n)$ を数列とします.
    高校だと {$x_n$} と中括弧で書くと思いますが大学の数学的な都合で小括弧で書くことにします.
    等比数列がどんな数列だったか,
    つまりどんなルールで与えられているのかというとこんなルールです:
    \begin{align}
    x_{n+1}
    =
    \alpha x_{n}.
    \end{align}
    $x_{n}$ が 0 のときでも成り立つ形で書きました.
    $x_{n}$ が 0 ではないならもちろん $x_{n+1} / x_{n} = \alpha$ と書けます.
    どんな $n$ に対しても比が一定な数列だから等比数列というのでした.

    等比数列の一般項

    このルールを使って一般項を出します.
    $x_{n+1} = \alpha x_{n}$ で,
    $x_{n} = \alpha x_{n-1}$ で,
    $x_{n-1} = \alpha x_{n-2}$ で,
    と延々続きます.
    これが終わるのは $x_1 = \alpha x_{0}$ です.
    順々に代入していくと $x_{n} = x_{0} \alpha^{n}$ になりますね.

    この計算, はまる人ははまります.
    実際私も高校の頃, 時々わけがわからなくなることがありました.
    試験で慌てていてパニクったことを思い出します.
    それはそれとして,
    ここで細々とした計算を詳しく解説してると本筋を見失うので省略します.
    数列をどんなところでどう使うのか,
    それをちゃんと見てからきちんとやった方がよさそうですしね.

    微分方程式と等比数列

    これがどこで出てきたかというと微分方程式を近似した方程式で出てきたのでした.
    放射性物質の崩壊で出てきたのはこんなやつです.
    \begin{align}
    \frac{x_{n+1} – x_{n}}{h}
    =
    – x_{n}.
    \end{align}
    これを整理すると $x_{n+1} = (1 – h) x_{n}$ で,
    $1-h = \alpha$ とすればまさにさっき書いた等比数列ですね.
    さっき書いたように等比数列は指数関数で書けます.
    で, 元の微分方程式の解も実は指数関数で書けます.

    もう一度: 数列も関数

    ちょっと先走ったことも書きました.
    何はともあれ, 数列は関数で, 関数は数と数の割り当てルールのことで,
    その割り当てルールとして等比数列っていう割り当て方があるんだ,
    ということです.
    微分方程式もその割り当てルールと見なせるっていうのが今書いたことです.
    微分方程式と揃えた言い方として上の式は差分方程式と言うこともあります.

    漸化式

    ちなみに上で紹介した割り当てルールはもっと一般的にできます.
    これ, 要は次の項をその前の項で決まると言っているだけなので,
    $a_{n+1} = a_{n}^2$ みたいにしてもいいし,
    $a_{n+1} = a_{n} + a_{n-1}$ みたいに 2 つ前と関係あるようにしてもいいです.
    この一般的なやつがいわゆる漸化式です.

    2 つの数列に対する漸化式

    もちろんこれ以外にもいろいろな漸化式があります.
    あとで使うのでもう 1 つ漸化式を紹介します.
    そのために 2 つの数列 $(x_{n})$, $(v_{n})$ を用意します.
    で, 次のような漸化式を考えます.
    \begin{align}
    x_{n+1} – x_{n}
    =
    v_{n}, \quad
    v_{n+1} – v_{n}
    =
    – \alpha^2 x_{n}.
    \end{align}
    2 つの数列が絡まりあう漸化式です.
    高校でやる言葉を使うなら,
    互いに階差数列を取ってると思ってもいいかもしれません.

    単振動の漸化式

    これがどこから出てくるかというと次の微分方程式からです.
    \begin{align}
    \frac{d^2 x(t)}{dt^2}
    =
    – \alpha^{2} x(t).
    \end{align}

    この微分方程式も物理でよく出てくる方程式で,
    いわゆる単振動の運動方程式です.
    さっきの漸化式を出すには $v(t) = x'(t)$ としてもとの方程式も $v’ = – \alpha^2 x$ とすればいいです.
    これを差分近似すると次の式が出てきます.
    \begin{align}
    x_{n+1}
    =
    x_{n} + h v_{n}, \quad
    v_{n+1}
    =
    v_{n} – h \alpha^{2} x_{n}.
    \end{align}
    $h$ が入っていて係数が微妙に違いますが,
    あまり気にしないでください.

    この辺の計算プログラムも準備してあります.
    この記事後半のプログラムパートを確認してください.


    アンケートの回答をお願いします

    今回もアンケートがあります.
    改善につなげるためぜひ回答をお願いします.

    ではまた次回をお楽しみに!

    プログラミングパート

    1 階の線型常微分方程式

    復習でもう 1 回.

    一番単純でしかも実際に使われる微分方程式としてまずは 1 階の線型常微分方程式を考えよう.
    ちょっと不吉な例であるが放射性物質の崩壊の方程式を紹介する.
    導出をしたければちゃんと物理を勉強してもらう必要がある.
    ここでは物理は省略して数学に集中する.

    \begin{align}
    \frac{dx}{dt} = – c x.
    \end{align}
    厳密解は $x = C_0 e^{-ct}$ だ.
    初期値を設定すれば $C_0$ はそこから決まる.

    微分を単純に離散化すると次のようになる.

    \begin{align}
    \frac{x_{n+1} – x_{n}}{\Delta t}
    =
    -c x_{n}.
    \end{align}

    $\Delta t$ は $h$ と書くこともある.
    整理すると次の通り.

    \begin{align}
    x_{n+1}
    =
    x_{n} – c (\Delta t) x_{n}.
    \end{align}

    これに沿って計算したのがいわゆるオイラー法.
    次のセルではこれをコードに落としている.

    1 階の方程式のプログラム

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    
    def radioactive_euler(nt, init = 10):
        dt = 2 / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        x[0] = init
    
        for i in range(1, nt):
            x[i] = x[i-1] - c * dt * x[i-1]
    
        # ベクトル計算で書き直したい
    
        return x
    
     # 近似解
    c = 1
    nt = 101
    init = 5
    x1 = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x1)
    
     # 厳密解
    t = np.linspace(0, 2, nt)
    x2 = init * np.exp(- c * t)
    plt.plot(t, x2)
    
     # 凡例
    plt.legend(['approximation', 'exact'])
    plt.show()
     # 描画
    #plt.savefig(img_file)
    #b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
    #img_str = "" % (b64)
    #print(img_str)
    

    RESULT

    2 階の常微分方程式: 単振動

    こちらは単振動.

    次に 2 階の常微分方程式を紹介しよう.
    高校の物理で出てくるばねの振動(単振動)がまさにこの例だ.
    項を増やすと減衰振動になったり、外力をつけたりといろいろなケースがある.
    まずは一番単純な式を考えよう.

    \begin{align}
    \frac{d^2 x}{dt^2}
    =
    – \omega^2 x.
    \end{align}

    オイラー法やルンゲ-クッタ法は
    1 階の方程式に対する計算法なので直接は使えない.

    今は中間処理として $v = dx/dt$ を置いて計算すればいい.
    これは単なる数値計算の便法ではない.

    速度の意味もあるから,
    という表面的な理由ではなくもっと深く解析力学の文脈で物理としても大事な視点だ.
    もっといえばシンプレクティック計算法などもっといい計算法にも発展する.

    オイラー法

    とりあえずオイラー法で計算したい.
    まずは微分方程式自体を書き直す.

    \begin{align}
    \frac{dx}{dt}
    =
    v, \quad
    \frac{dv}{dt}
    =
    – \omega^2 x.
    \end{align}

    これをオイラー法で近似しよう.

    \begin{align}
    x_{n+1}
    =
    x_{n} + h v_{n}, \quad
    v_{n+1}
    =
    v_{n} – h \omega^{2} x_{n}.
    \end{align}

    オイラー法をコードに落とす.

    オイラー法のプログラム

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib import animation
    
    def harmonic_euler(nt, init = (5, 0)):
        dt = t_range / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        v = np.zeros(nt)
        x[0] = init[0]
        v[0] = init[1]
    
        for i in range(1, nt):
            x[i] = x[i-1] + dt * v[i-1]
            v[i] = v[i-1] - dt * (omega ** 2) * x[i-1]
    
        # ベクトル計算で書き直したい
    
        return (x, v)
    
    
    omega = 2 * np.pi
    nt = 101
    t_range = 2
    init = (5, 0)
    harm = harmonic_euler(nt, init)
    t = np.linspace(0, 2, nt)
    
     # 厳密解
    x_exact = init[0] * np.cos(- omega * t)
    v_exact = - omega * init[0] * np.sin(omega * t)
    
     # グラフ描画
    plt.subplot(3, 1, 1)
    plt.title('x-t graph')
    plt.plot(np.linspace(0, 2, nt), harm[0])
    plt.plot(t, x_exact)
    plt.legend(['x approximation', 'x exact'])
    
    plt.subplot(3, 1, 2)
    plt.title('v-t graph')
    plt.plot(np.linspace(0, 2, nt), harm[1])
    plt.plot(t, v_exact)
    plt.legend(['v approximation', 'v exact'])
    
    plt.subplot(3, 1, 3)
    plt.title('phase space')
    plt.plot(harm[0], harm[1])
    
     # 描画
    plt.tight_layout()
    plt.show()
    #plt.savefig(img_file)
    #b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
    #img_str = "" % (b64)
    #print(img_str)
    

    RESULT

    オイラー法では近似の精度が悪い

    見ての通り時間が進むごとに誤差が大きくなる.
    nt を大きくすると少しはましになる.
    実際に上のコードで nt を大きくして再計算してみてほしい.

    あとまずいのは phase space の図だ.
    この系はエネルギーが保存する系だから相空間内の軌道が閉じてほしいのにそうなっていない.
    シンプレクティックにやれば解消できるようだが, とにかくここではよろしくない.

    この方程式でオイラー法はよろしくないことがわかった.
    とりあえずルンゲ-クッタでやってみよう.

    ルンゲ-クッタ法

    とりあえず近似式を書く.

    \begin{align}
    x_{n+1}
    &=
    x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}),
    \end{align}
    \begin{align}
    t_{n+1}
    &=
    t_{n} + h,
    \end{align}
    \begin{align}
    k_{1}
    &=
    v_{n},
    \end{align}
    \begin{align}
    k_{2}
    &=
    v_{n} + \frac{h}{2} k_{1},
    \end{align}
    \begin{align}
    k_{3}
    &=
    v_{n} + \frac{h}{2} k_{2},
    \end{align}
    \begin{align}
    k_{4}
    &=
    v_{n} + h k_{3}.
    \end{align}

    次が $v$ の式.

    \begin{align}
    v_{n+1}
    &=
    v_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}),
    \end{align}
    \begin{align}
    t_{n+1}
    &=
    t_{n} + h,
    \end{align}
    \begin{align}
    k_{1}
    &=
    – \omega^2 x_{n},
    \end{align}
    \begin{align}
    k_{2}
    &=
    x_{n} – \frac{h}{2} \omega^2 k_{1},
    \end{align}
    \begin{align}
    k_{3}
    &=
    x_{n} – \frac{h}{2} \omega^2 k_{2},
    \end{align}
    \begin{align}
    k_{4}
    &=
    x_{n} – h \omega^2 k_{3}.
    \end{align}

    ルンゲ-クッタ法のプログラム

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib import animation
    
    def harmonic_rk(nt, init = 10):
        dt = t_range / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        v = np.zeros(nt)
        x[0] = init[0]
        v[0] = init[1]
    
        def fx(t, x, v):
            return v
    
        def fv(t, x, v):
            return - (omega ** 2) * x
    
        # ベクトル計算で書き直したい
        for i in range(1, nt):
            xk1 = fx(dt * (i - 1), x[i-1], v[i-1])
            vk1 = fv(dt * (i - 1), x[i-1], v[i-1])
    
            xk2 = fx(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)
            vk2 = fv(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)
    
            xk3 = fx(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)
            vk3 = fv(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)
    
            xk4 = fx(dt * (i - 1),   x[i-1] + xk3 * dt,     v[i-1] + vk3 * dt)
            vk4 = fv(dt * (i - 1),   x[i-1] + xk3 * dt,     v[i-1] + vk3 * dt)
    
            x[i] = x[i-1] + dt / 6 * (xk1 + 2 * xk2 + 2 * xk3 + xk4)
            v[i] = v[i-1] + dt / 6 * (vk1 + 2 * vk2 + 2 * vk3 + vk4)
    
        return (x, v)
    
    
    omega = 2 * np.pi
    nt = 101
    t_range = 2
    init = (5, 0)
    harm = harmonic_rk(nt, init)
    t = np.linspace(0, 2, nt)
    
    
    x_exact = init[0] * np.cos(- omega * t)
    v_exact = - omega * init[0] * np.sin(omega * t)
    
    
    plt.subplot(3, 1, 1)
    plt.title('x-t graph')
    plt.plot(np.linspace(0, 2, nt), harm[0])
    plt.plot(t, x_exact)
    plt.legend(['x approximation', 'x exact'])
    
    plt.subplot(3, 1, 2)
    plt.title('v-t graph')
    plt.plot(np.linspace(0, 2, nt), harm[1])
    plt.plot(t, v_exact)
    plt.legend(['v approximation', 'v exact'])
    
    plt.subplot(3, 1, 3)
    plt.title('phase space')
    plt.plot(harm[0], harm[1])
    
     # 描画
    plt.tight_layout()
    plt.show()
    #plt.savefig(img_file)
    #b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
    #img_str = "" % (b64)
    #print(img_str)
    

    RESULT

    ルンゲ-クッタ短評

    今度の一致具合はなかなかよさそう.
    あくまで見た目の感じではあるけれども.

    良くなったのは一番下の図, 相空間軌道だ: ちゃんと閉じてくれた.

  • 関数とは何か? 微分方程式のシミュレーションの観点から/中高数学駆け込み寺 第 6 回

    こちらに PDF があります.
    サイトでは見づらい方, コンテンツを手元に置いておきたい方は
    ダウンロードしてご覧ください.

    ベクトルの復習

    前回はベクトルを紹介しました.
    ベクトルは始点と終点をつなぐ矢印であり,
    その矢印を辿って折れ線を作って曲線を近似するんだ,
    そういう視点からの紹介です.

    今回は関数です.
    そして数列, 漸化式, 微分, 微分方程式とつないでいく予定です.

    関数は対応ルール

    で, 関数.
    とりあえず数と数を適当なルールで結びつける装置のことだと思ってください.
    関数のよくある説明は自動販売機ですね.
    自動販売機にお金を入れて何かボタンを押すと,
    押したボタンに対応した飲み物が出てきます.
    この対応がまさに関数です1.

    ここで数学として大事なのは押したボタンに応じていつでも同じ飲み物が出てくることです.
    同じボタンを押したのに違う飲み物が出てこられたらルール壊れてますからね.

    数学で具体例を出すなら $f(x) = x^2$.
    $x$ を入れたら必ず $x^2$ を出してくれる機械を関数と呼んでいるわけです.
    この関数を支配するルールはもちろん $x$ を $x^2$ に変換することです.

    微分方程式も対応ルール

    このルールは本当は何でもいいです.
    この講座で大事なのはルールが微分方程式で決まる場合です.
    かといっていきなり微分にかっとぶのも大変なのでもうちょっとクッションを入れます.

    あなたは関数は実数に対して実数をあてるルールと思っているかもしれません.
    実際にはもっとゆるくて何かと何かを適当に割り当てるルールくらいに思っておいてください.
    例えば自然数2に対して実数を割り当てる関数だってあります.
    これはふつう数列と呼ばれています.
    つまり数列も関数なんですね.
    あなたが数列を勉強したことがあるなら,
    数列は「ある規則にしたがって数を書き並べたもの」みたいな説明を受けなかったでしょうか?
    この説明, まさに関数と同じです!

    数列も対応ルール

    高校だと数列は特にルールっていう感じで出てきます.
    等差数列とか等比数列とか階差数列とかもそうだし,
    漸化式なんてまさに数列を作っていくルールのことですしね.

    最初にシミュレーションでやったことを思い出してもらうと,
    プログラムを書いて計算していったのはまさに漸化式です.

    ベクトルも対応ルール

    あとあなたが大学の数学やその他の学問に興味があるなら次のようなことまで知っておくと楽しいかもしれません:
    数列は自然数全体を適当なルールで実数に割り当てる関数でした.
    ここで自然数全体じゃなくて $1$ から $n$ までを実数に割り当てる関数を考えてみましょう.
    実はこれはベクトルです.
    $n$ 次元のベクトルってやつです.
    ベクトルは他にいろいろ要件があるのであんまり適当なこと言うのも微妙なんですが.

    実在する高次元ベクトル

    ちなみに「$3$ より高い次元のベクトルなんてどこに使うの?」みたいなこと,
    工学部の大学生が言っているのを聞いたことがあります.
    でも実はシミュレーションでこういう高次元のベクトルをバリバリ使っています.

    今回のポイント

    最後のベクトルの話はわからなくても問題ないです.
    関数は数と数の対応ルールであること,
    数列は特殊な関数であること,
    この 2 つを覚えておいてください.

    今回は概念の説明だけだったのでちょっと難しかったでしょう.
    次回の数列では高校の復習しながら具体的な話もします.
    今回の話がわからなくてもあまり気にしないでください.

    今回はここまでです.
    お疲れ様でした.

    アンケートをお願いします

    今回もアンケートがあります.
    改善につなげるためぜひ回答をお願いします.

    ではまた次回をお楽しみに!


    1. うるさいことを言うと,
      売り切れているときにその飲み物が出てこないので数学の関数とは正確には一致していません.
      他にも問題はあって,
      例えば現実問題としては時々全然違う飲み物が出てくることだってありますね?
      今はそういう細かいことは無視して理想的な自動販売機を考えてください.
      ちなみにこういううるさい細々としたことが気になるあなた,
      実はとても数学に向いていますよ.
      大学の数学だとこういう話を山ほどやります.
      大学の数学もちょっと紹介する予定なので楽しみにしててください.
      ちなみに完全によけいなコメントなんですが,
      自動販売機はオブジェクト指向プログラミングでのインスタンスメソッドと言った方が正確な気がします. 
    2. この講座ではふつうの中高の数学の教科書とは違って $0$ は自然数だとします.
      Python のリスト (配列) は $0$ からはじまるからプログラムを組むときにその方が都合がいいんです. 
  • ベクトルとは何か? 微分方程式のシミュレーションの観点から/中高数学駆け込み寺 第 5 回

    こちらに PDF があります.
    サイトでは見づらい方, コンテンツを手元に置いておきたい方は
    ダウンロードしてご覧ください.

    ここまでの復習

    ここまでは微分方程式が自然現象のシミュレーションに使えることを紹介し,
    実際の微分方程式をいくつか見てきました.
    実際のシミュレーションでは本当に微分方程式を厳密に解いているわけじゃなくて,
    近似計算をしているんだ,
    そしてそれは四則演算なんだというのを大雑把に紹介してきました.

    近似といってもけっこうきちんと計算できますよ,
    ということを見るためにもとにかくいくつかの実例を紹介することを優先にしてきました.
    我慢強いあなたも「もういい加減きちんと数学の紹介して!」
    と思ってらっしゃるはずです.

    ベクトルからはじめます

    というわけで数学の紹介も簡単にやっていきましょう.
    微分方程式なんだからまずは微分じゃないの?
    と思う方もいらっしゃるでしょう.
    でも今回の話の流れでは微分の話をする前にベクトルをやった方がスムーズなので,
    ベクトルからいきます.

    ベクトルは矢印

    高校でベクトルは向きと大きさがある矢印のことと言われます.
    それはそれで大事な見方ですが,
    ここでは始点と終点がある矢印と思ってみてください.
    始点と終点をつないだたくさんのベクトルで作った折れ線を方程式の近似解とみなすのが今回の話のキモだからです.

    図のように曲線を適当な間隔でわけましょう.

    ベクトルをつないだ折れ線で曲線を近似.
    ベクトルをつないだ折れ線で曲線を近似.

    わけた点に $A_0, A_1, \dots, A_n$ と順に名前をつけていきます.
    そして $A_0$ と $A_1$, $A_1$ と $A_2$ というふうに始点と終点を入れかえながら隣の点どうしを結んでいきます.
    $A_0$ から $A_1$ に向かって進むんだ,
    という気持ちを表したいのでそれを $\overrightarrow{A_0 A_1}$ と書くことにします.
    記号の上の右向き矢印で点 $A_0$ から点 $A_1$ に向かっている気分もはっきり書いています.

    ここで何度も始点と終点をつないだ矢印を継ぎ足しています.
    この継ぎ足しという図形の操作がベクトルの足し算です.
    例を 1 つだけ書いておくと $\overrightarrow{A_0 A_1} + \overrightarrow{A_1 A_2} = \overrightarrow{A_0 A_2}$ です.

    微分方程式のシミュレーションから見たベクトル

    高校だとベクトルに関してもっといろいろなことが出てきます.
    定数倍したり内積を取ったり長さを調べたりなどなど.
    もちろんベクトルを使っていろいろやるならやっておかないと困ります.
    でも微分方程式から数学を眺めてみる立場ではまずここさえ乗り越えればどうにかなります.
    他のことは他のことをやるときに絡めてやってみてください.

    折れ線をつないでいって曲線を近似する様子をもっときちんと見てみましょう.
    具体的な曲線としては円を取り,
    それを等分して折れ線近似した様子をプログラムで描きました.
    次のページを見てください.
    最後に近似がよくなっていくアニメーションもつけています.

    基本的なプログラムはこの記事の後半にも載せています。
    アニメーションは Jupyter との連携でやっているので,
    この記事ではうまく動かせません.
    上のリンクのページから確認してください.

    今回はここまでです.
    お疲れ様でした.

    今回もアンケートがあります.
    改善につなげるためぜひ回答をお願いします.

    ではまた次回をお楽しみに!

    プログラミングパート

    円を折れ線近似する関数

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    
    div_num = 9
    def draw_circle_by_polygonal_line(div_num):
        thetas = np.linspace(0, 2 * np.pi, div_num)
        xs = np.cos(thetas)
        ys = np.sin(thetas)
    
        plt.plot(xs, ys)
        plt.gca().set_aspect('equal', adjustable='box')
    

    ベクトルでつないだ折れ線で円を近似する

    円周上の点を $n$ 分割し折れ線で分点をつないでいった.
    一旦個別に図示している.

    
    for i in range(4, 12):
        pnum = i - 3
        plt.subplot(2, 4, pnum)
        plt.xlim([-1.2, 1.2])
        plt.ylim([-1.2, 1.2])
        plt.title("div = " + str(i-1))
        draw_circle_by_polygonal_line(i)
        draw_circle_by_polygonal_line(1000)
    
    plt.tight_layout()
    

    RESULT

    ベクトルをつないだ折れ線で円を近似する

    9 分割までを重ねて 1 つの図にしてみた.

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    
    for div_num in range(10):
        draw_circle_by_polygonal_line(div_num)
    
    draw_circle_by_polygonal_line(1000)
    plt.xlim([-1.2, 1.2])
    plt.ylim([-1.2, 1.2])
    plt.show()
    

    RESULT

    アニメーションにしてみる

    分割を細かくしていくことで精度が上がる様子をアニメーションで見てみよう.

    下のプログラムを Jupyter 上で実行すると,
    図の両端にある「+」「-」ボタンでアニメーションのスピードを変えられる.
    適当に変えて遊んでみよう.

    注意
    次のプログラムは Jupyter 上で実行することを前提にしている.
    他の実行環境で動くかどうかは確認していない.
    ここを見ると Jupyter のアニメーションが見られる.

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib import animation
    from JSAnimation.IPython_display import display_animation
    
    def circle_data(div_num):
        thetas = np.linspace(0, 2 * np.pi, div_num)
        xs = np.cos(thetas)
        ys = np.sin(thetas)
        return (xs, ys)
    
    
    circles = []
    for i in range(4, 31):
        circles.append(circle_data(i))
    
    
    #fig = plt.figure(figsize=(10, 10));
    fig = plt.figure();
    
    
    ax = plt.axes(xlim=(-1.2, 1.2), ylim=(-1.2, 1.2));
    line, = ax.plot([],[],lw=2);
    
    
    plt.gca().set_aspect('equal', adjustable='box')
    
    def animate(data):
        x = data[0]
        y = data[1]
        line.set_data(x, y)
        return line
    
    anim = animation.FuncAnimation(fig, animate, frames=circles, interval=100)
    display_animation(anim, default_mode='reflect')
    
  • 微分方程式に関するここまでのまとめ/中高数学駆け込み寺 第 4 回

    こちらに PDF があります.
    サイトでは見づらい方, コンテンツを手元に置いておきたい方は
    ダウンロードしてご覧ください.

    次回からもう少し数学的に詳しく踏み込んでいきます.
    その前にここまでの話を簡単にまとめておきましょう.

    微分方程式全般

    • 微分方程式というのがあって理工系の大学生は必ず勉強する.
    • 自然法則は微分方程式で表現されることが多い.
    • テレビでもよく見かける洪水のシミュレーションや震源予測でも微分方程式が使われている.
    • テレビゲームのリアルな CG でも使われる.
    • 数学と縁遠そうな生物や文系の経済でも使われる.

    中高数学とシミュレーション

    ここでは微分方程式をそのまま扱うよりも差分で近似してシミュレーションで遊ぶことを目的にしています.
    そしてそのシミュレーションための差分計算で中高の数学を総動員するのでした.
    いくつか箇条書きでまとめておきましょう.

    • 微分方程式は差分でよく近似できる.
    • 差分にすれば四則演算しか使わない.
    • 微分方程式を差分で近似すると数列の漸化式が出てくる.
    • シミュレーションで本当の解 (厳密解) をかなりよく近似できる.
    • 本当の解として指数関数や指数関数が入った分数関数が出てくる.

    後で出す微分方程式の厳密解やその解を導く間に三角関数や対数関数も出てきます.
    当然これの微分積分も必要です.
    シミュレーションの気持ちを知るためにはベクトルが役に立ちますし,
    連立 1 次方程式を解くのに行列を知っていると便利です.
    微分方程式の厳密解を出すのに複素数が使えると便利な場面も多いし,
    実際に大学だとよく使いますね.
    電気回路の理論をやると必ず出てきます.
    ちなみに理工系で実験やろうと思うと必ずどこかしらに電気回路が出てきますから,
    理工系なら複素数使えないとまずいというのも思い知らされます.

    複素数と行列は高校の課程に入ったり入らなかったりするので,
    あなたが高校生ならなおのことどちらかまたは両方を知らないかもしれません.
    この中高数学駆け込み寺でも,
    全体像を知ってモチベーションを上げてもらうという目的があるのであまり深入りはしません.
    興味がある方のために最後に参考文献も含めて今後の勉強の指針をお伝えします.
    それまで少しの間待っていてください.

    次回からは個別の数学

    次回からは数学に関して少しずつ踏み込んでいきます.
    メインの微分や微分方程式の前にベクトルからはじめます.
    シミュレーションを考えるときにも大事なんです.

    アンケートをお願いします.

    今回もアンケートがあります.
    ぜひ回答をお願いします.

    ではまた次回をお楽しみに!

  • 経済や生物で使う微分方程式/中高数学駆け込み寺 第 3 回

    こちらに PDF が置いてあります.
    サイトでは見づらい方は PDF をご覧ください.

    マルサスの『人口論』と微分方程式

    あなたが高校生以上なら社会の授業で次のような言葉を聞いたことがあるかもしれません.

    人口は制限されなければ幾何級数的に増加するが生活資源は算術級数的にしか増加しないので生活資源は必ず不足する.

    これは有名なマルサスの『人口論』で議論された話です.
    幾何級数はいわゆる等比数列の和で,
    算術級数は等差数列の和です.
    マルサスは経済学者なので経済の話です.
    経済の話でこういう数学ネタが出てくるわけですね.

    ちなみに大学の経済学だと経済の話で微分積分が本当に出てきます.
    ときどき経済学部の入試で数学が受験科目に入っていることがあります.
    実際に大学に入ってから使うからです.
    それも高校の理系のレベルを越えた数学です.
    統計データをいろいろいじらないといけないので統計学が必要で,
    そっちでも割といろいろ数学が必要です.

    話を元に戻しましょう.
    マルサスの人口論に合わせて人口の増え方を考えます.
    ここではもっと一般的に生物の集団だと思いましょう.
    まずは生物の種類が 1 種類だとして考えていきます.
    生物っぽく人の数というよりも一般に個体数と呼ぶことにして,
    時刻 $t$ での個体数を $x(t)$ とします1.
    個体数の増加率を出生率の死亡率の差として定義します.
    微分方程式的にいうとこれは $x'(t) / x(t)$ です.
    増加率が一定値 $\alpha$ に等しいとすると $x'(t) / x(t) = \alpha$ だから次の微分方程式が出てきます.

    \begin{align}
    \frac{d x(t)}{dt}
    =
    \alpha x(t).
    \end{align}

    この微分方程式をマルサスモデルと言います.
    これを解くと $x(t) = x(0) e^{\alpha t}$ です.
    これで指数関数が出てきました.
    等比数列は一般項が $a_n = 2^{n}$ のように指数関数で書けている数列だから,
    これを足し上げていけばまさに幾何級数になりますね.
    放射性物質の崩壊とは指数の肩の符号が変わっているだけで,
    それ以外は同じ式です.
    放射性物質の崩壊と生物の個体数の変化を追う方程式が同じ形をしているわけです.

    シミュレーション

    シミュレーションの結果は放射性物質の崩壊のときと基本的に同じです.
    いちおうきちんとやっておきましょうか.
    まず微分方程式を次のように近似します.

    \begin{align}\label{math_refuge00005}
    \frac{x_{n+1} – x_{n}}{h}
    =
    \alpha x_n, \quad
    x_{n+1}
    =
    \alpha x_n + h x_n.
    \end{align}

    これにしたがって解いてみましょう.
    プログラムと結果の図はこのページの後半を確認してください.


    少し現実的な設定にしてみる

    もう少し現実的な設定にしてみます.
    個体数 $x(t)$ が大きくなると食物が足りなかったり,
    衛生状態の悪化で病気にかかりやすくなったりして増加率が減るはずです.
    そこで増加率 $\alpha$ を $\alpha – \beta x$ ($\alpha$, $\beta > 0$) に変えてみます.
    \begin{align}
    \frac{dx}{dt}
    =
    (\alpha – \beta x) x.
    \end{align}
    これにはロジスティック方程式という名前がついています.
    悪化を $\beta x$ と書くことに必然性はありません.
    とりあえずやってみただけです.

    ロジスティック方程式の厳密解

    解き方はさておき解は次のようになります.
    あなたがもしこの関数を微分できるなら,
    厳密解になることを計算して確認してみてください.
    \begin{align}
    x(t)
    =
    \frac{\alpha}{\beta} \frac{1}{1 + e^{- \alpha t}}.
    \end{align}

    これを近似すると次のようになります.
    \begin{align}
    \frac{x_{n+1} – x_{n}}{h}
    =
    (\alpha – \beta x_{n}) x_{n}, \quad
    x_{n+1}
    =
    x_n + h (\alpha – \beta x_{n}) x_{n}.
    \end{align}

    これも厳密解と比較しながら計算してみます.
    プログラムや図はこの記事の後半を参照してください.

    近似の精度が気になります

    これもよく近似できていることは認めてもらえるんじゃないでしょうか?
    でもうるさい人は「近似の精度はどう決めるの?」みたいなことを考えているはずです.

    もちろんいたって真っ当な指摘です.
    ただけっこう難しい問題がたくさんあって,
    なかなかスカっと綺麗に回答できません.
    あとで少し考えることにしてここではこのまま進みます.

    生物学として正しいの?

    あなたは数値計算の精度とかそれ以前の問題として,
    次のように思っているかもしれません.

    • この微分方程式で実際の生物の増減にあてはめられるの?
    • もっと現実的なことを考える必要があるんじゃないか?

    もちろんこれも正しい指摘です.

    例えばマルサスモデルでもロジスティック方程式でも右辺が気になります.
    実際には成長しないと子作りできないですからね.
    その成長が必要だという情報が方程式に盛り込まれていません.

    生物で考えるなら食物連鎖があるわけで食べられて個体数が減ることだってあります.
    天敵が増えたらその個体の数はあおりを受けて激減するはずです.

    もちろんこんなことは折り込み済みで,
    例えば遅延型微分方程式,
    ロトカ-ボルテラ方程式と呼ばれる方程式が対応しています.
    これを調べるのも面白いんですが,
    今回はこのくらいにしておきましょう.

    今回のポイント

    簡単にまとめると,
    今回は近似の精度が気になるあなたのために,
    実際どのくらいよく本当の解を近似できているのかを調べました.
    経済学や生物のようにあんまり数学との関係がなさそうな分野と数学の関係も紹介しています.

    数学的ポイント: 数列と漸化式

    長くなってきましたが最後にちょっと大事な話.
    式 \eqref{math_refuge00005} で $\alpha = 1$, $h = 1$ とすると $x_{n+1} = 2 x_n$ になります.
    この漸化式は高校でも出てくるやつで, もちろん公比 $2$ の等比数列の漸化式です.
    底が $2$ か $e$ かの違いはあっても指数関数で書けることは同じです.
    これはたまたまではなくて数学的に意味があることです.

    $x_{n+1} – x_{n}$ のように数列の隣の項の差を差分と言うことがあります.
    もちろん微分との関係を強く意識した言葉です.
    数列の問題, もっと言えば漸化式の問題は微分方程式とも深い関係があります2.
    あまり深入りはしませんが数学や物理に興味があるなら覚えておくといろいろ楽しいことがあります.

    アンケート

    今回もアンケートがあります.
    改善に繋げたいのでぜひ積極的に回答をお願いします.

    ではまた次回をお楽しみに!

    プログラミングパート

    マルサスモデル

    基本は放射性物質の崩壊と同じだからあっさり行こう.

    \begin{align}
    \frac{dx}{dt}
    =
    \alpha u.
    \end{align}
    厳密解は $x(t) = C_0 e^{\alpha t}$ だ.
    初期値を設定すれば $C_0$ はそこから決まる.

    微分を単純に離散化すると次のようになる.

    \begin{align}
    \frac{x_{n+1} – x_{n}}{h}
    =
    \alpha x_{n}.
    \end{align}

    整理すると次の通り.

    \begin{align}
    x_{n+1}
    =
    x_{n} + \alpha h x_{n}.
    \end{align}

    オイラー法でコードに落とそう.

    マルサスモデルの近似解プログラム

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    
    def malthus_euler(nt, init = 10):
        dt = 2 / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        x[0] = init
    
        for i in range(1, nt):
            x[i] = x[i-1] + c * dt * x[i-1]
    
        return x
    
    
    c = 1
    init = 1
    nt = 101
    x_approx = malthus_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT

    放射性物質の崩壊の時との違い

    t が大きくなると厳密解の方が少し大きくなる.
    指数の肩の符号が違うだけなので,
    放射性物質の崩壊で形式的に時間を負の方に伸ばせば同じようにズレが出てくるはずだ.
    私の今の感覚だとこのズレが今回のように大きくなるように出るのかどうかまではわからない.

    何はともあれ区間の分割数 nt を大きくしてみたのが次の結果:
    具体的には 101 から 1001 にした.

    パラメータを修正して再計算

    ライブラリが読み込まれていたり関数 malthus_euler が定義されている前提のコード片なので,
    これだけで実行できるわけではない.
    IPython か Jupyter なら順にセルを読み込んでいけばそのまま実行できる.

    
    
    c = 1
    init = 1
    nt = 1001
    x_approx = malthus_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT

    注意

    $nt = 101$ よりも近似の精度が良くなった.
    この辺は単純に振る舞ってくれるようだ.

    調和振動子だとオイラー法自体がうまく動かないから注意したい.

    ロジスティック方程式

    方程式は次の通り.

    \begin{align}
    \frac{dx}{dt}
    =
    (\alpha – \beta x) x.
    \end{align}

    厳密解は次の通り.

    \begin{align}
    x(t)
    =
    \frac{\alpha}{\beta} \frac{1}{1 + e^{- \alpha t}}.
    \end{align}

    近似すると次の通り.

    \begin{align}
    \frac{x_{n+1} – x_{n}}{h}
    =
    (\alpha – \beta x_{n}) x_{n}, \quad
    x_{n+1}
    =
    x_n + h (\alpha – \beta x_{n}) x_{n}.
    \end{align}

    では数値的に解いてみよう.
    とりあえず $\alpha = 2$, $\beta = 1$ で計算している.

    ロジスティック方程式の近似解プログラム

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    
    def logistics_euler(nt, init = 10):
        dt = T / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        x[0] = init
    
        # ベクトル計算で一気に計算したい
        for i in range(1, nt):
            x[i] = x[i-1] + dt * (alpha - beta * x[i-1]) * x[i-1]
    
        return x
    
    
    alpha = 2
    beta = 1
    init = 1
    nt = 101
    T = 5 # 時間変化を見る最大値
    x_approx = logistics_euler(nt, init)
    plt.plot(np.linspace(0, T, nt), x_approx)
    
    
    t = np.linspace(0, T, nt)
    x_exact = (alpha / beta) / (1 + np.exp(- alpha * t))
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT

    時間の刻みを変えてみる

    近似解と厳密解にちょっとズレが見える.
    そこで nt = 101 から nt = 1001 にしてみた.

    パラメータを修正して再計算

    
    
    alpha = 2
    beta = 1
    init = 1
    nt = 1001
    T = 5 # 時間変化を見る最大値
    x_approx = logistics_euler(nt, init)
    plt.plot(np.linspace(0, T, nt), x_approx)
    
    
    t = np.linspace(0, T, nt)
    x_exact = (alpha / beta) / (1 + np.exp(- alpha * t))
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT


    1. 本当は $x(t)$ は整数なんですがここではいったん実数だと思っておきます.
      この辺の処理は物理や化学でもよく出てきます.
      実はこの近似というか合理化もきちんとやっておかないといけません.
      ここではとてもやりきれませんけど.
      本当は時間についてもいろいろ議論があります. 
    2. もっと広く力学系と呼ばれる分野にまで広がります.
      力学系まで行くととんでもなく難しくなるのでとてもここで紹介はできません.
      でも幾何学とも関係してきたり整数論とも関係してきたり,
      射程範囲は広いです. 
  • 厳密解との比較: 放射性物質の崩壊/中高数学駆け込み寺 第 2 回

    こちらに PDF が置いてあります.
    サイトでは見づらい方はこちらをご覧ください.

    近似とは?

    今回アニメーションはありませんが引き続きシミュレーションの話をします.
    今回のテーマは近似についてです.

    シミュレーションというのはいいし,
    近似というのもいいけど,
    実際どのくらいよく近似できてるの?

    こういう問題を考えてみましょう.
    ちゃんと使えばすごく精度はいいですよ,
    というのが今回の主張です.

    • 微分方程式というのはいいけど近似なんて雑なことしたくない!
    • 近似なんてやって意味あるの? どうせ大雑把で使いものにならないんじゃないの?
    • 百歩譲って微分方程式は役に立つとしてもその近似計算は本当に役に立つの?

    こんな疑問に答えるのが今回のテーマです.

    放射性物質の崩壊に関する微分方程式

    で, ちょっと不吉な例ではあります.
    しかしというか残念ながらというか,
    役に立つ例になってしまったので放射性物質の崩壊に関する微分方程式を考えてみましょう.
    いちおう微分方程式から書いておきますね.
    \begin{align}\label{math_refuge00002}
    \frac{dx}{dt}
    =
    – c x.
    \end{align}
    気分だけ説明しておくと放射性物質が単位時間あたりに崩壊するスピードはそのときの放射性物質の量に比例する,
    というのを式で書いた形です.
    $c$ には物理的に大事な意味がありますが式や計算を軽くするためにここでは $c = 1$ にして話を進めます.
    そして前回のようにこの微分方程式を近似した式を書いてみます.
    \begin{align}\label{math_refuge00003}
    \frac{x_{n+1} – x_{n}}{h}
    =
    – x_{n}.
    \end{align}
    これは時間間隔 $h$ が十分小さければ時間 $h$ だけ離れた放射性物質の量 $x_{n+1}$ と $x_{n}$ の関係が上の式で書けることを意味しています.

    近似計算はどのくらい正確なの?

    物理的な話はこのくらいにして,
    この式にしたがって近似計算 (シミュレーション) したとき,
    本当の答えと近似した解がどのくらい近いのかを考えてみます.
    何でこれにしたかというともとの方程式の答えが厳密に書けるからです.
    あなたが指数関数, 特に自然対数の底 $e$ をご存知なら,
    式 \eqref{math_refuge00002} の答えは $x(t) = C e^{- t}$ と書けます.
    指数関数の微分をご存知ならこれを \eqref{math_refuge00002} に代入してみて本当に方程式の答えになっていることを確認してみてください.

    この答えと式 \eqref{math_refuge00003} にしたがって計算した答えがどのくらい一致するか,
    コンピュータに計算させてみましょう.
    結果はこの記事の後半,
    プログラムパートに載せてあります.

    近似なんていい加減ことしたくない!

    あなたは「近似なんて雑なことしていいの?」なんて思っているかもしれません.
    でも記事の後半にある「比較のために重ねる」節のグラフを見るとわかるように,
    厳密な様子ともよく合っているのでちゃんと使えば問題ありません.

    「ちゃんと使う」の「ちゃんと」が難しいと言われればそれはもちろんそうなんですが,
    それは本当に難しい話です.

    せっかくなので分割をどんどん大きくしていって近似の精度が上がっていく様子もグラフにしておきました.
    プログラムが書けるとこういうのもささっとやれます.
    遊びやすくなるのは本当にありがたいですね.
    中高生の頃にこれを知っていたら数学や物理だけじゃなくてプログラミングにもはまっていたと思います.

    やってみるとわかるんですが,
    書いたプログラムで意図通りの結果が出るとけっこう嬉しいです.
    パラメータをいじったり方程式を変えたり,
    ちょろっといじって遊べる要素も増えて遊べる範囲が広がるのもいいですね.

    Python プログラミングに関する資料,
    数値計算・シミュレーションに関する資料は GitHub に上げてありますし,
    今後も講座の進展に合わせて少しずつ増やしていくのでぜひあなたも遊んで遊んでみてください.
    このサイトにもいろいろな情報がありますよ.

    数学的なポイントまとめ

    あと今回の中高数学上のポイントを復習しておきます.
    指数関数とか自然対数の底とか,
    はたまたその微分が出てきました.
    放射性物質の崩壊をちゃんと調べるにはこういう数学が必要なんです.

    次回は経済や生物ネタです!

    次回はまた別の微分方程式を紹介します.
    経済学生物学で出てくる微分方程式ですよ.
    これが終わったらもう少し数学的にちゃんとした話をしましょう.
    まずは微分方程式の射程距離の長さを知ってほしいからです.
    もっとちゃんと数学の説明してほしいというあなた,
    もうちょっと我慢してください.

    今回もアンケートがあるのでぜひ回答をお願いします.

    ではまた次回をお楽しみに!

    プログラミングパート

    放射性物質の崩壊

    一番単純でしかも実際に使われる微分方程式として,
    まずは 1 階の線型常微分方程式を考えよう.

    ちょっと不吉な例であるが放射性物質の崩壊の方程式を紹介する.
    導出をしたければちゃんと物理を勉強してもらう必要がある.
    ここでは物理は省略して数学に集中する.

    \begin{align}
    \frac{dx}{dt}
    =
    – c u.
    \end{align}
    厳密解は $x(t) = C_0 e^{-ct}$ だ.
    初期値を設定すれば $C_0$ はそこから決まる.

    微分を単純に離散化すると次のようになる.

    \begin{align}
    \frac{x_{n+1} – x_{n}}{\Delta t}
    =
    -c x_{n}.
    \end{align}

    時間っぽく見えるよう $\Delta t$ と書いてみた.
    タイプが面倒なので $h$ と書くこともある.
    整理すると次の通り.

    \begin{align}
    x_{n+1}
    =
    x_{n} – c (\Delta t) x_{n}.
    \end{align}

    これに沿って計算したのがいわゆるオイラー法.
    コードに落としていこう.

    厳密解のグラフ

    まずは厳密解 $x(t) = e^{-t}$ をグラフに描こう.
    初期値は 1, 上の定数も $c = 1$ としている.

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    
    
    c = 1
    init = 1
    nt = 100
    
    
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(- c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['exact'])
    
    plt.show()
    

    RESULT

    近次解

    今度は近次解を計算してグラフにしてみよう.

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    
    def radioactive_euler(nt, init = 10):
        dt = 2 / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        x[0] = init
    
        for i in range(1, nt):
            x[i] = x[i-1] - c * dt * x[i-1]
    
        # ベクトル計算で書き直したい
    
        return x
    
    
    c = 1
    init = 1
    nt = 101
    x_approx = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    plt.legend(['approximation'])
    
    plt.show()
    

    RESULT

    比較のために重ねる

    これだけではよくわからないので重ねて描いてみる.

    
    
    x_approx = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(- c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT

    ほぼピッタリ重なって見える.
    シミュレーションの精度はけっこう良さそうだ.

    もちろん細かいことを言えばいろいろあるけれども,
    それはもっと進んだお話なので一旦不問にしておこう.

    細かく刻んだ方がいい近似になるはずだ

    上のグラフは区間を 100 分割して計算した.
    nt = 101 と言うところで分割を指定している.
    101 と言うのは分点数の指定なので実際の分割としては 100 等分なのだ.

    ここで刻みが荒いといい近似にならないことを大雑把に調べておこう.

    2 等分

    
    
    nt = 3
    x_approx = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    nt = 101
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(- c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT

    4 等分

    何となくそれっぽい形にはなっている.
    もちろん精度は大したことない.

    
    
    nt = 5
    x_approx = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    nt = 101
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(- c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULTS

    10 等分

    あまり適当なことを言うのも良くないが, けっこう滑らかに見える.
    精度はまだまだ.

    
    
    nt = 11
    x_approx = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    nt = 101
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(- c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT

    50 等分

    まだ差が目で見える.

    
    
    nt = 51
    x_approx = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    nt = 101
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(- c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT

    1000 等分

    100 は最初にやってあるから一気に大きくしてみた.

    
    
    nt = 1001
    x_approx = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x_approx)
    
    
    nt = 101
    t = np.linspace(0, 2, nt)
    x_exact = init * np.exp(- c * t)
    plt.plot(t, x_exact)
    
    
    plt.legend(['approximation', 'exact'])
    
    plt.show()
    

    RESULT

    参考

    コンピュータで計算するとき,
    いろいろな都合があって分割数を多くすれば単純に精度が上がっていくと言うわけではない.
    桁落ちや丸め誤差のような問題がある.
    それはそれで別途検討する必要があるのだ.

  • はじめに: 全体の大枠を掴もう/中高数学駆け込み寺 第一回

    こちらに PDF が置いてあります.
    サイトでは見づらい方は PDF をご覧ください.

    講座の目的: 微分方程式のシミュレーション

    今回はイントロとして説明することがいくつかあります.
    私が作っている他の講座と同じように,
    この中高数学駆け込み寺では細かいことを詳しくやっていくよりも,
    数学を勉強するモチベーションになるように数学の大きな姿をお見せしていきます.
    で, 実際に何をするのかというと,
    微分方程式のシミュレーションを通じて中高数学の大枠を掴んでいきます.
    もっと具体的に言うと微分方程式のシミュレーションには中高でやる数学がほとんど全て出てきます.

    • この数学はこんなところでこういう風に使うのか!

    使われているシーンを具体的に見てもらってモチベーションにしてもらおう,
    そんな講座です.
    この講座では文字計算には慣れていることを前提にしています.

    微分方程式って何?

    • 微分方程式はどこでどんな風に使われているの?
    • この講座では具体的にどんなことをするの?

    あなたが中高生なら微分は聞いたこともないかもしれません.
    あなたが中高の数学の復習をしたいと思っている大人なら,
    「自分にはまだ微分なんて早いんじゃ?」と思っているかもしれません.
    メインの微分方程式が何だかわからないんじゃ読み進めるのはきついですよね?
    というわけでまずは微分方程式の説明をします.

    微分方程式は理工系の基礎です.
    いわゆる自然法則は微分方程式で表現されることが多いのです.
    これがどこで使われるかというと,
    例えば洪水が起きたときの被害予測に使われます.
    どんな規模で起きた洪水が市街地のどこまで進入してくるか?
    こういうシミュレーションをテレビで見たことがないでしょうか.
    このシミュレーションに使われているのが微分方程式です.

    微分方程式は何に使うの?

    他にも微分方程式はありとあらゆるところで使われています.
    ゲームの CG を自然に見せるためにはまさに自然法則に則って風や水の動きを表現しないといけません.
    だから微分方程式がその背後に隠れています.
    天気予報はいまの気象条件から未来の様子を予測する必要があります.
    この予測にも微分方程式を使っています.
    機械を動かしていると熱くなることがありますね?
    あまりにも熱くなりすぎると機械が暴走してしまうので熱を効率よく逃がす必要があります.
    そのためには熱の流れを考えてその流体の動きをきちんと制御する必要があります.
    ここでは熱の流れの微分方程式を考えてそれをシミュレーションします.
    車を作るとき空気抵抗を調べるためにも流体力学が必要で,
    シミュレーションも使ってモノづくりにも活かしています.
    他にもいちいち挙げきれないくらい身の回りに微分方程式のシミュレーションを使っているモノがあります.

    微分方程式自体をきちんと調べようと思うと大学の数学が必要です.
    でもシミュレーションを中心に考えれば中高の数学でかなりいろいろなことがわかります.
    わかるだけじゃなくて実際に微分方程式を解いて遊んでみることだってできます.
    この自分でも遊べるところまで持っていけるのが大事じゃないかと思っていて.
    それが微分方程式のシミュレーションを選んだ理由の 1 つです.

    さて, ここまでで微分方程式が何で大事かはわかってもらえたでしょうか?
    これをやれば数学をいろいろ遊び倒せそうと思ってもらえたら嬉しいです.
    ちゃんとがんばれば自分でゲームを作ったりもできますしね.
    次はもう少し数学的な内容に踏み込みましょう.

    何はともあれ微分に関する話をやるわけです.
    そして微分は高校でやる内容です.
    あなたは微分に対して嫌な思い出を持っているかもしれません.
    あなたは「中学レベルからやり直したいのにそんなの無理だ!」と思っているかもしれません.
    もっと簡単なところからやってほしいと思っているかもしれません.

    でもこの講座では中高数学の大枠を掴んでもらう講座です.
    そして今回はこの講座の大枠を掴んでもらう講座です.
    どうかもうしばらく辛抱して読み進めてください.

    シミュレーションの実際

    まず何をするのかを見てもらいたいので,
    次のページを開いて 1 番下までスクロールしてください.
    動画の再生ボタンがありますから動画を再生してみましょう.

    これは 1 次元移流方程式のシミュレーション結果です.
    プログラムを書いてコンピュータに計算させ,
    その結果をアニメーションさせています.
    移流方程式は微分方程式なのでそこで微分を使っています.

    コンピュータは数学的に厳密に微分を計算できません.
    極限を取れないからです.
    そうかといって全く何もできないわけじゃなくて.
    実は微分の定義にしたがって微分を近似して計算しています.
    微分を近似すると実はただの引き算 (と割り算) になります.

    さっきのページの上の方,
    式 (2)-(4) を見てください.
    こっちにも同じ式を書いておきますね.

    \begin{align}
    \frac{\partial u}{\partial x}
    \approx
    \frac{u(x+\Delta x)-u(x)}{\Delta x}.
    \end{align}

    \begin{align}
    \frac{u_i^{n+1}-u_i^n}{\Delta t} + c \frac{u_i^n – u_{i-1}^n}{\Delta x}
    =
    0.
    \end{align}

    \begin{align}
    u_i^{n+1}
    =
    u_i^n – c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n). \label{math_refuge00001}
    \end{align}

    見ての通り加減乗除の四則演算しかしていません.
    特に上のページの (4) にあたる式 \eqref{math_refuge00001} に注目してください.
    実はこの式にしたがって計算した結果をアニメーションでお見せしたのがさっきの動画です.

    細かいことは追々やっていくとして何をしているか大雑把にいうと,
    ベクトルなり数列なりの計算です.
    また式 \eqref{math_refuge00001} を見てください.
    右上にある添字 $n$ に注目すれば左辺は $n+1$ で右辺は $n$ です.
    これは数列の漸化式とも思えますね?

    数列も高校でやる内容です.
    あなたが中学生ならちんぷんかんぷんでしょう.
    でも実際に高校でやることがモノづくりなり何なりで
    役に立っていることが感じてもらえればとりあえず OK です.
    ちなみに移流方程式は上でもちょっと出てきた流体力学で出てくる微分方程式です.

    もっと言うなら文字式の計算がゴリゴリ出てきていることも注意した方がいいでしょうか?
    文字式できないとこの計算全く追えないですからね.
    細かいことを言い出すときりがないので,
    いったん今回はこのくらいにしておきましょう.
    まとめると今回は次のポイントを掴んでもらえば十分です.

    今回のポイント

    • 微分方程式のシミュレーションをやっていく.
    • シミュレーションは実社会でいろいろな応用がある
    • 微分といっても近似を使うから結局は加減乗除しか使わない.
    • シミュレーションでは文字式の計算, ベクトル, 数列など中高の数学をバリバリ使っている.

    今回くらいのレベルでもベクトルで言うと
    100 次元ベクトルとかそういうレベルの計算が必要です.
    手計算ではやってられないのでプログラムを書いて計算して,
    さらにプログラムを書いて図にしたりアニメーションにしたりしています.
    このプログラミングについてはどこまで深掘りするかは未定です.
    要望が多いなら別にきちんとやろうかとも思っています.

    今回お見せしたのは変数が 2 つある偏微分方程式でちょっと難しいです.
    変数が 2 つあるとややこしいので次回以降,
    しばらく変数が 1 つしかない常微分方程式というのを見ていく予定です.

    最後にアンケートをお願いしています.
    改善に繋げていきたいのでぜひコメントをお願いします.

    では次回をお楽しみに!

  • 数列とは何か?プログラム解説 中高数学駆け込み寺

    1 階の線型常微分方程式

    復習でもう 1 回.

    一番単純でしかも実際に使われる微分方程式としてまずは 1 階の線型常微分方程式を考えよう.
    ちょっと不吉な例であるが放射性物質の崩壊の方程式を紹介する.
    導出をしたければちゃんと物理を勉強してもらう必要がある.
    ここでは物理は省略して数学に集中する.

    \begin{align}
    \frac{dx}{dt} = – c x.
    \end{align}
    厳密解は $x = C_0 e^{-ct}$ だ.
    初期値を設定すれば $C_0$ はそこから決まる.

    微分を単純に離散化すると次のようになる.

    \begin{align}
    \frac{x_{n+1} – x_{n}}{\Delta t}
    =
    -c x_{n}.
    \end{align}

    $\Delta t$ は $h$ と書くこともある.
    整理すると次の通り.

    \begin{align}
    x_{n+1}
    =
    x_{n} – c (\Delta t) x_{n}.
    \end{align}

    これに沿って計算したのがいわゆるオイラー法.
    次のセルではこれをコードに落としている.

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    
    import base64
    img_file = "tmp/image.tmp.png"
    
    def radioactive_euler(nt, init = 10):
        dt = 2 / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        x[0] = init
    
        for i in range(1, nt):
            x[i] = x[i-1] - c * dt * x[i-1]
    
        # ベクトル計算で書き直したい
    
        return x
    
     # 近似解
    c = 1
    nt = 101
    init = 5
    x1 = radioactive_euler(nt, init)
    plt.plot(np.linspace(0, 2, nt), x1)
    
     # 厳密解
    t = np.linspace(0, 2, nt)
    x2 = init * np.exp(- c * t)
    plt.plot(t, x2)
    
     # 凡例
    plt.legend(['approximation', 'exact'])
    plt.show()
     # 描画
    #plt.savefig(img_file)
    #b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
    #img_str = "" % (b64)
    #print(img_str)
    

    RESULT

    2 階の常微分方程式

    こちらが単振動.

    次に 2 階の常微分方程式を紹介しよう.
    高校の物理で出てくるばねの振動(単振動)がまさにこの例だ.
    項を増やすと減衰振動になったり、外力をつけたりといろいろなケースがある.
    まずは一番単純な式を考えよう.

    \begin{align}
    \frac{d^2 x}{dt^2}
    =
    – \omega^2 x.
    \end{align}

    さっきのオイラー法なりルンゲ-クッタ法なりは 1 階の方程式に対する計算法なので直接は使えない.
    この場合は中間処理として $v = dx/dt$ を置いて計算すればいい.
    これは単なる数値計算の便法ではない.
    速度の意味もあるから, という表面的な理由ではなくもっと深く解析力学の文脈で物理としても大事な視点だ.
    もっといえばシンプレクティック計算法などもっといい計算法にも発展する.

    とりあえずオイラー法で計算したい.
    まずは微分方程式自体を書き直す.

    \begin{align}
    \frac{dx}{dt}
    =
    v, \quad
    \frac{dv}{dt}
    =
    – \omega^2 x.
    \end{align}

    これをオイラー法で近似しよう.

    \begin{align}
    x_{n+1}
    =
    x_{n} + h v_{n}, \quad
    v_{n+1}
    =
    v_{n} – h \omega^{2} x_{n}.
    \end{align}

    オイラー法をコードに落とす.

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib import animation
    
    def harmonic_euler(nt, init = (5, 0)):
        dt = t_range / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        v = np.zeros(nt)
        x[0] = init[0]
        v[0] = init[1]
    
        for i in range(1, nt):
            x[i] = x[i-1] + dt * v[i-1]
            v[i] = v[i-1] - dt * (omega ** 2) * x[i-1]
    
        # ベクトル計算で書き直したい
    
        return (x, v)
    
    
    omega = 2 * np.pi
    nt = 101
    t_range = 2
    init = (5, 0)
    harm = harmonic_euler(nt, init)
    t = np.linspace(0, 2, nt)
    
     # 厳密解
    x_exact = init[0] * np.cos(- omega * t)
    v_exact = - omega * init[0] * np.sin(omega * t)
    
     # グラフ描画
    plt.subplot(3, 1, 1)
    plt.title('x-t graph')
    plt.plot(np.linspace(0, 2, nt), harm[0])
    plt.plot(t, x_exact)
    plt.legend(['x approximation', 'x exact'])
    
    plt.subplot(3, 1, 2)
    plt.title('v-t graph')
    plt.plot(np.linspace(0, 2, nt), harm[1])
    plt.plot(t, v_exact)
    plt.legend(['v approximation', 'v exact'])
    
    plt.subplot(3, 1, 3)
    plt.title('phase space')
    plt.plot(harm[0], harm[1])
    
     # 描画
    plt.tight_layout()
    plt.show()
    #plt.savefig(img_file)
    #b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
    #img_str = "" % (b64)
    #print(img_str)
    

    RESULT

    オイラー法では近似の精度が悪い

    見ての通り時間が進むごとに誤差が大きくなる.
    nt を大きくすると少しはましになる.
    実際に上のコードで nt を大きくして再計算してみてほしい.

    あとまずいのは phase space の図だ.
    この系はエネルギーが保存する系だから相空間内の軌道が閉じてほしいのにそうなっていない.
    シンプレクティックにやれば解消できるようだが, とにかくここではよろしくない.

    この方程式でオイラー法はよろしくないことがわかった.
    とりあえずルンゲ-クッタでやってみよう.

    ルンゲ-クッタ法

    とりあえず近似式を書く.

    \begin{align}
    x_{n+1}
    &=
    x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), \
    t_{n+1}
    &=
    t_{n} + h, \
    k_{1}
    &=
    v_{n}, \
    k_{2}
    &=
    v_{n} + \frac{h}{2} k_{1}, \
    k_{3}
    &=
    v_{n} + \frac{h}{2} k_{2}, \
    k_{4}
    &=
    v_{n} + h k_{3}.
    \end{align}

    次が $v$ の式.

    \begin{align}
    v_{n+1}
    &=
    v_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), \
    t_{n+1}
    &=
    t_{n} + h, \
    k_{1}
    &=
    – \omega^2 x_{n}, \
    k_{2}
    &=
    x_{n} – \frac{h}{2} \omega^2 k_{1}, \
    k_{3}
    &=
    x_{n} – \frac{h}{2} \omega^2 k_{2}, \
    k_{4}
    &=
    x_{n} – h \omega^2 k_{3}.
    \end{align}

    
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib import animation
    
    def harmonic_rk(nt, init = 10):
        dt = t_range / (nt - 1)
        # 初期条件設定
        x = np.zeros(nt)
        v = np.zeros(nt)
        x[0] = init[0]
        v[0] = init[1]
    
        def fx(t, x, v):
            return v
    
        def fv(t, x, v):
            return - (omega ** 2) * x
    
        # ベクトル計算で書き直したい
        for i in range(1, nt):
            xk1 = fx(dt * (i - 1), x[i-1], v[i-1])
            vk1 = fv(dt * (i - 1), x[i-1], v[i-1])
    
            xk2 = fx(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)
            vk2 = fv(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)
    
            xk3 = fx(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)
            vk3 = fv(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)
    
            xk4 = fx(dt * (i - 1),   x[i-1] + xk3 * dt,     v[i-1] + vk3 * dt)
            vk4 = fv(dt * (i - 1),   x[i-1] + xk3 * dt,     v[i-1] + vk3 * dt)
    
            x[i] = x[i-1] + dt / 6 * (xk1 + 2 * xk2 + 2 * xk3 + xk4)
            v[i] = v[i-1] + dt / 6 * (vk1 + 2 * vk2 + 2 * vk3 + vk4)
    
        return (x, v)
    
    
    omega = 2 * np.pi
    nt = 101
    t_range = 2
    init = (5, 0)
    harm = harmonic_rk(nt, init)
    t = np.linspace(0, 2, nt)
    
    
    x_exact = init[0] * np.cos(- omega * t)
    v_exact = - omega * init[0] * np.sin(omega * t)
    
    
    plt.subplot(3, 1, 1)
    plt.title('x-t graph')
    plt.plot(np.linspace(0, 2, nt), harm[0])
    plt.plot(t, x_exact)
    plt.legend(['x approximation', 'x exact'])
    
    plt.subplot(3, 1, 2)
    plt.title('v-t graph')
    plt.plot(np.linspace(0, 2, nt), harm[1])
    plt.plot(t, v_exact)
    plt.legend(['v approximation', 'v exact'])
    
    plt.subplot(3, 1, 3)
    plt.title('phase space')
    plt.plot(harm[0], harm[1])
    
     # 描画
    plt.tight_layout()
    plt.show()
    #plt.savefig(img_file)
    #b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
    #img_str = "" % (b64)
    #print(img_str)
    

    RESULT

    ルンゲ-クッタ短評

    今度の一致具合はなかなかよさそう.
    あくまで見た目の感じではあるけれども.

    良くなったのは一番下の図, 相空間軌道だ: ちゃんと閉じてくれた.