Gnuplotでカーブフィッティング

実験で得られた数値スペクトルデータがあるとする。
これには、2つのガウシアン型のピークを持つ成分が含まれている。
ガウシアン型は $ y = A_1 \times \exp{[A_2 (x-A_3)^2]} $ 型。
この実験結果から成分を抽出する方法を解説します。

Gnuplotでのフィッティングコマンドは

# 関数定義
f(x) = a*x + b  # 任意の関数・変数
    a = 1.0     # 初期値
    b = -5.0    # 

# フィッティング
fit [xmin:xmax] f(x) "解析データ" u 1:2 via a, b

Step. 1 とりあえずプロットする

fort.100 というのが実験データ( 1列目 : x, 3列目: y )とする。

plot "fort.100" u 1:3 

見ると、バックグラウンドと $x = 8, 13$あたりにピークがありそう。
振幅はどちらも2-4くらいに見える(およその目安把握が重要)。

Step. 2 バックグラウンドを計算して、結果から差し引く

# y = ax + b 型のバックグラウンド処理

B(x) = B1*x + B2
    B1 = 0.3
    B2 = 1.0

fit [17:19] B(x) "fort.100" u 1:3 via B1, B2

plot \
"fort.100" u 1:3 w l lc "blue" t "Original",\
B(x) w l lc "black" dt (5,3) t "Back ground",\
"fort.100" u 1:(($3)-(B1*($1)+B2)) w l lc "red" lw 2 t "subtract back ground"

バックグラウンドは直線に乗ってそうな $x = 17-19$ でフィッティング。
実行結果は以下の通り。
問題なくバックグラウンド処理はできていそう。

Step. 3 ガウシアンフィッティング

# 2つのガウシアン波形を含む h(x) でフィッティング

h(x) = f1*exp( f2*(x-f3)**2 ) + g1*exp( g2*(x-g3)**2 )
    f1 = 1.0        # 振幅 : 目視レベルで入力
    f2 = -1.0       # 横幅 : ガウシアンなので負の値入力
    f3 = 8.0        # 位置 : 目視レベルで入力
    g1 = 1.0        
    g2 = -1.0
    g3 = 13.0

fit [4:15] h(x) "fort.100" u 1:(($3)-(B1*($1)+B2)) via f1, f2, f3, g1, g2, g3

plot \
"fort.100" u 1:(($3)-(B1*($1)+B2)) w l lc "black" lw 8,\
h(x) w l lw 2 lc "yellow" dt (10,2) t "fitted curve",\

結果は次の通り

B1              = 0.0999337
B2              = 2.00122

f1              = 2.99945
f2              = -0.250139
f3              = 8.00014
g1              = 0.749658
g2              = -0.500211
g3              = 13

実際は $y = 3 e^{-0.25(x-8)^2} + 0.75 e^{-0.5(x-13)^2} + 0.1x + 2$ を計算していたので、精度よくフィッティングができたといえる。

注意点

初期値は適切な値を入れること

記事一覧

数値計算の備忘録

数値計算で忘れそうなことを不定期でメモします。

PCからのほうが見やすいですが、スマホの場合は横画面にすると見やすいです。

<記事一覧>

1. オイラー法
2. ルンゲクッタ法
3. 逆行列の計算
4. LAPACKのDGESVで連立方程式を解く
5. FORTRANで一様乱数生成
6. FORTRANの関数のメモ
7. FORTRANの行列・配列関連のメモ
8. Gnuplotでカーブフィッティング

<プラズマ関連>

1. マクスウェル分布
2. 単位 eV で表したマクスウェル分布
3. 単位 eV のEEDFを m/s に変換する

<原子分子過程>

1. 導入
2. 断面積の取得
3. 反応速度係数の計算

<その他・雑記事>

1. LatexのBeamerでパワポ
2. 原子質量単位のメモ
3. ビットとバイト
4. ポアソン方程式

---

プログラムは FORTRAN で、Ubuntu OS のパソコンでやっています。

このサイトは Google Blogger で、Markdown で書いたものをHTML変換して作っています。

何か連絡がありましたらよろしくお願いします。

Beamerでフォントを埋め込む

Ubuntu上のBeamerでプレゼン資料(PDF)を作って、Windows機から発表したい。

ただし、WindowsにPDFを移して見ると、フォントが太くなってしまう。

こういうときの対処方法は、Ubuntu側でPDFにフォントの埋め込み設定をする。

• Beamerで作ったファイル : resume.pdf
• windowsに持っていくファイル : presen.pdf

Ubuntu上で以下を実行;

gs -q -dNOPAUSE -dBATCH -dPDFSETTINGS=/prepress -sDEVICE=pdfwrite -sOutputFile=presen.pdf resume.pdf

ファイルサイズが大きくなっていれば多分大丈夫。

問題 2-1

(問題) マクスウェル方程式から連続の式を導け

---

(解答)

M. EQ. は

$$ \tag{1} \nabla \times \boldsymbol{E} = - \mu_0 \frac{\partial \boldsymbol{H}}{\partial t} $$

