2020-04-19_introduction オンラインプログラミング勉強会

この記事は14分で読めます

このサイトは学部では早稲田で物理を, 修士では東大で数学を専攻し, 今も非アカデミックの立場で数学や物理と向き合っている一市民の奮闘の記録です. 運営者情報および運営理念についてはこちらをご覧ください.

中高の数学の復習から専門的な数学・物理までいろいろな情報を発信しています.
中高数学に関しては自然を再現しよう役に立つ中高数学 中高数学お散歩コース
大学数学に関しては現代数学観光ツアーなどの無料の通信講座があります.
その他にも無料の通信講座はこちらのページにまとまっています.
ご興味のある方はぜひお気軽にご登録ください!


オンライン数学勉強会用イントロ

事前準備:sympy インストールとチェック

from sympy import *
from IPython.display import display

def custom_latex_printer(expr, **options):
    from IPython.display import Math, HTML
    from google.colab.output._publish import javascript
    url = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/MathJax.js?config=TeX-AMS_CHTML"
    javascript(content="""window.MathJax = {
        tex2jax: {
            inlineMath: [ ['$','$'] ],
            processEscapes: true
        }
        };""")
    javascript(url=url)
    return latex(expr, **options)

init_printing(use_latex="mathjax", latex_printer=custom_latex_printer)

x = Symbol('x')
expr = x**2-12*x+8
display(expr)
<IPython.core.display.Javascript object>

$$x^{2} – 12 x + 8$$

はじめの注意

  • 吃音があるので聞きにくいと思います。
  • 言いにくい場合は文字で書きます。
  • いろいろ試験的な取り組みです。
  • この資料はあとで配ります。

大方針

  • Python でプログラミング
  • 将来に向けたいろいろなテストを兼ねる
  • 中高生または中高数学復習用に作ったコンテンツのテスト・流用
  • 試験的に Google Colabolratory を利用
  • いま使っている「これ」
  • ローカルの環境構築がいらない
  • TeX の形式を使うと式もきれいに書ける:$$\int_\Omega f(x) \, dx$$
  • 何回くらい・どのくらいの期間になるかは不明
  • なるべくゆっくり・ゆるく進める
  • 日々いろいろなタスクがある大人向けの勉強会なのでゴリゴリにやると疲れ果てて続かない
  • その代わり長期戦でのんびりやる
  • 予習はやりたければやる
  • 無理がない範囲で、長続きさせることを第一に
  • 復習は勉強会の中でやる
    • その分ペースも遅くなる
    • その遅さが気になるなら予習をしよう
  • 復習を自力でやれるくらいならはじめから独学できていると思う

勉強のコツ

  • 「すぐに理解する・できる」という幻想を捨てる
  • 同じことを何度も繰り返して、少しずつ慣れる
  • 究極的には「独学」が必要
  • そのサポートをするのがこの勉強会の目的
  • 細部が詳しい本はたくさんあり、いろいろなレベルの本がたくさんある
  • この勉強会の当面の主目的は概要を掴むこと
  • 細かいことが気になったら都度調べるなり質問するなりしてもらう
  • そしてすぐにわからなくても気にしない

進め方の方針

  • 細かいことはさておき, プログラミングでできることをゴリゴリ紹介する
  • まずはお絵描き中心:微分方程式を解く
  • 後の内容も「中高数学のネタを視覚化する」方向で進める。
  • 微分方程式に関して
    • 高校数学の最終目的地の微分・積分に直結する
    • 微分・積分は統計学をはじめとしたその他いろいろな応用でも主力
    • 離散化すれば四則演算でしかないので微分・積分をダイレクトに数値計算する形で進める
    • 今回、イメージづくりで波動方程式・拡散方程式を見せる
    • 波と拡散というイメージと視覚がマッチした対象だから
    • まずは概要を掴むのが目的なので本質的な数学部分が難しくなるのは気にしない
    • 何をするにも慣れが必要なので, とにかく浴びて慣れてもらうのが目的
  • お絵描き+数値計算で Python プログラミングの気分を掴む
  • 少しずつプログラムの詳細を見ていく
    • 暴力的な量の四則演算を進めるだけなので「数学」の話は最低限に留めながら進める
  • ある程度慣れてきた段階でプログラムと並行して数学の話をする

