gnuplot上で行う任意関数のフィッティング

gnuplotには、”fit”というコマンドがあります。
これを利用すると、データ列に対して、任意の自由度を含んだ関数\(f(x)\)でフィッティングを行うことが出来ます。
手法はMarquardt-Levenberg法に基づいた非線形最小二乗法によるフィッティングです。

gnuplot ver4.6の場合

データ列”data.d”を\(f(x)=ax^2+bx+c\)でフィッティングする\(a,b,c\)を求めたい場合、gnuplot上で

f(x)=ax**2+bx+c
fit f(x) "data.d" u 1:2 via a,b,c

※gnuplot ver5.0以降の場合、エラーバーで重みづけしたフィッティングが出来るようです。
http://gnuplot.sourceforge.net/demo_5.0/fit.html

説明


以下のデータ列(ファイル名”data.d”)に対してgnuplot上でフィッティングを行います。

データは、関数
\(
\displaystyle f(x)=\frac{1}{1+e^{x-15}}
\)

に乱数を与えて生成したデータです。

このデータ列に対して、
\(
\displaystyle f(x)=\frac{a}{1+e^{(x-b)/c}}
\)

でフィッティングを行い、未知の定数\(a,b,c\)を決めたいと思います。

gnuplot上で

f(x)=a/(1+exp((x-b)/c))
a=1e0
b=10e0
c=0.5e0
fit f(x) "data.d" u 1:2 via a,b,c

と打つと、\(a,b,c\)にフィッティングした後の定数が入ります。
ここで、フィッティング前に\(a,b,c\)の値を入れているのは初期値です。
初期値があまりにも違うとうまくフィッティングが失敗するので、出来るだけ近い値を入れましょう。

実際に動かすと

> fit f(x) "data.d" u 1:2 via a,b,c
...

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 0.995618         +/- 0.005679     (0.5704%)
b               = 15.3268          +/- 0.051        (0.3327%)
c               = 0.970465         +/- 0.04439      (4.574%)

...
>

という文章が出力されます。

フィッティングの結果、
\(a=0.995618,b=15.3268,c=0.970465\)と求められました。
乱数を与える前のデータは\(a=1,b=15,c=1\)ですので、元の関数に近いことが分かります。
fitコマンドを動かした後は既に\(a,b,c\)に値が代入されているので、そのまま

plot f(x)

とすれば、フィッティング結果をプロットすることが出来ます。

データと共に載せれば以下のようになります。

フィッティングする範囲を制限したければ、

fit [10:20] f(x) "data.d" u 1:2 via a,b,c

とすると、\(x\)の範囲を\([10,20]\)に制限することが出来ます。

2変数関数のフィッティング