$$ \tag{2} \nabla \times \boldsymbol{H} = \epsilon _0 \frac{\partial \boldsymbol{E}}{\partial t} + \boldsymbol{J} $$

$$ \tag{3} \epsilon_0 \nabla \cdot \boldsymbol{E} = \rho $$

$$ \tag{4} \mu_0 \nabla \cdot \boldsymbol{H} = 0 $$

式(2)より、

$$ \tag{5} \epsilon _0 \frac{\partial \boldsymbol{E}}{\partial t} = \nabla \times \boldsymbol{H} - \boldsymbol{J} $$

式(3)を時間 $t$ で微分する

$$ \frac{\partial \rho}{\partial t} = \frac{\partial}{\partial t} (\epsilon_0 \nabla \cdot \boldsymbol{E}) = \nabla \cdot \left(\epsilon _0 \frac{\partial \boldsymbol{E}}{\partial t}\right) $$

(5)を上の式に代入する

$$ \frac{\partial \rho}{\partial t} = \nabla \cdot (\nabla \times \boldsymbol{H} - \boldsymbol{J}) $$

このうち、前半の項 $\nabla \cdot \nabla \times \boldsymbol{H}$ (回転の発散)は 0 となるので、

$$ \frac{\partial \rho}{\partial t} = - \nabla \cdot \boldsymbol{{J}} $$

変形すると連続の式になる。

ポアソン方程式

ポアソン方程式の例題

(例題) 長さ $20 cm$ の両端がアースされた平行平板を考える。 平板間に均一な密度 $10^{17} m^{-3}$のイオンが分布する場合、装置中央の電位は何$V$か？ 電子はないものとする。

(解答) ポアソンの方程式 $\nabla^2 \phi = -\rho/\epsilon_0$ を用いる。

境界条件は両端で電位 $\phi=0$、装置中央で$d\phi/dx = 0$ である。

(1) イオン密度を電荷密度に変換 : $\rho = e n_i$

(2) 立式 : $$ \frac{d^2\phi}{dx^2} = - \frac{en_i}{\epsilon_0} $$

(3) 1回積分 : $$ \frac{d\phi}{dx} = - \frac{en_i}{\epsilon_0}x + A $$

境界条件より、$A=0$

(4) 2回積分 : $$ \phi = - \frac{en_i}{2\epsilon_0}x^2 + B $$

境界条件より、$$ B = \frac{en_i}{2\epsilon_0} \left(\frac{L}{2}\right)^2 $$

$$ \therefore \phi = \frac{en_i}{2\epsilon_0}\left( \left( \frac{L}{2} \right)^2 -x^2 \right) $$

装置中央の電位は

$$ \phi(x=0) = \frac{en_i}{2\epsilon_0} \left( \frac{L}{2} \right)^2 = 9 \times 10^{16} (V) $$

となる。

ビットとバイト

 

・1 byte = 8 bit

FORTRANでは

・REAL(4) は単精度 = 4byte = 32 bit で小数 6-7桁まで計算可能

・REAL(8)は倍精度 = 8 byte = 64 bit で小数 15-17桁まで計算可能

FORTRANの行列・配列関連のメモ

FORTRANの行列関連の計算

行列のチュートリアル

行列の初期化

! DO文は不要で「:」を使って一括でできる
A(:) = 0.d0     ! VECTOR
A(:,:) = 0.d0   ! MATRIX

行列のイコール

INTEGER :: A(2,2), B(2,2), i

A(1,1) = 11 ; A(1,2) = 12
A(2,1) = 21 ; A(2,2) = 22

B=A     ! イコールはこれでよし

DO i = 1,2
    PRINT*, B(i,:)
END DO

行列の積

標準ライブラリの MATMUL を使う

PROGRAM TEST

REAL(8) :: A(2,2), B(2,2), C(2,2)

A(1,1) = 4.d0  ;  B(1,1) = 3.d0
A(1,2) = 2.d0  ;  B(1,2) = 4.d0
A(2,1) = 1.d0  ;  B(2,1) = 0.d0
A(2,2) = 3.d0  ;  B(2,2) = 5.d0

C = MATMUL(A,B)  ! 積 C = AB の計算

DO i = 1, 2
    PRINT*, C(i,:)  ! 表示
END DO

END PROGRAM TEST

積が計算できれば行列xベクトルも計算できる

PROGRAM TEST2

REAL(8) :: A(2,2), B(2), C(2)

A(1,1) = 4.d0
A(1,2) = 1.d0
A(2,1) = 2.d0
A(2,2) = 9.d0

B(1) =  8.d0
B(2) = -2.d0

C = MATMUL(A,B)

PRINT*, C(1)
PRINT*, C(2)

END PROGRAM TEST2

内積 dot_product(A,B)、転置 transpose(A) も使用可。 逆行列等は他の記事を参考に。

配列関係

配列の和

配列 A(5) のとき

• SUM(A) : 総和
• PROCUCT(A) : 総乗
• MINVAL(A) : 最小値
• MAXVAL(A) : 最大値
• SIZE(A) : 要素数 (A(2,4)だと2x4=8が出力)