コンテンツへスキップ

研究などでpythonを使って計算したりグラフ描画したりしたい人のための基礎事項まとめ本」でした。

「研究などでpythonを使って計算したりグラフ描画したりしたい人のための基礎事項まとめ本」でした。

科学技術計算のためのPython入門――開発基礎,必須ライブラリ,高速化 技術評論社 中久喜 健司 (著)

期間: 2019年9月4週目~10月1週目くらい

内容(個人的まとめ)

  1. pythonプログラミングの例
  2. IPython,Spyderなど
  3. pythonの基礎事項
  4. NumPy
  5. pandas
  6. SciPy
  7. Matplotlib
  8. 高速化

特徴

科学技術計算でPythonを使用したいと考えている方に向けた本になっていました。特定の何かをしようというものではなく、例が豊富で実行しながらpythonでできることを一通り網羅的に紹介して貰える本です。

感想

(参考として)読む前:pythonのフロー制御とモジュールを使って簡単な計算を自動化するプログラムが作成できる程度の知識がありました。一方で、Open CVなどその時々に必要なモジュールがないかググり、Qiita等に書いてある範囲内で使用することしかできませんでした。

全体

読む以前から多少の知識がある状態で、基本事項を改めて確認し直したい私にはとてもちょうどいい本でした。もしかしたらpythonを書いたことのない方には難しいかもしれませんが、例が豊富なのでひとまず読むことはできるかと思います。

公式リファレンスを読むのはまだ重い...という私と同じような方におすすめです。

1. pythonプログラミングの例

プログラムを書き、テストし、高速化して利用するという流れを一通り見せてくれようとしている章でした。作業フローを教えてくれる本はあまりであったことがなかったので貴重です。

例としてロケットシミュレーターを挙げており、公開されているコードを落としてこれるので実行しながら読み進められるはずです。が!地図をプロットするmatplotlib.basemapがエラー吐きまくりだったのでコードの実行はあきらめました。

2. IPython,Spyderなど

私はJupyter Notebookで普段書いています。Jupyter Notebookで書くとグラフや画像がその場で表示されるしTab補完もできるのでで重宝しています。

IPythonとはどんなものか、Spyderとは何かといった、知らなくてもどうにかなるが知っていると便利、でもわざわざ調べないことが書いてあり有難かったです。

3. pythonの基礎事項

一応知っていた内容が多い部分でした。「インデキシングとスライシングっていう名前の上での分類があったんだー」というような、教科書改めて読むと面白いという感じの内容です。

4. NumPy

NumPyいつもimportするけどnp.piしか使わないじゃん、と思っていた俺を殴りたい。多次元配列ndarrayとはどのようなものか、何がいいのかということを、例を実行しながらイメージできました。

5. pandas

よく聞くpandasとNumPyの違いとは何か、どんな時にpandasを使えばいいのか、なんとなく分かった気になれます。主にデータ型Series,DataFrame,Panelsに関して書かれています。

ただし、Panelsに関しては、「将来的になくなるよ、使わない方がいいよ」と警告してもらったので(下記)読み飛ばしました。本の悪い部分は情報の更新がわからないことなので仕方がないですね。

Panel is deprecated and will be removed in a future version.
The recommended way to represent these types of 3-dimensional data are with a MultiIndex on a DataFrame, via the Panel.to_frame() method
Alternatively, you can use the xarray package http://xarray.pydata.org/en/stable/.

6. SciPy

SciPyの項目では、以下の項目について例を挙げていました。

  • 統計関数
  • 離散フーリエ解析
  • ボーデ線図
  • データの内挿(グラフに近似線を追加するなど)
  • デジタル信号フィルタの設計
  • 行列の分解

7. Matplotlib

使いたいときにググればいいと思い、読み飛ばしました。

8. 高速化

自分が今必要としていない部分だと思いましたので、読み飛ばしました。

次に読みたい

関連記事

CSVファイルを指定して極座標系でグラフを描き、画像に出力するまでをやります。

ファイルを指定して極座標系でグラフを描き、画像に出力するまでをやります。

目次

  • 作業環境
  • CSV
  • 極座標系
  • ファイルのインポート

作業環境

  • anaconda 1.7.2
  • python 3.7.3
  • vscodeとexcel

jupyter notebookというIPythonを使用します。

jupyter notebook

jupyter notebook はサーバを立ててwebブラウザ上で文章を表示させるものです。 テキストエディタと同じ感覚で使用でき、IPythonの機能であるコードの解説の表示や図を含む実行結果をまとめて表示できます。

