{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### 2023年度計算機演習A・B\n", "\n", "# 第14回：フラクタル1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**フラクタル（fractal）**とは、「自己相似性」（全体と部分が相似となる性質）を持つような図形のことです。\n", "\n", "今回と次回の授業では、プログラミングによっていくつかのフラクタルの描画を行います。\n", "\n", "※フラクタルは通常、ある特定の操作を無限回繰り返すことで得られますが、プログラミングでそれは実現できないため、以下では有限回繰り返した結果を描画することにします。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. コッホ曲線\n", "\n", "**コッホ曲線（Koch curve）**は特に有名なフラクタルであり、次のようにして得られます。\n", "\n", "
\n", "\n", "- ステップ1：与えられた線分を3等分し、中央の線分を1辺とする正三角形を外側に描いて内側の辺を消す（左上→右上）。\n", "\n", "- ステップ2：得られた4本の線分に対して同様の操作を行う（右上→左下）。\n", "\n", "- ステップ3：得られた16本の線分に対して同様の操作を行う（左下→右下）。\n", "\n", "- …\n", "\n", "以下では、3つの関数を定義することによって、コッホ曲線を描画するプログラムを作成します。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "まず、複数の点からなる行列`nodes`に対して、それらを順に結んだ線分を描画する関数`draw_segments(nodes)`を定義します。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def draw_segments(nodes):\n", " #nodes[0,:]はx座標の配列、nodes[1,:]はy座標の配列\n", " plt.plot(nodes[0,:],nodes[1,:],\"b\")\n", " plt.grid()\n", " plt.gca().set_aspect(\"equal\")\n", " plt.show()\n", "\n", "nodes = np.array([[0,0],[1,0]]).T #点ごとに入力してから転置すると簡単\n", "draw_segments(nodes) #関数を呼び出して線分を描画する" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "次に、2つの点からなる行列`two_nodes`に対して、その線分に必要な操作を行った結果の点を追加した行列を返す関数`add_nodes(two_nodes)`を定義します。\n", "\n", "ここで、元の行列を\n", "\n", "$$\n", "\\begin{pmatrix}\n", "p_0 & p_1\n", "\\end{pmatrix}\n", "$$\n", "\n", "とすれば、関数が返すべき行列は下図における\n", "\n", "$$\n", "\\begin{pmatrix}\n", "p_0 & p_2 & p_4 & p_3 & p_1\n", "\\end{pmatrix}\n", "$$\n", "\n", "です。\n", "\n", "
\n", "\n", "$p_2$ と $p_3$ は線分 $p_0p_1$ をそれぞれ $1:2$ と $2:1$ に内分する点であるため、\n", "\n", "$$\n", "p_2=\\frac{2}{3}p_0+\\frac{1}{3}p_1,\\quad p_3=\\frac{1}{3}p_0+\\frac{2}{3}p_1\n", "$$\n", "\n", "となり、$p_4$ は $p_2$ を中心に $p_3$ を反時計回りに60度回転させた点であるため、\n", "\n", "$$\n", "p_4=p_2+\n", "\\begin{pmatrix}\n", "\\cos\\theta & -\\sin\\theta\\\\\n", "\\sin\\theta & \\cos\\theta\n", "\\end{pmatrix}\n", "(p_3-p_2),\\quad \\theta=\\frac{60\\pi}{180}\n", "$$\n", "\n", "が成り立ちます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def add_nodes(two_nodes):\n", " p0 = two_nodes[:,0]\n", " p1 = two_nodes[:,1]\n", " p2 = 2/3*p0+1/3*p1 #内分\n", " p3 = 1/3*p0+2/3*p1\n", " theta = 60/180*np.pi\n", " A = np.array([[np.cos(theta),-np.sin(theta)],[np.sin(theta),np.cos(theta)]]) #回転行列\n", " p4 = p2+A@(p3-p2) #p0～p3は一次元配列だが、A@の計算は問題なくできる\n", " return np.array([p0,p2,p4,p3,p1]).T\n", "\n", "two_nodes = np.array([[0,0],[1,0]]).T\n", "nodes = add_nodes(two_nodes)\n", "draw_segments(nodes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "さらに、複数の点からなる行列`nodes`に対して、各線分に`add_nodes`を行った結果の行列を返す関数`update_nodes(nodes)`を定義します。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def update_nodes(nodes):\n", " new_nodes = [nodes[:,0]] #最初の点のみからなるリスト\n", " for i in range(nodes.shape[1]-1): #0からnodesの列数-2まで変化させる\n", " two_nodes = nodes[:,[i,i+1]] #nodesの点を2つずつ使用する\n", " tmp_nodes = add_nodes(two_nodes) #add_nodesを呼び出して行列を返す\n", " for j in range(1,5):\n", " new_nodes.append(tmp_nodes[:,j]) #リストにtmp_nodesの最初以外の点を追加する\n", " return np.array(new_nodes).T\n", " #隣り合う線分は片方の端点を共有しているため、点を重複させないために上のような処理になる\n", "\n", "nodes = np.array([[0,0],[1,0]]).T\n", "nodes = update_nodes(nodes) #1回目の繰り返し\n", "nodes = update_nodes(nodes) #2回目の繰り返し\n", "draw_segments(nodes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "コードをまとめると、次のようになります。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def draw_segments(nodes):\n", " plt.plot(nodes[0,:],nodes[1,:],\"b\")\n", " plt.grid()\n", " plt.gca().set_aspect(\"equal\")\n", " plt.show()\n", "\n", "def add_nodes(two_nodes):\n", " p0 = two_nodes[:,0]\n", " p1 = two_nodes[:,1]\n", " p2 = 2/3*p0+1/3*p1\n", " p3 = 1/3*p0+2/3*p1\n", " theta = 60/180*np.pi\n", " A = np.array([[np.cos(theta),-np.sin(theta)],[np.sin(theta),np.cos(theta)]])\n", " p4 = p2+A@(p3-p2)\n", " return np.array([p0,p2,p4,p3,p1]).T\n", "\n", "def update_nodes(nodes):\n", " new_nodes = [nodes[:,0]]\n", " for i in range(nodes.shape[1]-1):\n", " two_nodes = nodes[:,[i,i+1]]\n", " tmp_nodes = add_nodes(two_nodes)\n", " for j in range(1,5):\n", " new_nodes.append(tmp_nodes[:,j])\n", " return np.array(new_nodes).T\n", "\n", "nodes = np.array([[0,0],[1,0]]).T\n", "n = 3 #繰り返し回数（値を変化させてみるとよい）\n", "for k in range(n):\n", " nodes = update_nodes(nodes)\n", "draw_segments(nodes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 演習1\n", "\n", "正三角形から得られるコッホ曲線を**コッホ雪片（Koch snowflake）**と呼びます。\n", "\n", "上記のコードを利用して、コッホ雪片を描画してください。ただし、繰り返し回数は5回とします。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#演習1のコード\n", "\n", "#正三角形の頂点の座標は各自で決める\n", "#点が時計回りになるようにnodesを定める必要がある" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 演習2\n", "\n", "コッホ曲線を得るための操作における分割点の位置や角度を変えることによって、自分なりの新しいフラクタルを作成して描画してください。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#演習2のコード" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. ミンコフスキーのソーセージ\n", "\n", "**ミンコフスキーのソーセージ（Minkowski sausage）**は、下図のように線分（黒）から3本の線分（赤）を作成する操作を繰り返すことによって得られるフラクタルです。図中の数値は、線分の比を表しています。\n", "\n", "※ミンコフスキーのソーセージにはいくつかのバリエーションがあり、ここではそのうちの一つを紹介します。\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "元の行列\n", "\n", "$$\n", "\\begin{pmatrix}\n", "p_0 & p_1\n", "\\end{pmatrix}\n", "$$\n", "\n", "から新たな行列\n", "\n", "$$\n", "\\begin{pmatrix}\n", "p_0 & p_2 & p_3 & p_1\n", "\\end{pmatrix}\n", "$$\n", "\n", "を求める方法は次の通りです。\n", "\n", "- 線分 $p_0p_1$ 上で $|p_0q|=\\frac{1}{\\sqrt{5}}|p_0p_1|$ を満たす $q$ を求めます。\n", "\n", "- $p_0$ を中心に $q$ を時計回りに $\\theta=\\sin^{-1}\\left(\\frac{1}{\\sqrt{5}}\\right)$ 回転させた点が $p_2$ です。\n", "\n", "- 線分 $p_0p_1$ 上で $|p_1r|=\\frac{1}{\\sqrt{5}}|p_0p_1|$ を満たす $r$ を求めます。\n", "\n", "- $p_1$ を中心に $r$ を時計回りに $\\theta=\\sin^{-1}\\left(\\frac{1}{\\sqrt{5}}\\right)$ 回転させた点が $p_3$ です。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 演習3（オプション）\n", "\n", "正方形から得られるミンコフスキーのソーセージを描画してください。ただし、繰り返し回数は5回とします。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#演習3のコード\n", "\n", "#sinの逆関数は「np.arcsin()」を使えばよい" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 第14回レポート課題\n", "\n", "演習1～演習3に取り組んでください。" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 4 }