SymPy による数式処理とグラフ作成 - 弘前大学 Home Sweet Home

相关文章推荐

聪明伶俐的皮带 · wxwidgets - C++ ...· 2 年前 ·

酷酷的蜡烛 · 广义线性混合模型GLMM_FarmerJoh ...· 3 年前 ·

知识渊博的跑步鞋 · 使用 LibreOffice ...· 3 年前 ·

好帅的饺子 · 服务器上的 XML | 菜鸟教程· 3 年前 ·

葛西真寿　弘前大学大学院理工学研究科

この Notebook では，「 wxMaxima による数式処理とグラフ作成」のテキストの内容を，Jupyter Notebook 上で Python の SymPy ライブラリを使って説明しています。

セクションの構成は，wxMaxima 版のテキストに準じています。

ですからこの Notebook は，一から SymPy を始めようとしている人だけでなく，Maxima は知っているが同じことを SymPy ではどのように計算するかということに興味を持っている人にも参考になるのではないかと思います。

Jupyter Notebook （の SymPy）で計算した結果を保存するには，「File」メニューから「Save as ... 」でファイル名をつけて保存します。

「File」メニューの真下に，なにやら四角で右上隅がかけている形のアイコンがありますが（これが「フロッピーディスク」というものであることは，20世紀の昔を知る人にしかわからないと思います），このアイコンをクリックしても保存されます。

また，「File」メニューから「Open... 」を選び，計算結果が保存された Notebook をまた利用することもできます。

（練習） ピタゴラス数

$a^2 + b^2 = c^2$ をみたす正の整数の組をピタゴラス数といいます。

ピタゴラス数は2つの任意の正の整数 $m, n$ から以下のようにして何組でもつくることができます。

$$ a = |m^2 - n^2|, \quad b=2mn, \quad c = \sqrt{a^2 + b^2} $$

このとき，$c$ は（平方根をとるのにもかかわらず）必ず整数になることを示しなさい。

（ヒント）

1等星と6等星の見かけの明るさをそれぞれ，$L_1,\ L_6$ とすると，

$$1 = -2.5\log_{10} L_1 + \mbox{定数}, \quad 6 = -2.5\log_{10} L_6 + \mbox{定数} $$

両辺を引くと，

$$ 1 - 6 = -2.5 (\log_{10} L_1 - \log_{10} L_6) $$

つまり，

$$ \log_{10} \frac{L_1}{L_6} = \frac{-5}{-2.5} = 2 $$

（練習） 屋根勾配

日本の大工さんは，家の屋根の勾配は角度ではなく，伝統的に3寸勾配，4寸勾配のように水平方向に1尺（10寸）に対して垂直方向に何寸立ち上がるかで示します。

（本当です。実際に家を建てた本人が言っているのだから間違いないです。）

下図のログハウスの屋根は8寸勾配です。傾斜角度にすると何度ですか？

（大工さんは8寸勾配と言いますが，照明器具を設置する電気屋さんは8寸勾配ではなく，傾斜角は何度だ？と聞きますので，この変換が必要でした。）

次の積分を計算しなさい。

$$\mbox{1.}\ \int e^x \sin x\ dx, \qquad\mbox{2.}\ \int \frac{x^2-2x-2}{x^3-1}\, dx

（練習） 鶴と亀と蟻

鶴と亀と蟻，個体数の合計は10。足の数は全部で34本。蟻は亀より1匹少ない。鶴，亀，蟻はそれぞれ何羽・何匹？

鶴・亀・蟻の個体数をそれぞれ $x, y, z$ として連立方程式を立てて solve します。蟻の足は何本かわかりますよね？

かつて線形代数の授業でこの問題を出したとき，真っ先に出た質問が「先生，蟻の足は何本ですか？」でした... orz