図を描くときは、実行してすぐ図を確認できるjupyterが楽です。anacondaをインストールするとまとめてインストールしてくれます。

CSV

csvって何

カンマで区切られた数値です。

VS code上で見やすく

csvはカンマで区切られているだけのデータなので見にくいです。そこで私が使っているvscode拡張機能を紹介します。

Rainbow CSV

Rainbow CSV

Excel VIewer

Excel viewer

ファイルを右クリックかcommand palletでpreviewが見れます。

極座標系の表示

=>matplotlib.pyplot.polar -matplotlib

import matplotlib.pyplot as plt
import numpy as np

r = 1
theta = 45
theta_rad = np.deg2rad(theta)# 角度をradに変換

ax=plt.subplot(projection="polar") #軸を表示
ax.scatter(theta_rad, r, c="red", s =15)

theta = np.arange(0.0, 4*2*np.pi, 0.01) #角度
r = 0.5*theta #半径
plt.polar(theta,r) # 極座標グラフのプロット

plt.savefig('PolarCoordinate.png')
plt.show()

subplotについて

fig.add_subplot(行,列,場所)を表します。(1,2,1)では列が(つまり縦方向に)半分になったグラフのうち1つめ(上側)が1つ表示されるらしいです。

import matplotlib.pyplot as plt

fig = plt.figure()
fig.add_subplot(1,2,1)
plt.show()
fig.add_subplot(1,2,1)
ax1 = fig.add_subplot(2,2,1)
ax2 = fig.add_subplot(2,2,2)
ax3 = fig.add_subplot(2,2,3)

ファイルのインポート

使っているデータは私が他の目的で作ったデータを適当に弄った意味のない数値列です。

import pandas as pd

df = pd.read_csv('fact_x2.csv')
print(df)
     turn[] length[mm]   total [mm]  radius[mm]    deg[deg]    \
0         0       0.00         0.00      0.0000    0.000000     NaN   
1         1      11.45        11.45      7.4000   88.653605     NaN   
2         1       0.45        11.90      7.4000   92.137808     NaN   
3         1      11.45        23.35      7.4000  180.791412     NaN   

~割愛~

117      17      16.10       854.45      9.4496   52.276096     NaN   
118      17       1.15       855.60      9.4496   59.248893     NaN   
119      17      16.10       871.70      9.4496  156.868057     NaN   

[120 rows x 10 columns]
import pandas as pd

df = pd.read_csv('fact_x2.csv')
print(df['turn[]'])
0       0
1       1
2       1
3       1

~割愛~

117    17
118    17
119    17
Name: turn[], Length: 120, dtype: int64

​指定できるのは行(Colum)。 横に伸びてる方です。

与える角度はdegじゃなくてradですので注意。

完成

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

df = pd.read_csv('fact_x2.csv')
r = df['Radius[mm]']
deg = df['Deg[deg]']
radian = np.deg2rad(deg)

fig = plt.figure(figsize=(10, 10))#新しい図面
ax = fig.add_subplot(111, projection='polar')#軸の表示
ax.plot(radian,r, '-o')
i=2
while i < 122:
    ax.plot([radian[i],radian[i-1]],[r[i], r[i-1]], 'or-')
    i +=2
print(np.pi)
plt.savefig('PolarCoordinate.png')
plt.show()

関数化

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

"""
CSVファイルから極座標系のグラフをプロットする
入力:(CSVファイル名, 出力画像ファイル名,半径、角度)
出力:出力画像ファイル名.png
"""
def polar(in_file,out_file,param_r,param_theta):
    df = pd.read_csv(in_file) #ファイル読み込み
    r = df[param_r]
    deg = df[param_theta]
    radian = np.deg2rad(deg) #deg をradに変換
    Len=len(r)#配列の長さを取得

    fig = plt.figure(figsize=(10, 10))#新しい図面
    ax = fig.add_subplot(111, projection='polar')#軸の表示
    ax.plot(radian,r, '-o') 
    i=2
    while i < Len: #1つ飛ばしどうしを赤線で結ぶ
        ax.plot([radian[i],radian[i-1]],[r[i], r[i-1]], 'xr-')
        i +=2
    plt.savefig(out_file)
    plt.show()

polar('test.csv','Polar_test.png','Radius[mm]','Deg[deg]')

雑記

進捗がない!!!!!!!!!!!
[追記](2019/9/19)「完成」のコードを変更しました

参考