基礎コンテンツの紹介

  • 中高数学をプログラミングを軸にまとめた。
  • この講座ではプログラミング、特に Python については深くは解説しない
  • 必要ならその時々の新しい本で勉強すること。
  • 古い本を選んでしまうと、最新の環境でその本に書かれたプログラムが動かない可能性がある。

講座の構成

  1. 数学に関わる Python の基礎
    • 重要なライブラリ numpy のまとめ
    • 重要なライブラリ matplotlib のまとめ
    • 重要なライブラリ sympy のまとめ
    • 上記ライブラリを使った数学プログラミングのまとめ
  2. 線型代数:ベクトルと行列
  3. 微分積分
  4. 確率
  5. 統計
  6. 常微分方程式
  7. 偏微分方程式

大まかな説明

  • 最初の第 1 章:Python を使ってどんなことができるかを説明
    • 特に Python 自体の基本的な機能とライブラリの使い方を説明
    • 必要なコードは各箇所に書かれているので、ここは飛ばしてすぐに本編に進んでも問題ない
  • 線型代数以降の各章から、数学・プログラミングが混然一体になった本編がはじまる。
  • 常微分方程式の章を先に読んでみるのもおすすめします。
  • 「大事なことは何度でも」の精神

注意・当面の進めるイメージ

  • 下の 2 つを先にざっと眺める
  • 既に無料講座として公開している常微分方程式の章を少しずつ詳しく進める
  • 微分積分をやる

環境構築

Google Colab

  • いわゆる gmail のアカウントがあれば使える。
  • 「google colaboratory」でググってトップに出てくるページを開けば開くはず
  • 公式のチュートリアルもあるので、それを見よう

Python のインストール

  • インストールのことを「環境構築」とも呼ぶことにする
  • 将来の対応を考えて、勉強会では Google Colaboratory で進めてみる
  • ローカル(自分の PC)にインストールしたいなら、適当にやる
  • 経験上、「素人」の状態で環境構築は本当に大変
  • プログラミング・環境構築に慣れていない状態で数学と絡めた環境構築がしたいなら、Anaconda でインストールするのが楽
  • Google Colab 前提なのでこれ以上詳しくは触れない
  • 本格的にプログラミングするならローカルに環境を作るのは必須

本題:数学とアニメーション

定積分のアニメーション

  • 定積分はグラフが囲む領域の面積
  • いわゆる区分求積法:次のように定義する

\begin{align}
\int_0^1 f(x) dx
&=
\lim_{n \to \infty} \frac{1}{n} \sum_{k=0}^{n-1} f \left( \frac{k}{n} \right), \
n = 2
&\Rightarrow
\frac{1}{2} \left(f\left(\frac02\right) + f\left(\frac12\right)\right) \
n = 3
&\Rightarrow
\frac{1}{3} \left(f\left(\frac03\right) + f\left(\frac13\right) + f\left(\frac23\right)\right)
\end{align}

  • 右辺に注目する
    • $\frac{1}{n}$ は単なる割り算
    • $\sum_{k=0}^{n-1}$ は $n$ 項足しているだけ
    • $\lim$ は $n$ をどんどん大きくしている
  • 数値計算でやること
    • 十分大きな $n$ で右辺を計算する
    • だいたい面積が近似できることを確認する

アニメーション用の基本関数

  • 次の関数でいろいろやっている:今回詳しくは触れない
import numpy as np
import matplotlib.pyplot as plt