2変数関数であろうとフィッティングは可能です。
元のデータ
(https://slpr.sakura.ne.jp/qp/supplement_data/data2.d)
は、関数
\(
g(x,y)=x\sin(xy+0.5)
\)

に適当な乱数を足して作られたデータです。グラフにすれば

となります。

このデータに対して、gnuplot上で関数
\(
g(x,y)=ax\sin(bxy+c)
\)

によるフィッティングを行います。コマンドは

g(x,y)=a*x*sin(b*x*y+c)
a=1.3e0
b=0.7e0
c=0.3e0
fit g(x,y) "data2.d" u 1:2:3 via a,b,c

です。
実行すると

>fit g(x,y) "data2.d" u 1:2:3 via a,b,c

...
Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 1.16313          +/- 0.004072     (0.3501%)
b               = 1.00073          +/- 0.001672     (0.1671%)
c               = 0.495134         +/- 0.003432     (0.6932%)
...

>

と得られますので、フィッティング結果と共にプロットすれば

となります。

初期値が答えに割と近いですが、これは失敗を防ぐためです。全ての初期値のパラメータを0で行った所失敗しました。

参考サイト


Gnuplotでの関数fitting
gnuplot demo script: fit.dem -ver5.0
gnuplot demo script: fit.dem -ver4.6
Fit (manual)
gnuplot で最小二乗フィッティングする -ゴルディアスの涙目

PV数の推移

シキノートは2014年12月に開始したので、現在(2018年8月)までに3年半経過しました。
おかげ様です。

自己満足のためにシキノートの閲覧データを公開します。

PV(ページビュー)数の推移



灰色の二つの線は線形補間と二次関数による補間です。
線形補間では負の切片を持ってしまったため、
恐らく二次関数的にPV数が増えている、というのが正しいでしょう。嬉しい限りです。

年ごとのPV数の推移は(1年は1月~12月)
2014年:64
2015年:23,993
2016年:67,647
2017年:125,982
2018年:141,141(8月17日まで)
となっています。

開始から43ヶ月
合計PV数:358,827
記事数:154
なので、1ヶ月当たり3.6記事です。

Analyticsのデータ


現在までの人気記事


タイトル 総PV数
gnuplotのカラーマップ 24,016
球体の空気抵抗と係数 17,803
Nexus5の購入とセルスタンバイ問題の解決 15,352
ルンゲクッタ法の説明と刻み幅制御 15,111
4次ルンゲ・クッタ法 10,836
レジスト系によるダメージ減少量 10,373
弾道計算(BB弾)の理論 10,280
gnuplotとアニメーション 9,260
離散フーリエ変換と高速フーリエ変換(fortran90) 8,299
三角波、のこぎり波、矩形波、その他の数式 7,768
弾道計算(BB弾)の結果 7,541
fortranコンパイラのコマンド 7,455
演算子の種類と説明 6,508
latexで数式全体を小さくする 6,109
乱数 5,849
正多角形とスピログラフの数式 5,834
fortranでのエラーメモ 5,543

二体衝突後に物体が持ち得る速度

2体が弾性衝突(エネルギー保存)を起こす時を考えます。
目的は、衝突後に2体が持つ速度ベクトルを知ることです。

2体を
物体1… 重さ\(m_1\)
物体2… 重さ\(m_2\)
衝突前の速度ベクトルを\(\mathbf{v}^{(i)}_1, \mathbf{v}^{(i)}_2\)、
衝突後の速度ベクトルを\(\mathbf{v}^{(f)}_1, \mathbf{v}^{(f)}_2\)
と書くと、衝突後に考えられる速度ベクトルは、
\(
\begin{align}
\mathbf{v}^{(f)}_1&=R(\theta)(\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{(f)}_2&=-R(\theta)(\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V} \\
\end{align}
\)

と表せられます。ここで、\(\mathbf{V}\)は重心の速度ベクトルで
\(
\begin{align}
\mathbf{V}=\frac{m_1\mathbf{v}^{(i)}_1+m_2\mathbf{v}^{(i)}_2}{m_1+m_2}
\end{align}
\)

を表し、\(R(\theta)\)はユニタリーな回転行列を表します。衝突が起こる際の角度はこの式からは求められません、あくまで、衝突によって起こり得る角度です。
ここで示しているのは、いろんな方向に衝突が起こることを考えた時の表記で、衝突によって変化することが出来る速度ベクトルを表しています。

特に二次元上での衝突を考える場合、回転行列は
\(
\begin{equation}
R(\theta)=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\end{equation}
\)

と書くことが出来ます。

二体の衝突


ある2つの物体が衝突することを考えます。
下添え字を物体を区別する番号、上添え字を衝突前\(i\) (initial), 衝突後\(f\) (final)として記述すると、
物体1… 重さ\(m_1\), 位置ベクトル\(\mathbf{r}_1\), 速度ベクトル\(\mathbf{v}_1\)
物体2… 重さ\(m_2\), 位置ベクトル\(\mathbf{r}_2\), 速度ベクトル\(\mathbf{v}_2\)
と書けます。但し、ここでは衝突前、後を表す上添え字\(i,f\)は省いています。

また、重心の位置,速度ベクトル\(\mathbf{R},\mathbf{V}\)は
\(
\begin{align}
\mathbf{R}&=\frac{m_1\mathbf{r}_1+m_2\mathbf{r}_2}{m_1+m_2} \\
\mathbf{V}&=\frac{m_1\mathbf{v}_1+m_2\mathbf{v}_2}{m_1+m_2}
\end{align}
\)

と書けます。

続いて重心から見た物体1,2の位置、速度ベクトルを\(‘\)を付けて表しますと、

物体1(重心系)… 重さ\(m_1\), 位置ベクトル\(\mathbf{r}’_1\), 速度ベクトル\(\mathbf{v}’_1\)
物体2(重心系)… 重さ\(m_2\), 位置ベクトル\(\mathbf{r}’_2\), 速度ベクトル\(\mathbf{v}’_2\)
と書くことが出来ます。
ここで、重心系における速度ベクトルは明示的に
\(
\begin{align}
\mathbf{v}’_1&=\mathbf{v}_1-\mathbf{V} \\
\mathbf{v}’_2&=\mathbf{v}_2-\mathbf{V}
\end{align}
\)

と表すことが出来ます。

重心系から見ると重心の位置、速度ベクトルはもちろんゼロとなります。
\(
\begin{align}
\mathbf{R}’&=\frac{m_1\mathbf{r}’_1+m_2\mathbf{r}’_2}{m_1+m_2}=\mathbf{0}\\
\mathbf{V}’&=\frac{m_1\mathbf{v}’_1+m_2\mathbf{v}’_2}{m_1+m_2}=\mathbf{0}
\end{align}
\)

重心系から見た系の全運動量\(\mathbf{p}’\)もゼロにならなければなりません。
このことから、
\(
\mathbf{p}’=m_1\mathbf{v}_1+m_2\mathbf{v}_2=\mathbf{0}
\)

という関係式が成り立ちます。
この式が言っているのは、重心系を考える限り、物体1,2の速度ベクトルは\(\mathbf{v}’_1,\mathbf{v}’_2\)は衝突が起こる起こらないにかかわらず、上式が満たさなければならない、ということです。
繰り返しますが、衝突が起こる起こらないにかかわらず、です。

衝突が起こり、どの角度に飛んでいくか?を与えられた情報からでは決めることは出来ません。
実験結果や、物体の形状を表す別の条件式が必要となります。

ここで出来ることは、起こりうる角度を知ることだけです。
衝突後に起こりうる角度を求めてみましょう。
重心系では系の合計の運動量\(\mathbf{p}’\)がゼロベクトルにならなければなりません。
このことは衝突前後でノルムが変わる衝突は許されないことを示しています。
これは、重心系においてユニタリーな変換が許されることを意味しています。

衝突後に重心系で取りうる角度は重心系で考えられる角度方向全てなので、
ユニタリーな回転行列\(R(\theta)\)を用いて、衝突後に重心系で速度ベクトルが
\(
\begin{align}
\mathbf{p}’_1&\to R(\theta) \mathbf{p}’_1\\
\mathbf{p}’_2&\to -R(\theta) \mathbf{p}’_2
\end{align}
\)

を満たす変化のみ許される、ということです。マイナスの符号は、\(\mathbf{p}’\)がゼロでなければならないという条件からきているものです

実験系での表現に戻せば、運動量ベクトルは
\(
\begin{align}
\mathbf{p}^{f}_1&=m_1R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+m_1\mathbf{V}\\
\mathbf{p}^{f}_2&=-m_2R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+m_2\mathbf{V}
\end{align}
\)

であり、速度ベクトルは
\(
\begin{align}
\mathbf{v}^{f}_1&=R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{f}_2&=-R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V}
\end{align}
\)

と書き表すことが出来ます。
特に二次元の場合、デカルト座標において回転行列は
\(
\begin{equation}
R(\theta)=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\end{equation}
\)

と書くことが出来るので、
\(
\begin{equation}
\left(
\begin{array}{r}
v^{(f)}_x\\
v^{(f)}_y
\end{array}
\right)
=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\left(
\begin{array}{r}
v^{(i)}_x-V_x\\
v^{(i)}_y-V_y
\end{array}
\right)
+
\left(
\begin{array}{r}
V_x\\
V_y
\end{array}
\right)
\end{equation}
\)

と表現することが出来ます。

衝突角度について


衝突角度を簡単に推測してみましょう。
衝突に際し重要なのは、物体と物体が接する点における法線ベクトル\(\mathbf{n}\)だと推測できます。
恐らく、実験系において法線ベクトルを軸として、反転させた入射ベクトル\(-\mathbf{v}^{(i)}_{1,2}\)を\(\pi\)回転させる角度になるでしょう、と予測できます。
これは私の直感なので間違っているかもしれません。

ともあれ、必要な関係式は\(\mathbf{v}^{(i)}_1,\mathbf{v}^{(i)’}_1 \)を軸とした時の、実験系の角度\(\theta_1\)と重心系の角度\(\theta’_1\)との間の関係式です。

衝突後の速度ベクトルは
\(
\begin{align}
\mathbf{v}^{(f)}_1&=R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{(f)}_2&=-R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V}
\end{align}
\)