無料磁場解析ソフトFEMMチュートリアルで使い方を確認しつつその時の記録を残しておきます。

チュートリアルで使い方を確認しつつその時の記録を残しておきます。

目次

  • FEMMとは
  • 動的磁場有限要素法とは
  • インストール
  • 磁場解析をしてみる:空心コイルソレノイドの磁場解析
  • 解析結果
  • まとめ
  • 雑記
  • 参考

FEMMとは

低周波数の磁気学と静電気学について,2次元平面と軸対象問題を解くためのソフトウェアです。Windows上で動作します。

ソルバという、複数の変数を含む数式において目標とする値を得るための最適な変数の値を求めることができる機能がついています。FEMMのソルバは 動磁場有限要素法を使用しているものみたいです。

Kelvin変換,Robin境界,SDI境界,IABC境界など多数のOpen boundary 境界が実装されている、みたいですが、解析方法とかよくわかりません。

チュートリアルにはLuaという言語と統合できるみたいなことが書いてありますが、PythonでもpyFEMMpyfemmを使って動かせるみたいです。 FEMM自体はC言語で動いているみたいです。

なお、この記事では磁場に関しては深く言及しません。

 ソルバー -weblio辞書

動的磁場有限要素法 とは

Finite Element Method (有限要素法)は、微分方程式を近似的に解くための数値解析方法です。複雑な形状。性質を持つ物体を小部分に分割することで近似し、全体の挙動を近似して計算機に計算させようとする方法です。要素の選択とメッシュの作り方を近似していきます。 構造力学や流体力学などの様々な分野で使用されているみたいです。

mesh

要素の円卓とメッシュの作り方を指定してあげます。 FEMMは固体要素を扱うことが多く、メッシュを正三角形で行っています。一般的に、メッシュを細かくすればするほど計算精度が良くなりますが、その分計算時間が掛かります。また、三角形よりも四角形、の方がj計算精度がいいです。

端的に言えば、各要素で編微分方程式を解いてるみたいです。(詳細についてしゃべるのを避けた)

インストール

Finite Element Method Magnetics :Download

からダウンロードします。

download画面

磁場解析をしてみる:空心コイルソレノイドの磁場解析

以下がバージョンが違いますが、チュートリアルです。

FEMM 3.4 Magnetostatic Tutorial1

一部変更していますが、重要な部分だけやります。

今回の解析対象はこれです。

解析対象

2Dの解析ソフトなので、FEMMで一回に解析するのはこれの断面になります。これは軸対称のモデルなので描くのは左側の四角になります。

新規モデルの作成

File/NEWでMagnetic problemを選択します。静磁気問題になります。次に問題の種類を指定していきます。axisymmetric(軸対称)を選択し、Length Unitをmm、Frequencyを0Hzにします。

平面で解析したい場合はPlanarを選択します。

境界線、コイルの描画

有限要素ソルバは対象のオブジェクトを含む空間の有限領域を指定して解析するみたいです(?)。なので、空間を指定していきます。 視界の設定をします。View/Keyboadから、Bottom:-4,Left:0,Top:4,Right:4に指定します。スクリーンは指定された領域を含む 最も小さい長方形に設計しなおされます。

次にOperate on nodesを選択します。小さいブロックボックスで左端にあるこれです。ノードを(0,4),(0,-4),(0,0)にマウスから設置します。Tabを使用すると、キーボードから入力できます。

次に、Operate on segmentsを選択します。左から2番目です。(0,4),(0,-4)をクリックすると線が引かれます。画像は見にくかったので変更しました(上記参照)。 次に、Operate on arc segmentsを選択します。

左から3番目です。(0,4),(0,-4)をクリックすると軸対称な円弧が描かれます。propertiesが現れますが、そのときArc Angleを180、Max segment Degreeを25にします。FEMMでは曲線は直線を繋げたものとして表されるので、解像度のように考えればよいでしょう。こおでは2.5で充分であるとします。

もし失敗した場合はEdit/Undoか、上の切り返す矢印マーククリックして選択した後にメニューのX印を押しましょう。

同じ手順で内側の四角を描きましょう。Operate on nodesを押して、(0.5,-1)(1.5,-1)(1.5,1)(0.5,1)にノードを置く。 その後Operate on segmentsで線で結びます。

材料や特性の関連付け

Operate on Block Labelsを選択します。コイル領域とコイル領域の外側の空気部にブロックラベルを置きます。

このプログラムは、問題を形状、材料、特性によって解くので、それらを関連付けていきます。