def scatter_integral_body(num_data, is_left, f, param):
    # 描画領域の指定
    fig = plt.figure(figsize=(param["figsize_x"], param["figsize_y"]))
    subplot = fig.add_subplot(1, 1, 1)
    subplot.set_xlim(param["x_min"], param["x_max"])
    subplot.set_ylim(param["y_min"], param["y_max"])

    # 参考: グラフを折れ線で描く. 折れ線近似の様子を見たい場合に使う.
    #linex = np.linspace(param["x_min"], param["x_max"], num_data) 
    #subplot.plot(linex, f(linex), color='blue')

    # 長方形の描画と面積の近似値計算
    area = 0
    step = (param["x_max"] - param["x_min"]) / num_data
    for x0 in np.arange(param["x_min"], param["x_max"], step):
        x = x0 if is_left else x0 + step
        rect = plt.Rectangle((x0, 0), step, f(x), alpha=param["alpha"])
        subplot.add_patch(rect)
        area += step * f(x)

    subplot.text(param["text_x"], param["text_y"], ('area = %f' % area))

    # 関数と三角関数の描画:上の supplot の処理の後に置かないと上の長方形が描かれない
    vf = np.vectorize(f)

    x_for_square = np.linspace(param["x_min"], param["x_max"], param["num_max"])
    y0 = vf(x_for_square)
    plt.plot(x_for_square, y0, color="red")

底辺と高さ 1 の直角三角形の面積

  • 底辺×高さ / 2 = 1/2
from ipywidgets import interact

def scatter(num_data, is_left):
    def f(x):
        return x

    param = {
        "x_min": 0.0,
        "x_max": 1.0,
        "y_min": 0.0,
        "y_max": 1.0,
        "figsize_x": 4,
        "figsize_y": 4,
        "alpha": 0.5,
        "text_x": 0.2,
        "text_y": 0.7,
        "num_max": 50
    }
    scatter_integral_body(num_data, is_left, f, param)

