対象者

前回の記事の続きです。
実は条件分岐とループ処理が書ければもう大体のことは実現できます。
本記事はいくつか課題を出題しますので、色々調べながら考えてみてくださいね！
初心者以外でも、Pythonまだあまり知らないこと多いな〜もっと知りたいな〜という方はぜひ！
何かいい課題思いついたら追加します。
(こんなのどう？ってアイデアあればぜひ教えてください)

目次

初級編
中級編
上級編

初級編

まずは簡単なものから。
numpyについて調べながらやってみてください。

課題1

$1 \sim 100$までの整数の和を求めるプログラムをできるだけ少ない行数、文字数で書いてください。
(いわゆるコードゴルフというやつです)

課題2

乱数を100個生成して全ての要素を出力してください。

課題3

課題2で生成した乱数の最大・最小・平均・分散・標準偏差を求めてください。

課題4

課題3について、matplotlibを使って箱ひげ図とヒストグラムを描いてください。

初級編解答例

あくまで例です。
もっとスマートなやり方があればぜひ教えてください🙏

解答1

ans1.py

import numpy as np


print(np.sum(np.arange(101)))

解答2

ans2_4.py

import numpy as np
import matplotlib.pyplot as plt


np.set_printoptions(threshold=np.inf)
x = np.random.rand(100)
print(x)

解答3

ans2_4.py

print("max = {}".format(np.max(x)))
print("min = {}".format(np.min(x)))
print("mean = {}".format(np.mean(x)))
print("std = {}".format(np.std(x)))
print("var = {}".format(np.var(x)))

解答4

ans2_4.py

plt.boxplot(x)
plt.show()

plt.hist(x)
plt.show()

中級編

中級編です。ちょっと難し目です。

課題5

ハート型グラフを描いてください。
(ヒント：ハート　グラフ　でググる)

課題6

カージオイドのグラフを描いてください。
ただしカージオイドは以下の数式で表されており、この式中の$a$の値について、$1 \sim 3$までを10等分した値で生成したグラフを全て重ねて表示してください。
($a=1$のグラフを生成、$a=1.222\ldots$のグラフを生成、...)

x = a (1 + \cos\theta) \cos\theta \\
y = a (1 + \cos\theta) \sin\theta \\
0 \le \theta \le 2 \pi \quad 1 \le a \le 3

中級編解答例

解答例です。
だんだん複雑になってきましたね。

解答5