メニューからProperties/Materials Library を選択しますこの中から、Air とCopper AVW Magnet wire/18 AGWを右のModel Materialsにドラック&ドロップで加えます。OKを押します。

コイルに回転属性を追加します。

properties/CircuitsからAdd Propertyを円卓し、new circuitに適当な名前をつけます。Seriesを選択して巻き線領域になるように指定します。これは時間調和問題に使用される設定です。 Circuit Currentに1を入力します。

円の中、四角の外のブロックラベルを右クリックで選択し、スペースか上のバーの指を指しているようなマークでpropertyを開きます。

Block typeを先ほどのAirにし、Let Triangle Choose Mesh Sizeのチェックを外して0.1にします。このとき、1辺最大0.1の正三角形で内部を満たそうとします。

四角の中はBlock type 18 AVW、Mesh size 0.1、In circuit Coil、Number of turn 1000を入れます。これで、時計回りに1000巻きの18材のコイルが指定されます。右ねじの方向が正。-1なら逆時計回りに1。

境界条件

境界はboundaryです。Prorerties/Boundaryを選択し、Add Propertyを選択し、nameをABCに、BC typeをMixesにします。ABCは(チュートリアルによると)境界の無いオープンスペースのイ ンピーダンスに近似した,-漸近的な境界条件. を作成している事を意味するらしいです。

c0 coefficient(係数) と c1 coefficientは以下で表されます。

$$ \frac{1}{\mu_r\mu_o}\frac{\partial A}{\partial n} + c_0 A + c_1 = 0$$

$$A:磁気ベクトルポテンシャル、\mu_r:境界に近い領域の比透磁率$$

$$\mu_o:真空の透磁率、n:境界に対するの標準の方向$$

ここでは、以下のようにします。

$$c_0=\frac{1}{\mu_r\mu_o R}$$ $$c_1=0$$

Rは球の外径です。今回の入力は、c0 coefficient に 1/(uo4mm)、c1 coefficientに0を入力します。uoに対する真空の透磁率の数値やmに当たる0.001を代入し、OK、OKをクリックします。 ちなみに1inch=0.0254mです。 空間の有限領域に着いて、無限境界空間にあるコイルによって作成されたフィールドをモデル化出来ました。

解析の開始

境界条件を指定するためにOperate on arc segmentsモードにします。外側の境界を意味する円弧を右クリックで選択し、スペースでpropertiesを開きます。Boundary condにABCを選択します。 バーの黄色のMesh generaterをクリックします。すると、このように塗りつぶされます。

FEA(有限要素法)による解析を実施します。turn the crankを押して解析を行ってください。 今回の解析は一瞬でしたが、プログレスバーが動いていない様子であった場合は境界条件の指定がうまくいっていないかも知れないので、条件をもう一度見直す必要があるかも知れません。

解析結果

眼鏡のアイコンをクリックすると解析結果がみれるタブが現れます。上のバーにポイント、等高線(コンター)、エリアを選択できます。

一点の値

左クリックで任意の点をクリックすることで、そのポイントの様々なフィールドプロパティがFEMM outputウィンドウに現れます。TABでより精密にポイントを選択できます。 消してしまった場合はView/outputで再度表示出来ます。上の方(0,4)真ん中右の方(3,0)を指定したときはこんな感じでした。

(0,4)
(3.0)

コイルの末端特性

コイルアイコンをクリックします。Coilプロパティはこのようになります。問題が線形で、一つの電流しかないので。Flux/Currentの結果はコイルのインピーダンスと解釈できます。Flux/Current = 0.897763 mH、Voltage/Current =R = 0.13155 Ωですね。

等高線(コンター)に沿ったフィールド値のプロット

コンターアイコンを押し、(0,4)(0,0)(0,-4)付近をクリックしてグラフアイコンをクリックします。デフォルトでは磁束密度の振幅がプロットされますが、ほかにも選択できます。 ポイントが選ばれた順番に、進行方向からみて↑右側を選択するようなプログラムになっています。

磁束密度のプロット

デフォルトでは白黒ですが、ツールバーでカラーにできます。

雑記

jMAGという有料のが研究室にあるので、使っていきたいですね。これはノートに入れられますが、jMAGは3D解析が出来るので。今後は下をやりたいですね。

  • pythonで様々なパターンを自動的に解析できないか探す
  • 他の解析ソフトを使用して今回の結果を検証

[追記2019/10/1]チュートリアルのリンクを変更しました

参考