.. toctree:: ===================== 10. numpyの基本と応用 ===================== numpyは、科学技術計算を容易にかつ高速に実現するために開発されたライブ ラリです。これまでに学んだ知識を使ってリストや多次元配列にデータをため て、for文等で配列のデータにアクセスできますが、データが巨大になると圧 倒的に遅くなります。numpyを使うことで、多次元配列の処理を簡潔に書けて 高速に処理できるようになります。またnumpyは、数学関数や乱数発生、線形 代数、統計処理、データの入出力等、科学技術計算で基本となる非常に豊富な 機能を提供しています。 https://docs.scipy.org/doc/numpy-1.13.0/reference/ 今回の学習目標は以下のとおりです。 * numpyでの配列操作(indexing, broadcasting)を理解する。 * numpyで使える代表的な関数を理解する。 * numpyで線形代数計算、統計処理を行う理解する。 なお慣習にしたがって、numpytは別名npとしてimportします。 .. code-block:: python :caption: example-01 import numpy as np # npとしてnumpyをimport (慣例) 多次元配列とデータ型 ==================== リスト変数をnp.array()を使って配列(array)オブジェクトを生成します。多 次元リストのようにサイズの異なる要素になることはできません。np.array() で生成したインスタンスは、要素数、形(行と列数)、型等は、クラス変数 (size, shpae, dtype等)として保持しています。 https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.ndarray.html .. code-block:: python :caption: example-02 import numpy as np # 変数(スカラー) s = 7 # 1次元配列 (全部整数なのでint64) (ベクトル) one = np.array([3, 4, 5]) # 2次元配列 = 行列 (1.1を1つ含むのでfloat64) 3(行)x3(列) two = np.array([[1.1, 2, 3], [5, 6, 7], [5, 6, 7]]) print(one.size, one.shape, one.dtype) # 3 (3,) int64 print(two.size, two.shape, two.dtype) # 9 (3, 3) float64 元素記号と座標のように、文字とfloatが混在したデータを保持したい場合は、 dtypeで型をobjectに指定します。またastype()メソッドで、生成したインス タンスのデータの型を後で変更することができます。objectタイプ場合、一部 のデータをfloat型へ変換することももできます。 .. literalinclude:: example-03.py :caption: example-03 .. note:: 複数の型を一つのarrayで取り扱うために、structured data-typeを定義する 方法もありますが、やや取扱いが面倒なのでここでは省略します。複雑な複 数の型のデータ操作を行う場合は、Panda (https://pandas.pydata.org/) モ ジュールを使った方が適切かもしれません。 numpyで指定できる正確な型の指定は以下のURLを参照してください。 https://docs.scipy.org/doc/numpy-1.13.0/user/basics.types.html numpyで生成した配列は、forやwhileによる繰り返し処理を行うこともできま す。ただし、forによる処理は非常に低速です。高速に計算させるためには、 forやwhileを使わずにnumpyで提供されている機能だけでプログラムを作成す る方が望ましいです。 .. code-block:: python :caption: example-04 :linenos: import numpy as np data = np.array([[1.1, 2, 3], [5, 6, 7], [5, 6, 7]]) # forでデータにアクセス for d in data: print(d[0], d[1], d[2]) 配列の生成と変形 ================ 配列を生成する方法と変形する方法を説明します。いくつか代表的なものの説 明に留めます。 https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.array-creation.html .. list-table:: 配列生成関数の一覧 :header-rows: 1 :widths: 20 40 40 :stub-columns: 0 * - 関数 - 引数 - 機能 * - np.linspace() - start, end, N - startからendまでをN分割した配列生成 (defaultでendは含む) * - np.arange() - start, end, step - startからendまでをstepで刻んだ配列生成 (endは含まない) * - np.ones() - (N, N, N...) - 多次元の要素が1の配列を生成 * - np.zeros() - (N, N, N...) - 要素が全て0の多次元の配列を生成 * - np.empy() - (N, N, N...) - 多次元の空配列を生成 * - np.random.random() - (N, N, N...) - 多次元のランダム[0, 1)の要素をもつ配列を生成 いくつか例を交えて説明します。 .. literalinclude:: example-04.py :caption: example-05 reshape()で、配列の形を自由に変えることができます。reshape()メソッドは、 適用する元のデータを変更しません。変更した配列は、返り値として受け取れ ます。reshapeは、メソッドとして用いても同じ結果が得られます。読みやす いように記述すれば良いでしょう。なおnumpyで使えるたいていの関数は、メ ソッドとしても利用できます。 .. literalinclude:: example-05.py :caption: example-06 ファイルからデータのload ======================== スペースで区切られたxとyのデータのように、データの並びが決まっていれば、 np.loadtxt()でファイルからデータを読んで配列を生成することがきます。デ フォルトで、np.loadtgxtha #から始まる行はスキップされ、スキップでデー タが区切られているととしてデータをloadします。デフォルトの挙動を変えた い時は、(skiprowsによる先頭行のスキップ、usecolsによる特定列の読み込み 等)マニュアルを参照してください。 https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.loadtxt.html 以下のような内容の :download:`test.dat <./test.dat>` をnp.loadtxt()で読み込んでみます。 .. literalinclude:: test.dat 2行文はskiprowsで読み飛ばします。#から始まる行は、自動でスキップされfloat型としてデータが読み込まれています。 .. literalinclude:: example-06.py :caption: example-07 実はnp.loadtxt()は非常に遅いです。文字として読んでnp.array()で配列に変換した方 が早いです。比較するプログラムは、以下のとおりです。 .. literalinclude:: example-07.py :caption: example-08 datetimeモジュールのdatetime.now()で現在時刻を取得し、計算時間を比較しています。 データの保存は、np.savetxt()を行うことができますが、ここでは取扱いませ ん。 Indexing ======== Indexingとは、 エレガントに配列をスライスしたり条件付きでデータを選択 できるnumpyの機能です。自由自在にIndexing操作ができるようになれば、 numpyマスターに一歩近づくことができるでしょう。詳細は以下のURLを参照し てください。 https://docs.scipy.org/doc/numpy-1.13.0/reference/arrays.indexing.html まず1次元の配列について例題を交えながら説明します。特に明示しない限り、 Indexingは浅いコピー(viewと呼ばれる)になります。新しい配列を深いコピー で生成したい場合、np.copy()メソッドでコピーを生成します。また要素番号を配列を渡して スライスしたり、TrueとFalseによるデータの選択は深いコピーになります。 配列に条件式を適用すると、同じサイズとTrueとFalseの配列が得られます。 TrueとFalseを使って行と列の選択を深いコピーで行うことができます。 .. literalinclude:: example-08.py :caption: example-09 2次元配列の例です。 **行**, **列** の順番でスライスする数字を指定する ことに注意してください。ベクトルでデータのスライスを指定すると深いコピー になります。またクロスした行と列の選択は、np.ix_()を使って選択できます。 こちらも深いコピーです。 .. literalinclude:: example-09.py :caption: example-10 3次元の場合、同じサイズの行列がいくつも重なっているイメージをもてば良 いでしょう。2次元の考え方に慣れればは同じです。 .. literalinclude:: example-11.py :caption: example-11 4次元以上に次元が増えても、3次元の考え方と同じです。 .. 講義後に追加 ============ 大窪がメモしたファイルを置いておきます。参考にしてください。 :download:`第12回講義メモ <2019-01-09/第12回講義メモ.pdf>` .. アンケートの統計 ---------------- .. image:: 2019-01-09/statistics.png :width: 800px :align: center ブロードキャスティング ====================== サイズの同じリスト(配列)に演算子を作用させると、リスト同士の演算できま す。サイズの異なるサイズの場合、小さい方の配列サイズを調整して演算して くれます。このような機能をブロードキャスティングと呼びます。。 https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html .. literalinclude:: example-12.py :caption: example-12 ブロードキャスティングを使うことで、リストの演算を直感的に理解しやすく なりますし、forでぐるぐる回すより圧倒的に高速です。numpyの配列による計 算とリスト変数での計算で、どれくらい計算時間が変わるかdatetimeモジュー ルを使って時間を計測して比較してみます。 .. literalinclude:: example-13.py :caption: example-13 計算結果は以下のとおりです(大窪のパソコンの場合)。numpy配列を使った方 が100倍以上早いです。 .. code-block:: none リスト変数: 0:00:00.324808 numpy配列: 0:00:00.003042 リスト変数: 0:00:00.025466 numpy配列: 0:00:00.000072 indexingとbroadcastingを組み合わせれば、柔軟なデータの操作や演算、条件 を満たすデータへのアクセスが自由自在にできます。indexingとbroadcasting による演算の例です。 .. literalinclude:: indexing-broadcasting.py :caption: example-14 メタンのxyzデータを使ってIndexingとbroadcastingの応用を考えてみます。 元素記号と座標のデータを一緒にした配列を作成した方が、データの取扱いも プログラムの見やすさも良くなりますので、object型でデータを保持してプログラムを作成してみます。 .. literalinclude:: CH4-indexing-broadcasting.py :caption: example-15 結合、繰り返し ============== 2つの配列を結合(np.concatenate())したり、1つの配列データを繰り返すここ とで新しい配列を生成(np.repeat()とnp.tile())することもできます。以下の URLにたくさんありますが、使いそうな例だけ示します。多次元配列の場合、 軸を指定することで、結合や繰り返し方向を指定することができます。 https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.array-manipulation.html .. literalinclude:: example-14.py :caption: example-16 数学関数 ======== 数学関数の一部を紹介します。引数に配列aを渡せば、配列aの全ての要素に関数が適用されます。 https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.math.html .. list-table:: 数学関数一覧 :header-rows: 1 :widths: 20 40 :stub-columns: 0 * - 関数 - 機能 * - np.sin(a) - :math:`\sin(\theta)` の計算。xはラジアンで与えることに注意。 * - np.cos(a) - :math:`\cos(\theta)` の計算。xはラジアンで与えることに注意。 * - np.tan(a) - :math:`\tan(\theta)` の計算。xはラジアンで与えることに注意。 * - np.arccos(a) - :math:`\cos^{-1}(x)` の計算。返り値はラジアンになることに注意。 * - np.arcsin(a) - :math:`\sin^{-1}(x)` の計算。返り値はラジアンになることに注意。 * - np.arctan(a) - :math:`\tan^{-1}(x)` の計算。返り値はラジアンになることに注意。 * - np.log(a) - 自然対数 :math:`\log(x)` の計算。 * - np.exp(a) - 指数対数 :math:`\exp(x)` の計算。 * - np.sqrt(a) - 平方根 :math:`\sqrt{x}` の計算。 * - np.log10(a) - 常用対数(底が10) :math:`\log(x)` の計算。 * - np.round(a, num) - aをnumの桁数で丸めます。 * - np.floor(a) - 切り捨て(床)を行います。 * - np.ceil(a) - 切り上げ(天井)を行います。 データを可視化するために、 :doc:`matplotlib <../12_Matplotlibの使い方/12>` を使って、データをプロットしてみます。 .. literalinclude:: example-15.py :caption: example-17 プログラムを実行すると次のような図が得られます。 .. image:: plot.png :width: 500px :align: center 注意点としてobject型の場合、numpyの数学関数に渡すためにastype()でfloat 型に変換せねばなりません。以下の例をみてください。np.sqrt()にobjectを 渡すとエラーになります。 .. literalinclude:: CH4-math.py :caption: example-18 .. 講義後に追加 ============ 大窪がメモしたファイルを置いておきます。参考にしてください。 :download:`第14回講義メモ <./第14回講義メモ.pdf>` .. .. image:: 2018-01-17/statistics.png :width: 800px :align: center 統計関数、ソートと検索 ====================== 合計や平均の計算等、統計関数の一部と、配列のソートと検索を行う関数を紹 介します。リスト変数で使えるsorted()関数やindex()メソッドよりも高速かつ柔軟な関 数がnumpyに準備されています。多次元の配列の場合、軸を選択した処理を行 うことができます。 統計処理を行う関数のマニュアルは、以下のURLにリストされていますが、い くつか代表的なものをピックアップして紹介します。 https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.statistics.html .. list-table:: 統計関数 :header-rows: 1 :widths: 20 40 :stub-columns: 0 * - 関数 - 機能 * - sum(a[, axis]) - 軸にそった足し算を行う。 * - mean(a[, axis]) - 軸にそった平均値を返す * - std(a[, axis]) - 軸にそった標準偏差を返す * - histogram(a, bins=n or ()) - ヒストグラムを生成します。binsが整数ならデータ範囲を等しくn分割 した区分とします。binsにリストを渡せば区分のエッジを指定するこ とになります。 統計関数を使う計算例を示します。ヒストグラムは、乱数で適当な数字を生成 してプロットしています。 .. literalinclude:: example-16.py :caption: example-19 ヒストグラムのプロットで得られる結果は以下のとおりです。 .. image:: histogram.png :width: 500px :align: center ソートや検索を行う処理のマニュアルは、以下のURLになります。いくつか代 表的なものをピックアップして紹介します。 https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.sort.html .. list-table:: ソート、検索、カウント関数 :header-rows: 1 :widths: 20 40 :stub-columns: 0 * - 関数 - 機能 * - sort(a[, axis]) - 軸にそってsortした結果を返す * - argsort(a[, axis]) - 軸にそってsortした要素番号の結果を返す * - max(a[, axis]) - 軸にそった最大値の要素を返す * - argmax(a[, axis]) - 軸にそった最大値の要素番号を返す * - min(a[, axis]) - 軸にそった最小値の要素を返す * - argmin(a[, axis]) - 軸にそった最小値の要素番号を返す ソートと検索を行う例を示します。 .. literalinclude:: example-17.py :caption: example-20 条件を満たす要素の数をカウントしたい場合、True/Falseの配列のsumを求め ればカウントできます(True:1, False:0)。 .. literalinclude:: example-18.py :caption: example-21 線形代数計算 ============ 線形代数で学んだ行列の計算を行う関数をピックアップして紹介します。ここ では、numpyの配列を行列と読み替えて考えてください。 https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.linalg.html .. list-table:: :header-rows: 1 :widths: 20 40 :stub-columns: 0 * - 関数 - 機能 * - np.dot(a, b) - 行列aと行列bの行列積の計算 @で代用可能 * - np.linalg.det(a) - 行列式の計算 * - np.linalg.inv(a) - 行列aの逆行列の計算 * - np.transpose(a) - 行列aの転置の計算 .Tでも同じ * - np.identify(n) - n次の単位行列を生成する。 * - np.linalg.eig(a) - aの固有値と固有ベクトルを返す。固有値(w)と、固有値に対応した固有ベクトル(v)(長さで規格化されている) * - np.linalg.solve(a, b) - 係数行列aとベクトルbを渡すとxを返す。 ax = b 行列積の計算は、行列の形に気をつけてください。A(m,n)*B(n,l)=C(m,l)のよ うにサイズが代わります。ここで、Aの列(n)とBの行(n)が揃ってなければ計算 できません。 .. literalinclude:: example-19.py :caption: example-22 具体的に数学の問題をいくつか解いてみます。np.solve()を使って以下の連 立方程式を解いてみます。 .. math:: 3x + 2y + 8z = 3 \\ 2y = 9z \\ -10x + y + 4 = 9z まず連立一次方程式を行列で書き直します。 .. math:: \begin{bmatrix} 3\\ 0\\ -4 \end{bmatrix} = \begin{bmatrix} 3 & 2 & 8\\ 0 & 2 & -9\\ -10 & 1 & -9\\ \end{bmatrix} \begin{bmatrix} x\\ y\\ z \end{bmatrix} 3つの連立一次方程式程度の計算は、一瞬でといてくれます。 .. literalinclude:: solve.py :caption: example-23 回転行列を使って分子を回転させる操作を行ってみましょう。細かい回転行列 はwikipediaでも参照してください。線形代数でもでてきたかと思います。任 意の座標を軸回りに :math:`\theta` 回転させることができます。 https://ja.wikipedia.org/wiki/%E5%9B%9E%E8%BB%A2%E8%A1%8C%E5%88%97 3次元の場合、 :math:`x, y, z` の回転行列 :math:`Rx, Ry, Rz` は以下のとおりです。 .. math:: Rx(\theta) = \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos\theta & -\sin\theta\\ 0 & \sin\theta & \cos\theta\\ \end{bmatrix} \quad Ry(\theta) = \begin{bmatrix} \cos\theta & 0 & \sin\theta\\ 0 & 1 & 0\\ -\sin\theta & 0 & \cos\theta\\ \end{bmatrix} \quad Rz(\theta) = \begin{bmatrix} \cos\theta & -\sin\theta & 0\\ \sin\theta & \cos\theta & 0\\ 0 & 0 & 1\\ \end{bmatrix} \quad もとの座標 :math:`(x, y, z)` を :math:`\theta` だけ回転させた座標 :math:`(x', y' ,z')` は、次のように計算できます(:math:`x` 軸回りに回転させる場合)。 .. math:: \begin{bmatrix} x'\\ y'\\ z'\\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos\theta & -\sin\theta\\ 0 & \sin\theta & \cos\theta\\ \end{bmatrix} \begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} CH4をx方向に10度ずつ回しながらx方向に2Ang.シフトさせたxyzファイルを10個 生成してみます。 .. literalinclude:: CH4-roatate.py :caption: example-24 プログラムを実行すると :download:`CH4-rotate.xyz <./CH4-rotate.xyz>` が生成されます。作成したxyzファイルをvestaで開いた結果を下にしめします。x軸に沿ってぐるぐる回っていることがわかるかと思います。 .. image:: CH4-rotate.png :width: 400px :align: center .. 講義後に追加 ============ 大窪がメモしたファイルを置いておきます。参考にしてください。 :download:`第15回講義メモ <./第15回講義メモ.pdf>` .. image:: 2018-01-24/statistics.png :width: 800px :align: center 多項式近似、内挿 ================ 実験データを下式のような多項式でフィッテングすることはよくあります。こ こで、 :math:`n` はdegree(次数)で自分でモデルを決めてフィッテングしま す。np.polyfit()を使えば、簡単に多項式フィットを行うことができます。 https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.polyfit.html .. math:: f(x) = a_0 + a_1x + a_2x^2 + \cdots + a_nx^n = \sum_{i=0}^n a_i x^i 適当な多項式の関数にノイズを加えたデータを作成し、多項式でフィッテング してみます。xとyがデータ、degを次数として、np.polyfit(x, y, degree)を 呼び出せば :math:`a_i` を戻り値として求めてくれます。np.polyfit()の戻 り値とxfをnp.polyvalに渡せば、多項式フィットしたyfを得ることができます。 以下の多項式でテストデータを作成し、np.random.random()でノイズを印加し て、1, 2, 3次でフィッテングしたプログラムを例に示します。 .. math:: &y = 0.1x^3 - 0.3x^2 - 10x + 5\\ &(a_3, a_2, a_1, a_0) = (0.1, -0.3, 10, 5) .. literalinclude:: example-20.py :caption: example-25 3次の多項式フィットで求めた結果は以下のとおりになります。ノイズを印加しているので、元の関数の係数を完全に再現できていませんが近い結果を得ています。 .. math:: (a_3, a_2, a_1, a_0) = (0.09852041, -0.30923446, -10.01156373, 10.70003027) 結果をプロットした図は以下のとおりです。 .. image:: polynominal.png :width: 500px :align: center データの内挿は、np.interp()で行うことができます。内挿したいxの配列、デー タ点の配列(xp, yp)でnp.interp(x, xp, yp)に渡せば内挿したyデータが得ら れます。 https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.interp.html .. literalinclude:: example-21.py :caption: example-26 結果をプロットした図は以下のとおりです。 .. image:: interpolate.png :width: 500px :align: center クイズ ====== Q1 -- linspace()とreshape()を使って以下の配列を作成する。 .. list-table:: :header-rows: 0 :widths: 100 :stub-columns: 0 * - .. image:: problem-1.png :scale: 80% :align: left 作成した配列のindexingにより黄色の配列を表示する。 .. list-table:: :header-rows: 0 :widths: 100 :stub-columns: 0 * - .. image:: problem-1-a.png :scale: 80% :align: left * - .. image:: problem-1-b.png :scale: 80% :align: left .. literalinclude:: problem-01.py Q2 -- 次の配列を直接入力して作成する。 .. list-table:: :header-rows: 0 :widths: 100 :stub-columns: 0 * - .. image:: problem-2.png :scale: 80% :align: left 作成した配列からTrue, False(boolean)を使って(深いコピー)次の配列を表示する。 (注意)何通りもやり方があります。 .. list-table:: :header-rows: 0 :widths: 100 :stub-columns: 0 * - .. image:: problem-2-a.png :scale: 80% :align: left * - .. image:: problem-2-b.png :scale: 80% :align: left .. literalinclude:: problem-02.py Q3 -- linspace()とreshape()を使って以下の配列を作成する。 .. list-table:: :header-rows: 0 :widths: 100 :stub-columns: 0 * - .. image:: problem-3.png :scale: 80% :align: left 作成した配列から、indexing、深いコピーによるスライスとbroadcastingを組 み合わせて黄色のデータを表示させる。1文でデータを表示させること(注意)何通りもやり方が あります。 .. list-table:: :header-rows: 0 :widths: 100 :stub-columns: 0 * - .. image:: problem-3-a.png :scale: 80% :align: left * - .. image:: problem-3-b.png :scale: 80% :align: left .. literalinclude:: problem-03.py Q4 -- xy平面に2Ang.離してCH4を10x10個しきつめたxyzファイルを生成するプログラ ムを作成する。以下のforループの中の「自分で考える」のところにプログラ ムを書き加えて完成させる。 .. code-block:: python import numpy as np CH4 = """ H 2.00000000 0.00000000 0.00000000 H 3.25515949 1.25515949 0.00000000 H 3.25515949 0.00000000 1.25515949 H 2.00000000 1.25515949 1.25515949 C 2.62757974 0.62757974 0.62757974 """ # Cを原点にする base = np.array(CH4.split(), dtype=object).reshape((5, 4)) base[:, 1:4] = base[:, 1:4].astype(float) base[:, 1:4] = base[:, 1:4] - base[base[:, 0] == "C", 1:4] data = np.empty((0, 4)) for i in range(10): for j in range(10): """自分で考える""" # xyzファイルの出力 outfile = "CH4-problem.xyz" o = open(outfile, "w") o.write("{}\n".format(len(data))) o.write("CH4 rotate and shift\n") for d in data: o.write("{} ".format(d[0])) o.write("{} ".format(d[1])) o.write("{} ".format(d[2])) o.write("{} ".format(d[3])) o.write("\n") print(outfile, "was created.") 生成するxyzファイルは下図のような構造です。 :download:`CH4-problem.xyz ` .. image:: problem-03.png :width: 500px :align: center .. literalinclude:: problem-04.py Q5 -- 以下の連立方程式をnp.solve()を使って解いて :math:`a_0, a_1, a_2, a_3` を求める。 .. math:: 3a_0 + 8.2a_1 = -2a_3\\ 2.3a_2 + 9.3a_3 = 0\\ -2.3a_0 + 1.8a_1 + 8a_2 + 7.2a_3 = 13\\ 0.8a_0 + 7.3a_1 = 5.3a_2 + 9a_3 .. literalinclude:: problem-05.py Q6 -- 300K, 600K, 1500KでのArの速度分布をマクスウェル速度分布の式に従ってプ ロットする。以下のプログラムのHOGE部分を自分で考えてプログラムを完成さ せる。プロットすると下図のようなグラフが得られます。 .. math:: &f(v) = 4\pi v^2\left(\frac{m}{2\pi R T}\right)^{3/2} \exp\left(-\frac{mv^2}{2RT}\right)\\ &m: {\rm Ar}質量\ {\rm kg/mol}\\ &R: ガス定数\quad 8.3144621\ {\rm J/K/mol}\\ &T: 温度\ {\rm K} .. code-block:: python import numpy as np import matplotlib.pyplot as plt R = 8.3144621 # J/K/mol m = 39.948*1e-3 # kg/mol fig, ax = plt.subplots() v = np.linspace(0, 2000, 1000) # m/s A = 4*np.pi*v**2 f = HOGE g = HOGE h = HOGE ax.plot(v, f, label="300K") ax.plot(v, g, label="600K") ax.plot(v, h, label="1500K") ax.legend() ax.set_xlabel("velocity m/s") ax.set_ylabel("probability") fig.savefig("problem-05.png") plt.show() .. image:: problem-05.png :width: 500px :align: center 作成したプログラムを改造して以下の計算を行う。 * 速度分布が規格化(積分すると1)されていることを確認する。300Kで得た分布の各点を長方形として、全領域の足し算で良いです。 * 300K, 600K, 1500Kでもっとも頻度の高い速度 * 300K, 600K, 1500Kの平均速度 * 300K, 600K, 1500Kで1000 m/s以上のArの割合(%) すべて数式として解析的に得ることができますが、プログラムを使ってデータから計算すること。 .. literalinclude:: problem-06.py Q7 -- 以下の図のような1辺が100 Ang.の立方体にArが10000個入った :download:`Ar.xyz <./Ar.xyz>` がある。 立方体の中心座標は、(x, y, z) = (0, 0, 0)になっている。 .. image:: problem-06-1.png :scale: 80% :align: center * x座標が-2から2の範囲に存在するArを抜き出してxyzファイルを作成する。 * 立方体の中心から40 Ang.以内に存在するArを抜き出してxyzファイルを作成する。 それぞれ以下のようなxyzファイルが得られます。 .. image:: problem-06-2.png :scale: 70% :align: center .. literalinclude:: problem-07.py Q8 -- 温度300Kの容器に存在するArガスの速度ベクトル(8000個のAr分)のデータ :download:`Ar_velocity.dat <./Ar_velocity.dat>` を読んで、np.histogram()により速度のヒストグラムを作成しプロットする。問題5と同じ分布関数になりますが、区分の幅(bins)で高さが変わります。速度範囲を20分割してヒストグラムを作成しプロットした結果を下に示します。 .. image:: problem-07.png :width: 500px :align: center .. literalinclude:: problem-08.py Q9 -- 反応次数が1(一次反応)の化学反応の場合、反応物の濃度 :math:`[A]` は時間の経過とともに初期濃度 :math:`[A]_0` から以下の式に従って減少する。 .. math:: &[A] = [A]_0\exp\left(-kt\right)\\ &k(T): 速度定数\ s^{-1}\\ &t: 時間\ s\\ 両辺の自然対数をとれば以下の一次の関係式が得られる。 .. math:: \ln([A]) = \ln([A]_0) - kt 以下の反応について、時間と :math:`{\rm [N_2O_5]}` のデータを実験から得た。 .. image:: problem-08-01.png :width: 260px :align: center 得られたデータは以下のとおり。np.polyfit()を用いて直線フィットし反応速度 :math:`k` と初期濃度 :math:`{\rm [N_2O_5]_0}` を求める。 .. code-block:: none Time/s [N2O5]/M 50 0.0717 100 0.0511 200 0.0250 300 0.0125 400 0.00685 実験データから以下のような図を作成すれば、 :math:`k` と :math:`{\rm [N_2O_5]_0}` が求まります。 .. image:: problem-08-02.png :width: 500px :align: center .. literalinclude:: problem-09.py