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

!pip install --upgrade sympy

import sympy as sp
from sympy.plotting import plot
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 sp.latex(expr, **options)

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

x = sp.Symbol('x')
expr = x**2-12*x+8
display(expr)

Collecting sympy
[?25l  Downloading https://files.pythonhosted.org/packages/ce/5b/acc12e3c0d0be685601fc2b2d20ed18dc0bf461380e763afc9d0a548deb0/sympy-1.5.1-py2.py3-none-any.whl (5.6MB)
     |████████████████████████████████| 5.6MB 3.3MB/s
[?25hRequirement already satisfied, skipping upgrade: mpmath>=0.19 in /usr/local/lib/python3.6/dist-packages (from sympy) (1.1.0)
Installing collected packages: sympy
  Found existing installation: sympy 1.1.1
    Uninstalling sympy-1.1.1:
      Successfully uninstalled sympy-1.1.1
Successfully installed sympy-1.5.1



<IPython.core.display.Javascript object>

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")

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>

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>

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'))

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'))

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'))