で表されていますので、実験系での角度、すなわち\(\mathbf{v}^{(i)}_1\)と\(\mathbf{v}^{(f)}_1\)のなす角度\(\theta_1\)は
\(
\begin{align}
\theta=\text{acos}
\left(
\frac{\mathbf{v}^{(i)}_1\cdot \mathbf{v}^{(f)}_1}{|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|}
\right)
\end{align}
\)

となります。
今、弾性散乱を考えているので、必要であれば\(|\mathbf{v}^{(i)}_1|=|\mathbf{v}^{(f)}_1|\)が成り立っています。
なので、
\(
\begin{align}
\theta=\text{acos}
\left(
\frac{\mathbf{v}^{(i)}\cdot [R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}]}{|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|}
\right)
\end{align}
\)

です。
変形すれば、
\(
|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|\cos\theta_1-\mathbf{v}^{(i)}_1\mathbf{V}=\mathbf{v}^{(i)}_1R(\theta’)(\mathbf{v}_1^{i}-\mathbf{V})
\)
という関係式が導けます。
この式は\(\theta\)と\(\theta’\)を結びつける式になっていますので、実験系での角度が分かれば、重心系での角度に焼直せるため、衝突した時の方向も決められるのです(簡単には求められなさそうですが・・・)。

ゼロ点を探す(二分法、挟み撃ち法、Anderson-Björk法、Brent法、Newton法、Steffensen法)