これはググればいろんな方程式が出てきますが、今回は[これ](https://enjoymath.pomb.org/?p=15)を使います。

x^2 + (y - \sqrt{|x|})^2 = 1 \\
-1 \le x \le 1

このままでは描けないので工夫します。

\begin{align}
  (y - \sqrt{|x|})^2 &= 1 - x^2 \\
  y - \sqrt{|x|} &= \pm \sqrt{1 - x^2} \\
  y &= \sqrt{|x|} \pm \sqrt{1 - x^2}
\end{align}

または極形式で

\begin{align}
  x &= \cos \theta \\
  y &= \sin \theta + \sqrt{|x|} \\
    &= \sin \theta + \sqrt{|\cos \theta|}
\end{align} \\
0 \le \theta \le 2 \pi

ans5.py

import numpy as np
import matplotlib.pyplot as plt


def f(x, sign):
    if sign == "plus":
        return np.sqrt(np.fabs(x)) + np.sqrt(1 - x ** 2)
    elif sign == "minus":
        return np.sqrt(np.fabs(x)) - np.sqrt(1 - x ** 2)
    else:
        raise ValueError("不明な符号: {}".format(sign))

def g(theta):
    return np.cos(theta), np.sin(theta) + np.sqrt(np.fabs(np.cos(theta)))


#x = np.linspace(-1, 1, 201)
theta = np.linspace(0, 2 * np.pi, 302)

fig, ax = plt.subplots(1)
#plt.plot(x, f(x, "plus"), color="r")
#plt.plot(x, f(x, "minus"), color="r")
x, y = g(theta)
plt.plot(x, y, color="r")
ax.set_aspect("equal")
plt.show()
plt.savefig("heart.png")

解答6

ans6.py

import numpy as np
import matplotlib.pyplot as plt

def r(a, theta):
    coef = a * (1 + np.cos(theta))
    return coef * np.cos(theta), coef * np.sin(theta)

theta = np.linspace(0, 2 * np.pi, 100)
a = np.linspace(1, 3, 10)

fig, ax = plt.subplots(1)
for i in a:
    x, y = r(i, theta)
    plt.plot(x, y)
ax.set_aspect("equal")
plt.show()
plt.savefig("cardioid.png")

上級編

さらに複雑になります。
色々ググりながらやってみよう。

課題7

サイクロイドのwikiにあるアニメーションを描いてください。
サイクロイド曲線だけでOKですが、可能なら動円と動径も描きましょう。

課題8

ある自然数$n \gt 2$を入力すると、その数までの素数のリスト(またはnumpy配列)を返す関数を作成してください。
この時、$n = 1e6 = 100000$を入力した時の実行時間が0.1秒を下回るようにしてください。

上級編解答例

解答7

サイクロイドは以下の設定で描きます。

x = r(\theta - \sin \theta) \\
y = r(1 - \cos \theta) \\
r = 1 \quad 0 \le \theta \le 2 \pi

オプションとして、動円は

(x - \theta)^2 + (y - r)^2 = r^2 \\
\Leftrightarrow \left\{
  \begin{align}
    x &= r \cos \psi + \theta \\
    y &= r \sin \psi + r
  \end{align}
\right. \\
0 \le \psi \le 2 \pi

動径はサイクロイド点を$(a, b)$、動円の中心を$(c, d)$とすると

\left\{
  \begin{array}{cc}
    y - d = \cfrac{d - b}{c - a} (x - c) & (x \lt |c - a|, a \ne c) \\
    x = c & (a = c, b \le y \le d \quad or \quad d \le y \le b)
  \end{array}
\right.

のように描けます。

ans7.py

%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation


def f(r, theta):
    return r * (theta - np.sin(theta)), r * (1 - np.cos(theta))

def circle(r, theta, n_point):
    psi = np.linspace(0, 2 * np.pi, n_point)
    return r * np.cos(psi) + theta, r * np.sin(psi) + r

def line(r, theta, n_point):
    a, b = f(r, theta)
    if theta - a == 0:
        x = np.full(n_point, theta)
        y = np.linspace(min(b, r), max(b, r), n_point)
        return x, y
    else:
        x = np.linspace(min(a, theta), max(a, theta), n_point)
        y = (r - b) / (theta - a) * (x - theta) + r
        return x, y


# パラメータ設定
n_point = 200
r = 1
theta = np.linspace(0, 2 * np.pi, n_point)

# グラフ設定など
fig = plt.figure()
ax = plt.axes()
plt.xlim(theta[0] - 0.5, theta[-1] + 0.5)
plt.ylim(0, 2.5 * r)
plt.grid()
interval = 5

# シーン作成
images = []
for i in range(len(theta)):
    # 動円
    x, y = circle(r, theta[i], n_point)
    im1, = plt.plot(x, y, color="g")
    # 動径
    x, y = line(r, theta[i], n_point)
    im2, = plt.plot(x, y, color="r")
    # サイクロイド曲線
    x, y = f(r, theta[:i])
    im3, = plt.plot(x, y, color="b")
    # サイクロイド点
    if i != 0:
        im4, = plt.plot(x[-1], y[-1], color="b", marker="o")
    else:
        im4, = plt.plot(x, y, color="b", marker="o")
    ax.set_aspect("equal")
    images.append([im1, im2, im3, im4])

# グラフを描く
ani = animation.ArtistAnimation(fig, images, interval=interval)
plt.show()
#ani.save("cycloid.gif", writer="imagemagick")

jupyter notebook上でgifを表示するには

ans7.py

%matplotlib nbagg

を最初に記述する必要があります。 ![cycloid.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/640911/64aaabcc-4afe-e609-7caf-fba813f7d663.gif) 同じ要領で[外サイクロイド](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A1%E3%82%A4%E3%83%AB:Epicycloid(5,2)_animated.gif)や[内サイクロイド](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A1%E3%82%A4%E3%83%AB:Hypocycloid(5,2)_animated.gif)なんかにチャレンジしてもいいですね！ **<<注意>>** [`matplotlib`のバグ](https://github.com/matplotlib/matplotlib/issues/17097)でgif保存できません。

上記のgifはスクリーンショット→動画をgifにするアプリを用いて生成したgifです...

解答8

ここでは**エラトステネスのふるい**というアルゴリズムを用います。アルゴリズムの内容自体は至極単純で、

フラグを$n + 1$個用意する(インデックス番号と自然数を対応させるため)
$0$番目と$1$番目のフラグは折る($=$自然数$0$と$1$は素数ではない)
$i = 2 \sim \sqrt{n} + 1$までループ
- $i$番目のフラグが立っていれば素数なので、それ以外の$i$で割り切れる数のフラグを折る
フラグが立っているインデックスを返す

という感じです。

ans8.py

import numpy as np
import time as time

np.set_printoptions(threshold=10, edgeitems=10)

# エラトステネスのふるい
def prime(n):
    numbers = np.ones(n + 1, dtype=np.bool)
    numbers[0 : 2] = False
    border = int(np.sqrt(len(numbers))) + 1
    for i in range(2, border):
        if numbers[i]:
            numbers[::i] = False
            numbers[i] = True
    return np.where(numbers)[0]

start = time.time()
primes = prime(int(1e6))
end = time.time()
print(primes)
print("time: {}".format(end - start))

参考

matplotlib アニメーション (Jupyter)

初めてのPythonシリーズ記事