y = Function('y')(x)
a = Symbol('a')
eq = Eq(Derivative(y, x), a*y)
x = Function('x')(t)
K = Symbol('K', positive = True)
eq = Eq(Derivative(x, t, 2), -K*x)
ans = dsolve(eq, x,
      ics = {x.subs(t, 0):x0,
             diff(x, t).subs(t, 0):v0}
x = Function('x')(t)
dsolve(Derivative(x, t, 2), x,
       ics = {x.subs(t, 0):x0,
              diff(x, t).subs(t, 0):v0}
（練習） 一様重力場中の投げ上げ運動
時刻 $t=0$ に高さ $0$ から初速度 $v_0$ で鉛直上方に投げ上げた物体の時刻 $t$ における高さ $y(t)$ を求めなさい。運動方程式は以下の通り。
$$\frac{d^2 y}{dt^2} = - g$$
ただし，$g$ は重力加速度の大きさで，定数である。
速度に比例する空気抵抗がある場合，運動方程式は以下のようになる。
$$\frac{d^2 y}{dt^2} = - g - \beta\frac{dt}{dt}$$
ここで，$\beta$ は定数。前問と同じ初期条件のときに $y(t)$ を求めなさい。
前問2.で求めた解について， $\beta \rightarrow 0$ のときに前問1. の答えに一致するかどうか，確かめなさい。
a1, a2, a3 = symbols('a1 a2 a3', real=True)
a = Matrix([a1, a2, a3])
（練習） ベクトルの内積，外積
$\displaystyle \boldsymbol{a} =
\left(-\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, \sqrt{3}\right), \quad
\boldsymbol{b} =
\left(-1, 1, -\sqrt{2}\right)$ について以下の量を計算しなさい。
$$1.\ |\boldsymbol{a}|, \quad 2.\ |\boldsymbol{b}|, \quad
3.\ \boldsymbol{a}\cdot \boldsymbol{b}, \quad
4.\ \boldsymbol{a}\times \boldsymbol{b} $$また，2つのベクトル $\boldsymbol{a}, \boldsymbol{b}$ のなす角を $\theta$ としたときの $\cos\theta$ および $\sin\theta$ を求めなさい。
スカラー３重積の以下の性質を確かめなさい。
$$\boldsymbol{a}\cdot (\boldsymbol{b}\times\boldsymbol{c}) =
\boldsymbol{b}\cdot (\boldsymbol{c}\times\boldsymbol{a}) =
\boldsymbol{c}\cdot (\boldsymbol{a}\times\boldsymbol{b}) $$
ベクトル３重積の以下の性質を確かめなさい。
$$\boldsymbol{a}\times (\boldsymbol{b}\times \boldsymbol{c}) =
(\boldsymbol{a}\cdot \boldsymbol{c})\boldsymbol{b} -
(\boldsymbol{a}\cdot \boldsymbol{b}) \boldsymbol{c}
O = Matrix([[cos(theta), -sin(theta), 0],
            [sin(theta),  cos(theta), 0],
            [0,           0,          1]])
（練習） 直交変換されたベクトルの大きさ
直交変換されたベクトル $\boldsymbol{r}^{\prime} = O \boldsymbol{r}$ の大きさは，もとのベクトル $\boldsymbol{r}$ の大きさと同じであることを示しなさい。
ヒント：例えば以下のように...
plot(cos(x), (x, -2*pi, 2*pi),
     axis_center=(-2*pi, -1),
     line_color='red',
     legend=True,
     xlabel='x',
     ylabel='y');
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] =8,6
plt.rcParams['font.size'] = 12
from IPython.display import set_matplotlib_formats
%matplotlib inline
set_matplotlib_formats('svg')
p= plot(x**2 - 1, 4*x - 5, (x, -5, 5),
        ylim=(-5, 10), legend=True, show=False);
p[0].line_color='b'
p[1].line_color='r'
p.show()
    return 4*x - 5
# f(x) = g(x) をみたす x を求める。
sol = solve(Eq(f(x), g(x)), x)
Python では，あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。
以下のような内容のファイル test.dat がカレント・ディレクトリ ./ にあるとします。
# x   y
  0   0
  1   1
  2   4
  3   9
  4   16
  5   25
# test.dat の x1, y1 は前節で定義済み。
plt.plot(x1, y1, 'bo', label='data');   # test.dat の x1, y1 を青丸で plot
# y = x**2 のデータをここで作成。
x = np.arange(0, 5.2, 0.2)
y = x**2
plt.plot(x,  y,  '-r', label='$y=x^2$'); # y = x**2 を赤線で plot
plt.legend();                            # label で定義しておいた判例を表示
半径1の円の方程式は $x^2 + y^2 = 1$ です。
この円を描くには，$y = \pm\sqrt{1-x^2}$ として $y = f(x)$ の形にするよりも，以下のような媒介変数表示にしたほうが簡単でしょう。
$$ x = \cos t, \quad y = \sin t, \quad(0 \le x \le 2\pi)$$このような媒介変数表示の2次元グラフを SymPy で描くには以下のようにします。
import matplotlib.pyplot as plt
# デフォルトでは plt.rcParams['figure.figsize'] = (6.0, 4.0) だった
plt.rcParams['figure.figsize'] = (5, 5)
plot_parametric(cos(t), sin(t), (t, 0, 2*pi), xlim=(-1.1, 1.1), ylim=(-1.1, 1.1));
（練習） 楕円のグラフ（楕円中心を原点として）
同様にして，楕円のグラフを描くこともできます。原点を中心とし，長半径 $a$，短半径 $b$ の楕円の式は
$$ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 $$です。媒介変数表示では，
$$ x = a \cos t, \quad y = b \sin t \quad (0 \le t \le 2\pi)$$と書けます。$a, b$ に適当な値を入れて楕円のグラフを描きなさい。
太陽からの万有引力を受けて運動する惑星（惑星だけでなく，準惑星や小天体も含みます）は，太陽を焦点とした楕円軌道を描きます。焦点を原点とし，長半径 $a$，離心率 $e$ の楕円の方程式は，極座標 $r, \phi$ を使って以下のように表すことができます。
$$ r = \frac{a(1-e^2)}{1 + e\cos \phi}$$さて，かつては第9惑星，現在では準惑星の一つである冥王星も楕円軌道を描きます。冥王星の軌道長半径 $a_P = 39.767 \,\mbox{au}$，離心率 $e_P = 0.254$ を使って冥王星の軌道を描きます。
まず，極座標表示の楕円の式を関数として定義します。
import matplotlib.pyplot as plt
# デフォルトでは plt.rcParams['figure.figsize'] = (6.0, 4.0) だった
plt.rcParams['figure.figsize'] = (6, 6)
aP = 39.767
eP = 0.254
plot_parametric(r(aP,eP,phi)*cos(phi), r(aP,eP,phi)*sin(phi),
                (phi, 0, 2*pi), xlim=(-50,50), ylim=(-50,50));
aN = 30.1104
eN = 0
p = plot_parametric((r(aP,eP,phi)*cos(phi), r(aP,eP,phi)*sin(phi)),
                    (r(aN,eN,phi)*cos(phi), r(aN,eN,phi)*sin(phi)),
                    (phi, 0, 2*pi), xlim=(-50, 50), ylim=(-50, 50), show=False)
p[0].line_color='b'
p[1].line_color='r'
p.show()
plt.rcParams['figure.figsize'] = (6, 4)
p = plot_parametric((r(aP,eP,phi)*cos(phi), r(aP,eP,phi)*sin(phi)),
                    (r(aN,eN,phi)*cos(phi), r(aN,eN,phi)*sin(phi)),
                    (phi, 0, 2*pi), xlim=(25, 35), ylim=(5, 20),
                    axis_center=(25, 5), show=False)
p[0].line_color='b'
p[1].line_color='r'
p.show()
plt.rcParams['figure.figsize'] = (6, 6)
p = plot_parametric((r(aP,eP,phi)*cos(phi), r(aP,eP,phi)*sin(phi)),
                    (r(aN,eN,phi)*cos(phi), r(aN,eN,phi)*sin(phi)),
                    (aN*phi/(2*pi)*cos(ans1), aN*phi/(2*pi)*sin(ans1)),
                    (aN*phi/(2*pi)*cos(ans2), aN*phi/(2*pi)*sin(ans2)),
                    (phi, 0, 2*pi), xlim=(-50, 50), ylim=(-50, 50),
                    axis_center=(0, 0), show=False)
p[0].line_color='b'
p[1].line_color='r'
p[2].line_color='black'
p[3].line_color='black'
p.show()
グラフをみると，月別平年気温は12ヶ月周期の正弦関数または余弦関数のように見えます。では，以下のように関数フィットをしてみましょう。
関数フィットするために，以下のように scipy.optimize を import します。以下のような関数でフィットしてみます。
$$f(x) = \beta_0 + \beta_1  \cos\left(\frac{2\pi x}{12}\right) + \beta_2 \sin\left(\frac{2\pi x}{12}\right)$$
def theoreticalValue(beta):  # sympy.cos ではなく，np.cos や np.pi を使う。
    f = beta[0] + beta[1] * np.cos(2*np.pi*x / 12) + beta[2] * np.sin(2*np.pi*x / 12)
    return f
x = np.array(month)
y = np.array(temp)
initialValue = np.array([10, 10, 1])
betaID = leastsq(objectiveFunction, initialValue)
betaID[0]
def ftemp(month):
    r = beta0 + beta1 * np.cos(2*np.pi*month/12) + beta2 * np.sin(2*np.pi*month/12)
    return r
month1 = np.arange(1, 12.1, 0.1)
temp1 = ftemp(month1)
plt.plot(x, y, 'r.')
plt.plot(month1, temp1, 'b');