ここでは、非線形方程式の解を数値的に求める方法である、
二分法、
挟み撃ち(False position)法、
Anderson-Björk法、
Brent法、
Newton法、
Steffensen法
等について紹介します。

  1. 問題設定
  2. まとめ
  3. 数値解法
  4. 複数のゼロ点を見つけたい場合

問題設定


区間\([a,b]\)で定義された実関数\(f(x)\)に
\(
f(x)=0
\)

を満たす\(x\)が1つある。\(x\)を求めよ。但し\(f(a),f(b)\ne 0\)。

まとめ


1変数の場合:Brent法(二分法と逆2次補間、挟み撃ち法の組み合わせ)
多変数の場合:Newton法
を使うのが良いでしょう。

数値解法


1変数の問題で、
教育上良く使われる方法は二分法
実用的な手法はBrent法
です。
またNewton法は多変数の場合にも複素数の場合でも容易に拡張できるので、その意味で実用的です。

区間\([a,b]\)にゼロ点が1つ存在する事を数学的に言えば、\(f(a)f(b)\lt 0\)ということです。

二分法


要点

・解が確実に見つかります。
・関数の振る舞いが非連続であっても、たちが悪い関数でも確実ですが、収束は遅いです。

計算方法

二分法は位置\([a,b]\)における関数の符号の情報だけを見てゼロ点を推測します。
なので、区間内にゼロ点がある事が分かっていれば、そこを境として符号が違うので確実に解を囲い込んでいくことが出来ます。
教育的に分かりやすい方法であり、根を探すにあたり失敗が無い方法ですが、
滑らかな関数であっても計算時間がかかるため、あまり実用的ではありません。

収束速度について

2分法では、1回の計算あたり範囲が\(1/2\)になるため、漸化式
\(
\displaystyle \varepsilon_n=\varepsilon_{n-1}\cdot 2^{-1}
\)

が成立するため、
\(
\varepsilon_n=\varepsilon_{0}\cdot 2^{-n}
\)

が導けます。ここで\(\varepsilon_{0}\)は初期のxの範囲(\(\varepsilon_0=|b-a|\))を表します。
収束精度と計算回数の関係を決めるためにはnについて解いて、
\(
\begin{align}
\varepsilon_n &=2^{-n} \varepsilon_{0} \\
2^n &=\frac{\varepsilon_0}{\varepsilon_n} \\
n &=\log_2\left(\frac{\varepsilon_0}{\varepsilon_n}\right) \\
\mbox{または}& \\
n &=\frac{\ln\left(\varepsilon_0/\varepsilon_n\right)}{\ln{2}}
\end{align}
\)
より、例えば初期の解の範囲(\(\varepsilon_0\))を1, 解を\(10^{-14}\)まで収束させようとすれば、大体\(n=47\)となり、47回、2分法を繰り返す必要があります。
逆に言えば、どんな状況でも47回繰り返せば相対誤差を14桁収束させることが出来るということです。

収束の判定は
要求精度\(\epsilon_{\text{tol}}\)、機械精度\(\epsilon_{\text{mac}}\)、解の存在範囲\(\delta\)、新たな解の推定値\(x\)、
とすると、
\(
\delta \lt 4\epsilon_{\text{mac}}|x|+\epsilon_{\text{tol}}
\)