num_max = 50
interact(scatter, num_data=(2, num_max, 1), is_left=True)
interactive(children=(IntSlider(value=26, description='num_data', max=50, min=2), Checkbox(value=True, descrip…





<function __main__.scatter>

参照用

\begin{align}
\int_0^1 f(x) dx
&=
\lim_{n \to \infty} \frac{1}{n} \sum_{k=0}^{n-1} f \left( \frac{k}{n} \right), \
n = 2
&\Rightarrow
\frac{1}{2} \left(f\left(\frac02\right) + f\left(\frac12\right)\right) \
n = 3
&\Rightarrow
\frac{1}{3} \left(f\left(\frac03\right) + f\left(\frac13\right) + f\left(\frac23\right)\right)
\end{align}

二次関数

  • $x^2$ を $[0,2]$ で積分
  • 面積は $\int_0^2 x^2 dx = 8/3 = 2.6666…$
from ipywidgets import interact

def scatter(num_data, is_left):
    def f(x):
        return x**2

    param = {
        "x_min": 0.0,
        "x_max": 2.0,
        "y_min": 0.0,
        "y_max": 4.0,
        "figsize_x": 4,
        "figsize_y": 4,
        "alpha": 0.5,
        "text_x": 0.2,
        "text_y": 3.0,
        "num_max": 50
    }
    scatter_integral_body(num_data, is_left, f, param)

num_max = 50
interact(scatter, num_data=(2, num_max, 1), is_left=True)
interactive(children=(IntSlider(value=26, description='num_data', max=50, min=2), Checkbox(value=True, descrip…





<function __main__.scatter>

定積分の大雑把なまとめ

  • どんどん細かくしていくと、細かさに応じて近似がよくなる。
  • 近似を上げていった最果てで厳密な値が得られるとみなす。
  • 視覚的に確認すると定義の意図が見えやすくなる
  • 定義の意図がわかっていないとプログラムを書くのも難しい

1 次元の波動方程式

  • ひもの左端を揺らしたときのひもの振動(波)の様子を見る
  • 波が右端まで行くと反射するところも見られる

偏微分方程式としての波動方程式

\begin{align}
u_{tt} = c^2 u_{xx}, \quad
\frac{\partial^2 u}{\partial t^2}
=
c^2 \frac{\partial^2 u}{\partial x^2}.
\end{align}

  • ここで $c$ は定数で物理的には波の速度。
  • 物理だとこれをどうやって導出するかも問われる。
  • ここではそうした議論はせず、とにかく数値的に解いてみて波を表している様子を眺める

\begin{align}
\frac{\partial u}{\partial t}
&=
\lim_{\Delta t \to 0}\frac{u(t + \Delta t, x) – u(t,x)}{\Delta t} \
&\simeq
\frac{u(t + \Delta t, x) – u(t,x)}{\Delta t}
\end{align}

\begin{align}

\frac{\partial^2 u}{\partial t^2}

\frac{\partial}{\partial t} \frac{\partial u}{\partial t}=
\frac{u_t(t+\Delta t,x) – u_t(t,x)}{\Delta t}
&\simeq
\frac{\frac{u(t + \Delta t, x) – u(t,x)}{\Delta t} – \frac{u(t , x) – u(t- \Delta t,x)}{\Delta t}}{\Delta t}
\end{align}

最終的に計算する式

\begin{align}
u(t + \Delta t, x ) =
2 u(t,x) – u(t – \Delta t, x) +
\left(\frac{c \Delta t}{\Delta x}\right)^2
\left(u(t, x + \Delta x) – 2 u(t,x) + u(t, x – \Delta x) \right).
\end{align}

  • 次の時刻 $t + \Delta t$ の値を計算するのに現在時刻 $t$ とひとつ前の時刻 $t – \Delta t$ の値を使っている。
  • 最終的には四則演算だけ
from IPython.display import HTML
from matplotlib import animation
from tqdm import tqdm
import numpy as np
import matplotlib.pyplot as plt

lx = 1  # 計算領域の長さ
nx = 41 # 領域の分割数
xs = np.linspace(0, lx, nx)

def wave(nx):
    dx = lx / (nx - 1) # 空間方向の刻み
    nt = 500           # 計算回数
    dt = 0.01          # 時間刻み
    c = 1              # 波の速度
    f = 2.0            # 強制振動の周波数

    # 初期条件
    u0 = np.zeros(nx)

    # 結果の配列と nt, nx の次元を持つ配列で初期化
    us = np.zeros((nt, len(u0)))

    # 初期条件を結果の配列の各行にコピー
    us[:,:] = u0.copy()

    alpha = (c * dt / dx)**2

    for i in tqdm(range(1, nt-2)):
        t = i * dt
        # 左端を強制振動
        us[i, 0] = np.sin(2 * np.pi * f * t)
        us[i+1, 0] = np.sin(2 * np.pi * f * t)

        # 時間発展
        us[i+1, 1:-1] = 2 * us[i, 1:-1] - us[i-1, 1:-1] \
                        + alpha * (us[i, 2:] - 2 * us[i, 1:-1] + us[i, :-2])
        # 右端の境界値指定
        us[i+1, nx-1] = 0

    return us

# 解
us = wave(nx)

fig = plt.figure();
# グラフの軸の設定
ax = plt.axes(xlim=(0, lx), ylim=(-5, 5));
line, = ax.plot([], [], lw=2);

def animate(u):
    line.set_data(xs, u)
    return line

anim = animation.FuncAnimation(fig, animate, frames=us, interval=50)
#anim.save('07-pde_03_wave_ex01.tmp.mp4', writer="ffmpeg")
plt.close(anim._fig)
HTML(anim.to_jshtml(default_mode='reflect'))
100%|██████████| 497/497 [00:00<00:00, 34649.84it/s]

1 次元の拡散方程式

  • 真ん中に置いておいた物質が周囲に拡散していく様子を見る

偏微分方程式としての拡散方程式

\begin{align}
\frac{\partial u}{\partial t}= \nu \frac{\partial^2 u}{\partial x^2}.
\end{align}

最終的に計算する式

\begin{align}
u(t+ \Delta t, x) =
u(t, x) + \nu \frac{\Delta t}{(\Delta x)^2} (u(t, x + \Delta x) – 2u(t,x) + u(t, x – \Delta x)).
\end{align}

  • これも四則演算だけ。
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import HTML
from matplotlib import animation

lx = 2  # 計算領域の長さ
nx = 41 # 領域の分割数
xs = np.linspace(0, lx, nx)

def diffusion(nx):
    dx = lx / (nx - 1) # 空間方向の刻み
    dt = 0.001         # 時間刻み
    nt = 200           # 計算回数
    nu = 1.0           # 拡散係数

    # 初期条件:中央以外すべて 0
    u0 = np.zeros(nx)
    # 真ん中に物質を置く
    x_center = int((nx - 1) / 2)
    u0[x_center-2: x_center+2] = 1

    # 結果の配列と nt, nx の次元を持つ配列で初期化
    us = np.zeros((nt, len(u0)))

    # 初期条件を結果の配列の各行にコピー
    us[:,:] = u0.copy()

    for i in range(0, nt-1):
        us[i+1, 1:-1] = us[i, 1:-1] + nu * dt / dx**2 *  (us[i, 2:] - 2 * us[i, 1:-1] + us[i, :-2])

    return us

us = diffusion(nx)

fig = plt.figure();
# グラフの軸の設定
ax = plt.axes(xlim=(0, lx), ylim=(0, 1.2));
line, = ax.plot([], [], lw=2);

def animate(u):
    line.set_data(xs, u)
    return line

anim = animation.FuncAnimation(fig, animate, frames=us, interval=50)
#anim.save('07-pde_04_diffusion_ex01.tmp.mp4', writer="ffmpeg")
plt.close(anim._fig)
HTML(anim.to_jshtml(default_mode='reflect'))

2 次元の波動方程式

  • イメージとしては膜の振動
  • 真ん中を振動させてその波の伝播を見る

偏微分方程式としての波動方程式

\begin{align}
u_{tt} = c^2 (u_{xx} + u_{yy}).
\end{align}

  • ここで $c$ は定数で物理的には波の速度
  • 1 次元のときとの違いは $u_{yy}$ の追加
  • 物理的には次元に依存する議論もあってそれほど簡単ではない

最終的に計算する式

\begin{align}
&u(t + \Delta t, x, y) \
&=
2 u(t,x,y) – u(t – \Delta t, x, y) \
\quad &+ \left(\frac{c \Delta t}{\Delta x}\right)^2 \left(u(t, x + \Delta x,y) – 2 u(t,x,y) + u(t, x – \Delta x,y) \right) \
\quad &+ \left(\frac{c \Delta t}{\Delta y}\right)^2 \left(u(t, x,y + \Delta y) – 2 u(t,x,y) + u(t, x,y – \Delta y) \right) \
\end{align}

  • 右辺の時刻は $t$ と $t – \Delta t$、つまり過去の時刻しか出てこない
  • 過去の時間での情報さえわかれば未来の挙動がわかる
from IPython.display import HTML
from matplotlib import animation
from mpl_toolkits.mplot3d import Axes3D
from tqdm import tqdm
import numpy as np
import matplotlib.pyplot as plt

lx = 1  # 計算領域の長さ
ly = 1  # 計算領域の長さ
nx = 11 # 領域の分割数
ny = 11 # 領域の分割数
nx_center = int((nx - 1) / 2)
ny_center = int((ny - 1) / 2)
xs = np.linspace(0, lx, nx)
ys = np.linspace(0, ly, ny)

x1d = np.linspace(0, lx, nx)
y1d = np.linspace(0, ly, ny)
xs, ys = np.meshgrid(x1d, y1d)

def wave(nx, ny):
    dx = lx / (nx - 1)
    dy = ly / (ny - 1)
    nt = 300           # 計算回数
    dt = 0.5 * dx * dy # 時間刻み
    c = 1              # 波の速度
    f = 2.0            # 強制振動の周波数

    # 初期条件
    u0 = np.zeros((nx, ny))

    # 結果の配列と nt, nx, ny の次元を持つ配列で初期化
    us = np.zeros((nt, nx, ny))

    # 初期条件を結果の配列の各行にコピー
    us[:,:,:] = u0.copy()

    alpha_x = (c * dt / dx)**2
    alpha_y = (c * dt / dy)**2

    for i in tqdm(range(1, nt-2)):
        t = i * dt
        # 中心を強制振動
        us[i, nx_center, ny_center] = np.sin(2 * np.pi * f * t)
        us[i+1, nx_center, ny_center] = np.sin(2 * np.pi * f * t)

        # 時間発展
        us[i+1, 1:-1, 1:-1] = 2 * us[i, 1:-1, 1:-1] - us[i-1, 1:-1, 1:-1] \
                            + alpha_x * (us[i, 2:, 1:-1] - 2 * us[i, 1:-1, 1:-1] + us[i, :-2, 1:-1]) \
                            + alpha_y * (us[i, 1:-1, 2:] - 2 * us[i, 1:-1, 1:-1] + us[i, 1:-1, :-2])
        # 境界値指定
        us[i+1, 0, :] = 0
        us[i+1, nx-1, :] = 0
        us[i+1, :, 0] = 0
        us[i+1, :, ny-1] = 0

    return us

# 解
us = wave(nx, ny)

fig = plt.figure();
ax = fig.add_subplot(111, projection='3d')
ax.set_xlim(0, lx)
ax.set_ylim(0, ly)
ax.set_zlim(-2, 2)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$u$')
surf = ax.plot_surface(xs, ys, us[0], rstride=1, cstride=1, linewidth=0, cmap='jet')

def animate(u):
    return ax.plot_surface(xs, ys, u, rstride=1, cstride=1, linewidth=0, cmap='jet')
    #surf.set_data(u)
    #return surf

anim = animation.FuncAnimation(fig, animate, frames=us, interval=50)
#anim.save('07-pde_06_wave_2dim_ex01.tmp.mp4', writer="ffmpeg")
plt.close(anim._fig)
HTML(anim.to_jshtml(default_mode='loop'))

微分方程式の大雑把なまとめ

  • (いろいろな都合によって)物理では無限に細かいところで方程式を立てている
  • そのままでは計算機に計算させられない
  • 微分係数・導関数を定義によって有限化:差と商にわける。
  • 有限に落とした部分は無限に細かくしていけば元の微分係数・導関数を近似できていると期待する

次回以降

  • 常微分方程式のところをのんびり進める
  • まずは微分方程式で何ができるのかをのんびり
  • 微分方程式で記述できる現象を説明したあと、近似計算の数理を追う
  • プログラムに落とし込む工夫を見る
  • 数値計算プログラムを大雑把に眺める
  • 講座・コンテンツ本体を眺める

質疑応答

  • 今の時点での目標が何かとか
  • こんな感じで進めてほしいとか
  • こんなネタを扱ってほしいとか

要望メモ

  • 小さな課題が欲しい。
  • 成功体験を積もう。
  • 予習の範囲を出す。
  • 式とグラフの対応。いろいろお絵描きしてみる。どうやってプログラムに落とすか確認する。

アンケート

毎回アンケートを取っています.
質問や要望がある場合もこちらにどうぞ.

アンケートは匿名なので気楽にコメントしてください.
直接返事してほしいことがあれば,
メールなど適当な手段で連絡してください.


中高の数学の復習から専門的な数学・物理までいろいろな情報を発信しています.
中高数学に関しては自然を再現しよう役に立つ中高数学 中高数学お散歩コース
大学数学に関しては現代数学観光ツアーなどの無料の通信講座があります.
その他にも無料の通信講座はこちらのページにまとまっています.
ご興味のある方はぜひお気軽にご登録ください!

  • このエントリーをはてなブックマークに追加
  • LINEで送る

関連記事

  • コメント (0)

  • トラックバックは利用できません。

  1. この記事へのコメントはありません。

このサイトについて

数学・物理の情報を中心にアカデミックな話題を発信しています。詳しいプロフィールはこちらから。通信講座を中心に数学や物理を独学しやすい環境づくりを目指して日々活動しています。
  • このエントリーをはてなブックマークに追加
  • LINEで送る

YouTube チャンネル登録

講義など動画を使った形式の方が良いコンテンツは動画にしています。ぜひチャンネル登録を!

メルマガ登録

メルマガ登録ページからご登録ください。 数学・物理の専門的な情報と大学受験向けのメルマガの 2 種類があります。

役に立つ・面白い記事があればクリックを!

記事の編集ページから「おすすめ記事」を複数選択してください。