もしくは
\(
f(x)=0
\)

で判定します。
この意味は、解の存在範囲を\(\pm 2\epsilon_{\text{mac}}\)の範囲まで特定するか、解の存在範囲が要求精度\(\epsilon_{\text{tol}}\)まで達するかのどちらかが満たされるかで判定します。


Fortranプログラムはこちら。

二分法でやると以下の通り20回必要になります。

$ gfortran bisection.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      5.2500000000000000
      2.8750000000000000      5.2500000000000000
      2.8750000000000000      4.0625000000000000
      2.8750000000000000      3.4687500000000000
      2.8750000000000000      3.1718750000000000
      2.8750000000000000      3.0234375000000000
      2.9492187500000000      3.0234375000000000
      2.9863281250000000      3.0234375000000000
      2.9863281250000000      3.0048828125000000
      2.9956054687500000      3.0048828125000000
      2.9956054687500000      3.0002441406250000
      2.9979248046875000      3.0002441406250000
      2.9990844726562500      3.0002441406250000
      2.9996643066406250      3.0002441406250000
      2.9999542236328125      3.0002441406250000
      2.9999542236328125      3.0000991821289062
      2.9999542236328125      3.0000267028808594
      2.9999904632568359      3.0000267028808594
      2.9999904632568359      3.0000085830688477
      2.9999995231628418      3.0000085830688477
      2.9999995231628418      3.0000040531158447
      2.9999995231628418      3.0000017881393433
      2.9999995231628418      3.0000006556510925
   3.0000000894069672    
$

6桁収束するために24回も評価回数が必要です。
計算時間を気にしないのならば良い方法ですが、やはり遅いです。

挟み撃ち法(False position法)



要点

・解が確実に見つかります。
・勢いよくゼロ点を横切っていない場合、収束が遅くなります。

計算方法

二分法との違いはゼロ点の推測に用いる情報量です。
二分法では位置\([a,b]\)における関数の符号だけを見てゼロ点を推測しますが、
挟み撃ち法では位置\([a,b]\)における関数のも見てゼロ点を推測しています。
この挟み撃ち法の推測は、点\((a,f(a)),(b,f(b))\)を結ぶ直線が通るゼロ点して値を推測するのです。

一つ、挟み撃ち法の説明を見ていた時に、解の囲い込み方法が2通りあるようです。
オリジナルのものは片方の点を絶対に変えない方法のようですが、二分法と同じように決めている物もあります。

ここでは、二分法と同じように決めていく後者の方法を載せておきます。
二分法の新しい点を決める所だけ変えればプログラムは完了です。

要求精度を6桁にすると、実行結果は

$ gfortran false.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      4.8581311942172727
      2.6310989812608572      4.8581311942172727
      2.6310989812608572      3.1997221817009311
      2.9967206978247356      3.1997221817009311
      2.9967206978247356      3.0000354381074144
      2.9999999998952847      3.0000354381074144
   3.0000000000000000    
$

となり、初めの1回は無視すると6回で収束に至っていることが分かります。早いですね。

Anderson-Björk法


アンダーソン・ビョーク法と呼ばれるアルゴリズムです。
詳しくは
求根アルゴリズム 挟み撃ち法 / イリノイ法 / アンダーソン・ビョーク法 [2/3] -サブロウ丸
のページが詳しいのでそちらを参照してください。

この方法は割線法が遅くなってしまう時の対処です。

program main
  implicit none
  double precision::a,b,tol,sol
  double precision,external::f
 
  a=10d0
  b=0.5d0

  tol=1d-6
  call AndersonBjork(a,b,f,tol,sol)
  write(6,*)sol

  stop
end program main

function f(x)
  implicit none
  double precision,intent(in)::x
  double precision::f

  f=2d0*(atan(x-3d0)+0.5e0*sin(x-3d0))
 
  return
end function f

subroutine AndersonBjork(ax,bx,f,tol,ans)
  implicit none
  double precision,intent(in)::ax,bx,tol
  double precision,intent(out)::ans
  double precision,external::f
 
  integer::i
  integer,parameter::itmax=100
  double precision::a,b,c,fa,fb,fc
  double precision::tol1,xm,m
  double precision,parameter::eps=epsilon(1d0)

  a=ax
  b=bx
  fa=f(a)
  fb=f(b)
 
  if ((fa .eq.0.0d0 .or. fb .eq. 0.0d0).or. &
       (fa * (fb/dabs(fb)) .le. 0.0d0))then
     do i=1,itmax
        !if(a.gt.b)write(6,'(2f24.16)')b,a
        !if(b.ge.a)write(6,'(2f24.16)')a,b

        c = a-fa*(a-b)/(fa-fb)

        fc=f(c)
        if(fa*fc.gt.0d0)then
           if(fc/fa.lt.1d0)then
              m=1d0-fc/fa
           else
              m=0.5d0
           endif
           fb=m*fb
           
           a=c
           fa=fc
        else
           if(fc/fb.lt.1d0)then
              m=1d0-fc/fb
           else
              m=0.5d0
           endif
           fa=m*fa
           
           b=c
           fb=fc
        endif

        tol1=2.0d0*eps*dabs(c)+0.5d0*tol
        xm = 0.5d0*(a-b)
        if ((dabs(xm).le.tol1).or.(fc.eq.0.0d0))exit
     enddo
     ans=c
  else
     write(6,*)'f(ax) and f(bx) do not have different signs,',&
          ' zeroin is aborting'
  end if

  return
end subroutine AndersonBjork

Brent法



要点

・解が確実に見つかります。
・基本的には逆二次補間(inverse quadratic interpolation)で根を探索します。
・逆二次補間が原理的に出来ない場合(二つ以上の関数値が同じ場合など)では挟み撃ち法を行います。
・”探索が失敗”した時は二分法に切り替えて根を探します。
・1変数の場合、実用的な方法として推奨されています。

計算方法

逆二次補間は3点を結ぶ逆二次関数(上下に凸ではなく、左右に凸をもつ関数)としてラグランジュ補間をして、その根をゼロ点の位置として用いる方法です。
しかし、逆二次補間を行う際に必要な3点の内、2点以上が同じ関数値を持つと解が発散する、という問題があるため、その時は挟み撃ち法に切り替えて計算を行います。

二分法は、逆二次補間や挟み撃ち法を用いた為に解の収束が遅くなる、と判定されたときに使われます。

実際に動かしてみると逆二次補間よりも挟み撃ち法が使われる頻度が多いようです。

詳細は(Van Wijngaarden-Dekker-Brent法, William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery著, 丹慶勝市ら訳(1993) 『Numerical Recipes in C [日本語版]C言語による数値計算のレシピ』 技術評論社 261pp.)に書いてありますので、こちらを参照してください。

プログラム自体はパブリックドメインのNetlibにあるもの(zeroin.f)を基本としましょう。ただし、Fortran77で掛かれているので並列計算を考えた時に厄介です。文番号を消したものを置いておきます。

要求精度を6桁にすると、実行結果は

$ gfortran main.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      4.8581311942172727
      2.6310989812608572      4.8581311942172727
      2.6310989812608572      3.5172292241144740
      2.6310989812608572      3.0069570368988154
      2.9997489497326733      3.0069570368988154
      2.9997489497326733      3.0000000032534229
      2.9999995032534215      3.0000000032534229
   3.0000000032534229    
$

となり、初めの1回は無視すると7回で収束に至っていることが分かります。挟み撃ち法とだいたい同じです。

Netlibがpublicdomainであることの記述は以下の本にも見ることが出来ます。
https://books.google.de/books?id=2GPNBQAAQBAJ&pg=PA109&lpg=PA109&dq=netlib+public+domain&source=bl&ots=d6Uw8D5JMv&sig=Dw3wXTAQThFKLVrL7pm_qnCaq9s&hl=de&sa=X&ved=0ahUKEwjO05Xxt_XKAhXCwxQKHUBlC5sQ6AEIXjAI#v=onepage&q=netlib%20public%20domain&f=false

ちなみにですが、wikipediaにあるBrent法の疑似コードを実装したのですが、うまく働きませんでした。
Brent’s method -wikipedia en
ブレント法 -wikipedia ja

Newton法



要点

・解近傍のただ1点から収束させるので必要な情報が少なく済みます。
・解近傍のただ1点上における関数の微分を用いて収束させ、一回の計算で一致桁数が倍になります。
・本当に解の近くから始めないと、失敗します。
・多変数への拡張が容易ですが、関数の微分を数値微分で代用すると効率が落ちます。

計算方法

ニュートン法について詳しくは
ニュートン法(1、2次元、多次元)
に書いたので、こちらを参照してください。


プログラムは以下のもので、

実際に実行すると

$ gfortran newton.f90
$ ./a.out
      4.0000000000000000
      2.4339000410871483
      3.0980975267667068
      2.9994762813770324
      3.0000000000829825
   3.0000000000000000
$

という結果を得ます。

番外編:Steffensen法(Steffensen’s method)


Steffensen法は自己無撞着方程式を解くために使われる方法で、Newton法に似ています。
特に良く使われるのは、方程式がたまたま
\(
x=g(x)
\)

と書けている特定の方程式を満たす\(x\)を探す場合に使われます。

ですが、
\(
f(x)=x-g(x)
\)

を考えれば
\(
f(x)=0
\)

を考える問題と同じ、ということになります。
ここでは、方程式
\(
x=x^2/2
\)

を満たす\(x\)を探す問題を考えます。

この問題は
\(
x-x^2/2=0
\)

と同じです。

計算手法は
htt://park.itc.u-tokyo.ac.jp/kato-yusuke-lab/nagai/note_070420_self.pdf
でまとめられている通り、

\(
\begin{align}
a&=g(x_n), \\
b&=g(a),\\
x_{n+1}&=x_n-\frac{(a-x_n)^2}{b-2a+x_n}
\end{align}
\)
の漸化式で求められます。

プログラムではこちら。

複数のゼロ点を見つけたい場合



複数のゼロ点を見つけたい場合、一つの方法として、

解のある区間の特定→収束

を繰り返すことで得られる、と考えられるでしょう。

具体的には計算区間を荒く区切り、符号が変わる範囲を特定、その符号が変わる範囲に対してのみゼロ点を探すのです。

上の図の縦線は等間隔に荒く探した時の、関数が正になる時に青、負の時に赤色に表示しています。
収束した結果が緑の点です。

fortranプログラムはこんな感じになります。

導関数を用いたゼロ点の探索


関数のゼロ点の間隔が非常に短い時に便利です。
等間隔に区切ってから探すという考えは同じです。
まず、導関数のゼロ点を探します。
ゼロ点が存在するならば関数のゼロ点は導関数のゼロ点に挟まれたところにあるので、その情報を元に探します。
導関数を解析的に与えていますが、もしも高精度にゼロ点を決めたい場合、二重数を用いた導関数の取得が良いかもしれません。


プログラムは例えば以下の通りです。

jpg,png画像からeps画像への変換

jpg, png画像などのラスタライズ画像をベクタ形式のeps画像に変換する方法です。

まとめ


一番良い方法は、LinuxでImagemagickのconvertコマンドを使って、

convert aa.jpg eps2:aa.eps

とする方法です。

1. Imagemagickを使う


Linuxが使える場合です。
Imagemagickのconvertコマンドで

convert aa.jpg eps2:aa.eps

とします。
eps2が無い

convert aa.jpg aa.eps

でも変換できますが、画像容量が大きくなってしまいます。

例えば、カラーマップの画像 668kBの”aa.jpg”

オリジナル画像

convert aa.jpg eps2:aa.eps
convert aa.jpg bb.eps

の二通りで変換した場合

668K aa.jpg   : 元画像
 17M bb.eps   : 通常のconvert
612K aa.eps   : eps2を付けたconvert

と大きな違いが出ます。

線のpng画像”rr.png”(元サイズ29.5kB)

であっても

convert rr.png eps2:rr.eps
convert rr.png rrn.eps

の二通りで変換した場合

 32K rr.png   : 元画像
5.5M rrn.eps  : 通常のconvert
104K rr.eps   : eps2を付けたconvert

となります。

eps画像をeps2によって再変換することもできますが、画質が許容できないくらい落ちてしまうのでやめたほうが良いと思います。
再変換に関して、恐らく何らかのオプションを付ければ大丈夫と思いますが、調べてはいません。

2. Gimpを使う


Linux, windowsが使える場合です。

Gimpで画像を開いて

ファイル→名前を付けてエクスポート

に進みます。
その後、ファイル名、拡張子を.epsにすれば、eps形式が得られます。

この方法では、
32Kの上の元画像が、4.5Mのeps画像になりました。

参考サイト

ImageMagickでeps形式に画像を変換するhttp://fatsumiyoshi.hatenablog.com/entry/2012/10/01/174354