sikino のすべての投稿

-1のn乗の計算速度

fortran90で、
\(a=(-1)^n, ~(n=0,1,\cdots, M)\)を計算する早いアルゴリズムは

a=1
if(mod(n,2).eq.1)a=-1

です。


物理、数学をやっていると至る所で\((-1)^n\)を見受けます。
これを数値計算する際に早い計算方法は何でしょうか?

試してみるのは以下のものです。


\(
(-1)^n
\)


\(
(-1)^n
\)
※①を倍精度型で計算


\(
e^{i\pi n}
\)
※\(i\)は虚数単位、\(\pi\)は円周率


\(
\begin{eqnarray}
\left\{
\begin{aligned}
& 1 ,~~(n=0,\mbox{偶数})\\
& -1 ,~~(n=\mbox{奇数})\\
\end{aligned}
\right.
\end{eqnarray}
\)


\(
a=-a
\)で逐次計算

使ったプログラムは以下のものです。

2次元時間依存しないシュレーディンガー方程式の数値解法

2次元の時間依存しないシュレーディンガー方程式を変分原理に基づいて解きます。

対象とする問題は以下の時間依存しないシュレーディンガー方程式です。

ここで、\(H_0\)は数値計算で用いる基底関数のハミルトニアン、\(V(\mathbf{r})\)はその他のポテンシャルです。

この系の固有関数\(\phi(\mathbf{r})\)を

として\(H_0\)の固有関数で展開します。ここで、\(\varphi_{\mathbf{r}}\)は

の固有値問題の解である、固有値\(E’_n\)に属する固有関数です。また、\(\varphi_{\mathbf{r}}\)は

として規格直交化されています。

式(2)を式(1)に代入した後、左から\(\varphi^*_{\mathbf{r}}\)を掛けて全空間で積分すれば

が得られます。ここで、\(V_{m,n}\)を

と置きました。式(7)を別の表現をすれば

と書くことが出来ます。なので、左辺のハミルトニアンの行列形式を対角化すれば、求めたい系の固有値、固有ベクトルが得られるわけです。固有ベクトルが得られたら式(2)に従って波動関数を構成すれば良いのです。

これが一般的な変分原理による解法です。
以降は2次元の問題に限り考えていきます。

2次元の場合


2次元の場合かつ\(H_0\)はxだけ含む部分とyだけ含む部分とで分離することが出来る場合を考えます。

この場合、固有関数は変数分離することが出来るので、2次元の基底関数を1次元系の基底関数の直積をとして考える方法を取ることが出来ます。すなわち、x,y方向の量子数\(n_x,n_y\)を用いて量子数\(n\)と以下のように関係して構成される、と考えるのです。

ここで、\(n\)と\(n_x,n_y\)には1対1の関係があれば順番は関係ないことに注意してください。それぞれの基底関数は

を満たします。もちろん、1次元系の基底関数がそれぞれ正規直交化されていれば二次元の場合でもその関係は保たれ、

が成立します。この基底関数で式(7)を書き換えてやれば

となります。

基底関数の対応


この問題はある長径、短径で指定される楕円内に存在する格子点の個数を数える問題と同じです。

nの順番を記録しておきさえすればどんな順番でも構いません。適当に求めればよいです。
ただし、ここではシュレーディンガー方程式を解いた結果、エネルギーの低い状態から順に欲しいので、\((n_x, n_y)\)の組み合わせでエネルギーが低い順に決定していきます。

今、エネルギーは\(x,y\)方向ともに単調に増加します。
\((n_x,n_y)\)の考えられる組み合わせで、
一番小さいエネルギーは\((1,1)\)なので、まず\(n=1\)は\((1,1)\)に対応させます。
次に小さいエネルギーで考えられるのは\((1,2),(2,1)\)のどちらかです。
ここで重要なのは現在より前に出てきた\((n_x,n_y)\)の組に1を足した組み合わせしか候補にあがりません。
なので、探索するのはこの条件を満たすものだけで充分です。

これをプログラムしたものは、下のtar.gz内に”qnumber”という名前でサブルーチンに入っています。

プログラム


こちらです。
http://slpr.sakura.ne.jp/qp/supplement_data/highprecision_tise2d.tar.gz

解凍すると、中にintegrate.f90 main.f90というファイルがはいっています。
Lapackを使うのでそれにリンクしてコンパイルしてください。intel®fortranならば、

ifort -mkl integrate.f90 main.f90

で良いと思います。
実行すると、

$ ./a.out
 Solving TISE...
     10 /    100
     20 /    100
     30 /    100
     40 /    100
     50 /    100
     60 /    100
     70 /    100
     80 /    100
     90 /    100
    100 /    100
     1.807[CPU sec]
                    1   1.0000086267300243    
                    2   2.0001412369101184    
                    3   2.0001696393476047    
                    4   3.0004659233630600    
                    5   3.0006690578629835    
                    6   3.0013176072105336    
                    7   4.0029644525394392    
                    8   4.0032488578583107    
                    9   4.0091516049507794
$

という結果が得られるでしょう。ここで、** / 100という数字はハミルトニアン行列要素の計算の経過を表しています。100は行列が\(100\times 100\)の行列であるという意味です。
1 1.0000086267300243
は、対角化後の固有エネルギーの一番低い状態から順に出力しています。
デフォルトでは二次元の調和振動子ですので、
\(
E_n=(n_x+1/2)+(n_y+1/2)
\)

の順になります。一番低いエネルギーから順に\(1,2,2,3,3,3,4,4…\)が解析解となります。

具体的な計算例


三角形型のポテンシャル

\(
\begin{equation}
V(x,y)=\left\{
\begin{array}{c}
0,~~~(y\lt x \cap x\lt a \cap y\gt -a) \\
10,~~~(otherwise)
\end{array}
\right.
\end{equation}
\)
ここで、\(a=3\).

2つの四角形のポテンシャル

\(
\begin{equation}
V(x,y)=\left\{
\begin{array}{c}
0,~~~(y\lt (x+b) \cap y \gt -(x+b) \cap y\lt -(x+b)+a \cap y \gt (x+b)-a) \\
0,~~~(y\gt (x-b) \cap y \lt -(x-b) \cap y\gt -(x-b)+a \cap y \lt (x-b)+a) \\
10,~~~(otherwise)
\end{array}
\right.
\end{equation}
\)
ここで、\(a=3\)

1次元時間依存シュレーディンガー方程式を解かねばならぬ

1次元時間依存シュレーディンガー方程式を確実に解くことを目的とします。

本稿の方針は
計算時間が掛かってもいいので、高精度に解く
ことを主眼においています。対象とする問題は、
ポテンシャルに空間的、時間的に不連続な点があっても高精度に解く
事を目的とします。

ポテンシャルが滑らかで、変なことが起こらないのならば、空間や時間を離散化して解く
クランク=ニコルソン法による時間発展が簡単で優秀な方法です。

本稿では導出する過程では\(N\)次元の問題として扱っていますが、プログラムは1次元の問題のみです。

方針


空間に関する積分を適応型数値積分、時間に関する積分を適応刻み幅ルンゲ=クッタ法で行います。
要は空間も時間も離散化しないで解く、ということです。

時間依存しないシュレーディンガー方程式は適応型数値積分であるQUADPACKを使うことで求められます。
詳しくはこちらへ。

解法


原子単位系における時間依存シュレーディンガー方程式は

として書けます。
波動関数を

のように展開します。ここで、基底関数\(\varphi_n(\mathbf{r})\)はハミルトニアン\(H_0\)の固有関数で、以下の固有値問題

の固有値\(E_n\)に属する固有関数です。
この固有関数は以下の通り規格化されているとします。

関数で展開するのはハミルトニアンに含まれる空間に対する二階微分を解析的に行い、消す為です。
式(1)に代入すると、

となり、左から\(\varphi^*_m(\mathbf{r})\)を掛けて全空間で積分すると

という式が得られます。ここで、\(U_{m,n}(t)\)を

と書きました。

この\(U_{m,n}(t)\)は

の性質を満たします。すなわち、ポテンシャル\(V(\mathbf{r},t)\)が実数値関数であれば、\(U_{m,n}(t)\)はエルミート行列となります。行列形式で書くとすれば、

と表せられます(補遺1)。

ここでは行列形式で書いても見やすくなるだけで、今回は使わないので、式(6)を高精度に解くことに集中します。

高精度に解くために

この問題を解くために高精度に求めなければならない重要な部分は

1. 空間に対する積分
\(\displaystyle
\int \varphi^*_m(\mathbf{r})f(\mathbf{r})\varphi_n(\mathbf{r})d\mathbf{r}
\)


2. 時間に対する積分(連立微分方程式)
\(\displaystyle
\frac{d c_m(t;c_1,~\cdots,c_{N})}{dt}=F_m(t;c_1,~\cdots,c_{N})
\)

という箇所を、不連続箇所があっても高精度に求める手法を使えばよい訳です。

空間に関してはQUADPACKによる適応型の数値積分, 時間に関しては刻み幅制御のルンゲ=クッタ法を用いることで対処します。
各々の詳細は以下のページをご覧ください。
最速・高精度の数値積分
ルンゲクッタ法の説明と刻み幅制御

ここまで確実に解けると仮定すると、精度の不安は基底関数の数だけに依存します。
基底関数の数を積めば積むほど精度が高い計算が出来ますが、計算時間が莫大になります。
ここは計算機の性能と折り合いをつけるしかありません。

ここから考える数値計算上の基底関数\(\varphi_n(x)\)は計算区間の端でゼロになるsin基底関数で、区間[x_a,x_b]で定義され、
\(\displaystyle
\varphi_n(x)=\sqrt{\frac{2}{x_b-x_a}}\sin\left(n\pi\frac{x-x_a}{x_b-x_a}\right)
\)

です。これは、固有値問題
\(\displaystyle
-\frac{1}{2}\frac{d^2}{dx^2}\varphi_n(x)=E_n\varphi_n(x),~~E_n = \frac{1}{2}\left(\frac{n\pi}{x_b-x_a}\right)^2
\)

境界条件
\(\varphi_n(x_a)=\varphi_n(x_b)=0\)
の元での解となっています。

さて、定式化は終わりましたので、あとはプログラムして解くだけです。

莫大な計算時間


問題は計算時間が膨大にかかりすぎる事です。

計算時間を評価してみます。

・次元数を\(D\)
・1次元辺りの基底関数の数は\(N\)
・合計\(N^D\)元1次連立微分方程式を解く
・1つの連立微分方程式の右辺を計算するのに\(N^D\)個の積分
・適応刻み幅ルンゲクッタ法を\(N_{t}\)回繰り返す
・適応刻み幅ルンゲクッタ法は時間1ステップ辺り\(5\)回評価
・積分は空間の分点数\(N_{\mathbf{r}}^D\)回に比例
・ポテンシャルが対称行列であることを利用すれば\(1/2\)倍
・CPU個数を\(N_{\mathrm{cpu}}\)とし、完璧に並列計算が行われれば\(1/N_{\mathrm{cpu}}\)倍
・関数を1回呼び出すのに必要な時間を\(t_f\)

であるので、全計算時間(関数が評価される回数に比例する、と考える)は
\(\frac{5}{2N_{\mathrm{cpu}}}\times N^D_{\mathbf{r}}N_{t}N^{2D}t_f\)と求められます。

現実的な量で計算時間を見積もってみます。
10[原子単位]秒の時間発展を区間20[原子単位]メートルで計算することを考えます。

\(D=1,~ N_{\mathbf{r}}=1000,~N_{t}=2000,~N=50,~N_{\mathrm{cpu}}=1\)とすると、\(1.25\times 10^{10}\)回関数を呼び出さなければなりません。

手持ちのノートパソコン(1.9GHz、1コア使用)で、\(2.5\times 10^8\)回の関数\(f(x)=\sin(x)\)の呼び出しに掛かるCPU時間を測ると約7[cpu秒]かかりました。
ということは、関数を1回呼び出すのにかかる時間\(t_f\)は約\(3\times 10^{-8}\)秒です。

よって、\(1.25\times 10^{10}\)回の関数の呼び出しだけで約400秒、すなわち8分かかります。

現実には他の処理も入るので、8分は絶対にかかる、という事です。
仮にCPUが10個あれば40秒、CPU速度が倍程度の物を使えば20秒位になるでしょうか。
まぁまぁ受け入れられる時間になります。

2次元を計算しようと思ったら基底関数が二乗、かつ空間の分点も二乗になるので、CPUが10個,4GHz近くのCPUを使えても\(20\times 1000\times 50=1000000\)秒程度、すなわち11日です。

11日の計算は現実的ではありません。
二乗で効いてくる基底関数の数を減らすか、空間の積分の精度を犠牲にする、CPUを増やす等の対策をしなければなりません。

また、ここで示した時間は最低かかる時間ですので、実際に動かすのであればこれ以上は確実にかかるでしょう。関数の中の関数を呼び出す等の操作があれば倍々に増えていくのです。

計算量を減らす工夫


ポテンシャルが変わらない場合、\(U_{m,n}(t)\)の積分は、前の時刻と値が変わらないのでいちいち計算するのは無駄です。

数値的に”ポテンシャルが変わっていない”という情報を取り出すことが出来れば、この無駄を排除することが出来ます。

”ポテンシャル関数が変わっていない”というのはどうすれば評価できるのでしょうか?
私が考えたのは、Quadpackの自動積分を行う時に評価された引数の点における関数の値が全て一致すれば関数が変わっていない、と判断するという評価方法です。

実際にこれはうまく行きまして、ポテンシャルが一度だけ切り替わる、という問題に対して計算時間が
1085[CPU sec] → 25[CPU sec]
にまで減少しました。約1/50の計算時間になりました。

先ほどの11日で~と言っていたのはこの工夫をし無ければ、ですので、ポテンシャルの変化が殆どない問題に限れば1時間程度で計算が終わる、ということです。

もちろん、ポテンシャルが変化し続ける場合も試しましたが、変化はほぼありません。
83[CPU sec] → 84[CPU sec]
程度になりました。
この場合は2次元の例で出した通り、11日間掛かる推測です。

初期状態の準備


式(6)を解こうとしても初期条件が与えられなければ数値的に解くことは出来ません。
初期条件で良く使われるのは、1)解析的な解を与える場合と 2)系の固有状態を与える場合の2種類があります。

  1. 解析的な解を与える場合

    この場合、時間発展を考える場合、上記の時間発展方法ではsin関数基底に射影して係数\(c_n(t)\)を求めなければなりません。
    時刻\(t=t_a\)で初期状態\(\psi(x,t_a)\)が解析的に与えられた場合、係数\(c_n(t)\)は、

    として与えられます。

  2. 系の固有状態を与える場合

    この場合は時刻\(t_a\)のポテンシャル中の時間依存しないシュレーディンガー方程式
    \(\displaystyle
    \left[H_0 + V(\mathbf{r},t_c)\right]\phi_n(\mathbf{r})=E’_n\phi_n(\mathbf{r})
    \)

    を実際に解かなければなりません。
    そのためには、時間発展に使う基底と同じ基底\(\varphi_n(\mathbf{r})\)で展開し、対角化するのが賢い方法です。
    すなわち、

    として展開し、ハミルトニアンを行列表示にし、対角化して固有値、固有ベクトルを得ればいいのです。

計算の解析


計算結果を解析するにあたり、ある時刻の固有関数系の基底にどの位存在確率があるかという情報が良く使われます。

初期状態で利用した固有状態の存在確率を調べましょう。
今、求めたいのは時刻\(t_c\)のポテンシャル中の固有関数で展開した時の係数\(a_m(t)\)を求める事です。

ここで、\(\phi_n(\mathbf{r})\)は時刻\(t_c\)のポテンシャル中の固有関数であり、

を満たしています。また、\(\varphi_n(\mathbf{r})\)は数値計算を行う際の基底関数であり、

を満たしています。
任意の時刻\(t\)における波動関数を

として表したいので、\(a_n(t)\)について求めていけば、

という関係式を導くことが出来ます。

数値計算プログラム


プログラムは以下のリンクにあります。
http://slpr.sakura.ne.jp/qp/supplement_data/highprecision_1Dtdsetise.tar.gz
そのままのプログラムでは、以下の設定となっています。
プログラムの利用に関してはhttp://slpr.sakura.ne.jp/qp/about/をお読みください。

解凍すると、4つのファイル(main.f90, integrate.f90,anime.plt,anime_project.plt)が入っています。計算パラメータの設定は、main.f90に入っているので各自調整してください。

繰り返しますが、全て原子単位系です。デフォルトでは、
————————————————–
計算時間の範囲\([ta,tb]\) … \([-5, 20]\)秒
出力時間枚数\(Nt\) … \(100\)枚
1ステップ当りの時間積分の精度\(\mathrm{RKeps}\) … \(10^{-5}\)
1ステップ当りの時間の最小刻み幅\(\mathrm{hmin}\) … \(10^{-4}\)
※理想は最小刻み幅は倍精度ならば\(10^{-14}\)程度が良いと思いますが、時間変化するポテンシャルを計算しようとする時、いつまでも最小刻み幅で進む時間がたまにあります。こういう時間はさっさと抜け出たいので、それを考慮して設定できるようにしています。

位置空間の範囲\([xa,xb]\) … \([-10, 10]\)
出力時の分点数\(Nx\) … 200点
時間積分の精度\(\mathrm{QPeps}\) … \(10^{-6}\)

初期状態を時間依存しないシュレーディンガー方程式を解くか、解析解で与えるかiname … “tise”
inameが”tise”の場合、どの時刻\(t_{ini}\)の固有状態を使うか … \(t_a\)

基底関数の数\(N\) … \(30\)個
計算に用いるスレッド数\(N{\mathrm{cpu}}\) … 1つ
空間積分の分点数を強制的に増やすか否かQPforce … 0

初期状態波動関数… 時間依存しないシュレーディンガー方程式の基底状態
ポテンシャル\(
\begin{eqnarray}
V(x,t)=\left\{
\begin{aligned}
0.5 x^2 (t\lt 0)\\
0.1 x^2 (t\ge 0)
\end{aligned}
\right.
\end{eqnarray}\)
————————————————–
で解くプログラムとなっています。

lapackと、必要があればopemMPを使います。

ifort -mkl -openmp integrate.f90 main.f90
./a.out

もしくは、gfortranでlapackにリンクしてintegrate.f90 main.f90を一緒にコンパイル、実行してください。
実行後、以下のファイルを生成します。

  • wf_***.d
    波動関数の各時刻におけるデータです。時間順にwf_0.d(初期状態), wf_1.d, wf_2.d,…という風に指定された等間隔の時刻で出力します。アウトプットのファイル数はNtで指定した数です。
  • project_ratio.d
    各時刻でのsin関数基底への射影した際の式(2)の係数\(c_n(t)\)の絶対値二乗を出力します。最も高い状態が振幅を持ってしまう場合、基底数が足りていないことを示しています。\( |c_n(t)|^2\)の事です。
  • unitarity.d
    ユニタリー性を確認します。要は数値計算の発散とかで全計算区間内の存在確率が変化していないか見るためのものです。\(\sum_{n=1}^N |c_n(t)|^2\)の事です。
  • timerange.d
    時間グリッドを確認します。物理的な意味は無く、どの時間で計算負荷が大きいかを見るためのものです。
  • timeindex_correspondence.d
    時間インデックス(wf.***.dの***の値と、原子単位系の時刻の変換)を出力します。

もしも初期状態を時間依存しないシュレーディンガー方程式を解いて準備していた場合、以下のファイルを追加で出力します。

  • initial_eigenvalue.d
    初期状態のエネルギー固有値を出力します。式(14)の\(E’_n\)です。
  • project_tdse_to_tisebasis.d
    各時刻における初期状態の固有状態に射影した時の係数の絶対値の二乗を出力します。式(13)の係数\(a_n(t)\)の絶対値二乗である\(|a_n(t)|^2\)を出力します。

波動関数をgnuplot上でアニメーションとして見たい場合、同封してあるファイル”anime.plt”をgnuplot上でloadすれば良いです。

gnuplot> load "anime.plt"

また、gifアニメが欲しい場合、

gnuplot> load "anime_project.plt"

とすれば、波動関数と時間依存しないシュレーディンガー方程式を解いた時の基底関数への係数の絶対値二乗を共に出力します。

実行例


デフォルトの設定で計算を行います。
実行すると

$ ./a.out
0 / 100        1.0000000000
10 / 100        1.0000001064
20 / 100        1.0000001082
30 / 100        1.0000012466
40 / 100        1.0000079960
50 / 100        1.0000082264
60 / 100        1.0000114047
70 / 100        1.0000154242
80 / 100        1.0000154993
90 / 100        1.0000193449
     1.948[CPU sec]
$

という文が出力されます。
左の***/100とは、今計算が100枚中***枚まで終わった、ということを示しており、枚数は変数Ntに対応しています。
右の1に近い数字は、全空間の存在確率を表しており、計算が発散しないかどうかを確認できます。これはルンゲクッタ法の要求精度\(RKeps\)に大きく依存します。
最後のCPUsecは計算に要したCPU時間であり、1コアではおおよそ実時間と同一です。なので、この計算は約2秒で終わる、ということです。ちなみにCPUは1.9GHzの物です。

anime_project.pltによってgifアニメを生成すると、以下のgifアニメが得られます。

その他の出力されるデータは以下の通り。

計算例


ポテンシャルによる反射


※この上図の縦軸は数値計算上の基底である\(\varphi_n(x)\)への射影した時の係数の絶対値二乗です。

ラビ周波数


基底状態と第一励起状態の間の周波数でポテンシャルを揺らしています。

鬼畜なポテンシャル


こんな変なポテンシャルでも計算できます、というデモです。

補遺1

もし仮に、\(U_{m,n}(t)\)が時間によって変化しないのであれば、
\(\displaystyle
\frac{\partial \mathbf{x}}{\partial t}=\mathbf{A}\mathbf{x}
\)

の形の微分方程式であるので、指数関数の解を仮定して\(\mathbf{A}\)を対角化、その固有値固有ベクトルを求めて、線形結合で表現すればokです(線形代数II/連立線形微分方程式等にあるのでそちらを参照してください)。

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変数関数であろうとフィッティングは可能です。
元のデータ
(http://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’\)を結びつける式になっていますので、実験系での角度が分かれば、重心系での角度に焼直せるため、衝突した時の方向も決められるのです(簡単には求められなさそうですが・・・)。

ゼロ点を探す(二分法、挟み撃ち法、Brent法、Newton法、Steffensen法)

ここでは、非線形方程式の解を数値的に求める方法である、
二分法、
挟み撃ち(False position)法、
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回で収束に至っていることが分かります。早いですね。

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

ニュートン法(1、2次元、多次元)

ニュートン法に関するお話です。
”ニュートン法”と呼ばれる方法は
2種類(1.初期値近傍の極小値を求めるニュートン法、2. ゼロ点を求めるためのニュートン法)
があるようですが、ここでは2. 関数のゼロ点を求める方のニュートン法についてのお話です。

  1. まとめ
  2. 1次元のニュートン法
    1. 初期値の推定
    2. 導関数の近似
  3. 1次元ニュートン法のプログラム
  4. 2次元ニュートン法
  5. 2次元ニュートン法のプログラム
  6. 多次元ニュートン法
  7. 高次のニュートン法?
  8. 補)ニュートン法が使われる問題
  9. 参考文献

まとめ

1次元の場合

\(\displaystyle
x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})},~~n=0,1,2\cdots,
\)

ここで、\(x_0\)は解の近傍であること。

2次元の場合

\(
\begin{align}
\left(
\begin{array}{c}
x_{n+1} \\
y_{n+1}
\end{array}
\right)
=
\left(
\begin{array}{c}
x_n \\
y_n
\end{array}
\right)
-\left(
\begin{array}{cc}
f_x(x_n,y_n) & f_y(x_n,y_n) \\
g_x(x_n,y_n) & g_y(x_n,y_n)
\end{array}
\right)^{-1}
\left(
\begin{array}{c}
f(x_n,y_n) \\
g(x_n,y_n)
\end{array}
\right)
\end{align}
\)

\(
\begin{align}
&\left(
\begin{array}{cc}
f_x(x_n,y_n) & f_y(x_n,y_n) \\
g_x(x_n,y_n) & g_y(x_n,y_n)
\end{array}
\right)^{-1} \\
&~~~~=
\frac{1}{f_x(x_n,y_n) g_y(x_n,y_n)-f_y(x_n,y_n)g_x(x_n,y_n)}
\left(
\begin{array}{cc}
g_y(x_n,y_n) & -f_y(x_n,y_n) \\
-g_x(x_n,y_n) & f_x(x_n,y_n)
\end{array}
\right)
\end{align}
\)

ここで、
\(
\begin{align}
f_x(x_n,y_n) &= \left.\frac{\partial f(x,y)}{\partial x}\right|_{x=x_n,~y=y_n} ,~~& f_y(x_n,y_n) = \left.\frac{\partial f(x,y)}{\partial y}\right|_{x=x_n,~y=y_n} \\
g_x(x_n,y_n) &= \left.\frac{\partial g(x,y)}{\partial x}\right|_{x=x_n,~y=y_n} ,~~& g_y(x_n,y_n) = \left.\frac{\partial g(x,y)}{\partial y}\right|_{x=x_n,~y=y_n}
\end{align}
\)

多次元の場合

\(
\mathbf{x}_{n+1}=\mathbf{x}_{n} – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

ここで、
\(
\begin{align}
\mathbf{x}_n=
\left(
\begin{array}{c}
x_1\\
x_2\\
\vdots \\
x_N
\end{array}
\right)_n
,~~
\mathbf{f}_n=
\left(
\begin{array}{c}
f_1(\mathbf{x}_n)\\
f_2(\mathbf{x}_n)\\
\vdots \\
f_N(\mathbf{x}_n)
\end{array}
\right)
,~~
\mathbf{J}_n=
\left.\left(
\begin{array}{cccc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_N} \\
\vdots&\vdots&\ddots&\vdots \\
\frac{\partial f_N}{\partial x_1} & \frac{\partial f_N}{\partial x_2} & \cdots &\frac{\partial f_N}{\partial x_1}
\end{array}
\right)\right|_{\mathbf{x}=\mathbf{x}_n}
\end{align}
\)

※\((~~)_n\)の添え字nはn回目の繰り返しにおけるベクトル\(\mathbf{x}\)を表しています。


スポンサーリンク


ここから本文

1次元のニュートン法


方程式
\(
f(x)=0
\)

を満たす\(x=a\)を求めることを考えます。

まず、この問題をグラフ
\(
y=f(x)
\)

のゼロ点(\(y=0\)となるような\(x=a\))を探す問題だ、という問いに置き換えます。

これを、解近傍で関数をテーラー展開で表し、その展開式のゼロ点を探すことで解を近似することで方程式を解きます。

関数\(f(x)\)が与えられたとき、任意の点\(x’\)周りにおけるテーラー展開は、
\(
f(x) = f(x’)+f'(x’)(x-x’) + O(\Delta^2),~~\Delta=|x-x’|
\)

と書くことが出来ます。もしも\(x\)と\(x’\)が近接していれば、\(\delta=|x-x’| \ll 1\)となるので、\(\Delta^2\)の項は無視できるほど小さくなると期待します。

さて、この式を導く際に用いた条件は、\(x\)と\(x’\)が近くにある、という条件のみです。なので、求めたい解\(a\)と\(a\)の近くの適当な点\(x_0\)を用いれば、
\(
\begin{align}
f(a) &= f(x_0)+f'(x_0)(a-x_0) +O(\Delta^2)~~~…(1)\\
f(x_0) &= f(a)+f'(a)(x_0-a) +O(\Delta^2) ~~~…(2)
\end{align}
\)

の2通りの式を考えることが出来ます。(2)に含まれる\(f'(a)\)は本当の\(a\)が分からないと評価するのが難しいので、(1)を選びます。

(1)の左辺は条件よりゼロです。残りの右辺の\(a\)について解けば、

\(
\begin{align}
0 = f(x_0)+f'(x_0)(a-x_0)+O(\Delta^2) \\
\to a = x_0-\frac{f(x_0)}{f'(x_0)}+O(\Delta^2)
\end{align}
\)


ここで、\(f(x_0)=f'(x_0)=O(\Delta^0)\)を想定しているため、\(f'(x_0)\)の割り算を行っても誤差のオーダーは変化せず、\(\Delta\)に対して2乗のままとなります。

\(O(\Delta^2)\)を無視すると、得られる解は近似値となります。この近似値を\(x_1\)と置くと
\(\displaystyle
x_1 = x_0-\frac{f(x_0)}{f'(x_0)}
\)

得られる解\(x_1\)は\(x_0\)よりも解\(a\)に近いため、\(O(\Delta^2)\)はさらに小さくなります。なので、もう一度計算します。この時の近似値を\(x_2\)と置くと
\(\displaystyle
x_2 = x_1-\frac{f(x_1)}{f'(x_1)}
\)

ここで得られた\(x_2\)は更に\(a\)に近いため、・・・と繰り返していきます。

すなわち、解\(a\)に近い初期値\(x_0\)を与えた時、漸化式
\(\displaystyle
x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})},~~n=0,1,2\cdots,
\)

を繰り返すことで解\(a\)に近づいていく、ということです。
この式から分かる通り、解近傍で導関数がゼロになる場合はゼロ割発生するかもしれないので危険です。

もしもこのステップが無限回繰り返されれば、無限回目のステップの近似値\(x_\infty\)は
\(
x_\infty=a
\)

に収束します。

以上のように導関数を用いて方程式の解を求める方法は、ニュートン法(またはニュートン=ラフソン法)と呼ばれます。
この方法は、解に近いに所の値とその導関数が分かりさえすればいいので、この方法は1次元、多次元でも使用する事ができ、また複素関数でさえ使用することが出来ます。
さらに、誤差が\(O(\Delta^2)\)で減衰していくため、一度繰り返すごとに厳密解に一致する桁数が2倍になっていきます。

しかし問題があり、
1. 初期値を解の近傍に取らないと失敗する
2. 導関数が分からないと使えない
という大きな問題を抱えています。

また、解近傍でテーラー展開可能でなければならず、解近傍で導関数がゼロになる点があってはいけません。

  • 初期値の推定


    ニュートン法は解の近傍で、テーラー展開可能な領域から始めなければ用いることは出来ません。
    これは工夫をするしかありません。

    解近傍の初期値を推定する方法として、以下の2つの方法が良く取られます。
    1. 解の範囲を限定し、その近傍からニュートン法を行う(1次元のみ)。
    2. 適当な項を落とせば方程式が解けてしまう場合、その解を出発点とし、残りの項を摂動的に加えながら逐次解く。

    1. 解の範囲を限定する

      これは1次元の場合に使える方法です。
      単純な考え方で、例えば広い区間\([x_a,x_b]\)にゼロ点があることだけが分かっているとします。
      この区間を\(N+1\)分割し、
      \(
      x^{(0)}=x_a,x^{(1)},\cdots,x^{(k)},x^{(k+1)},\cdots, x^{(N-1)},x^{(N)}=x_b
      \)

      とします。もしも\(x^{(k)}\)と\(x^{(k+1)}\)の符号が違う場合、この区間でゼロを横切っているということになりますので、この点からニュートン法を始めれば良いのです。

      問題は、解の範囲を狭めるために関数を\(N+1\)回評価しなければならない、という点で計算コストがかかってしまう、という点です。

    2. 既知の解から摂動的に逐次解く

      これは若干特殊な方法で、あてずっぽうな上の考えとは違います。
      本来解きたい式を
      \(
      f(x)=f_0(x)+g(x)
      \)

      と分けます。ここで、
      \(
      f_0(x)=0
      \)

      の解は簡単に解けて、\(x=a_0\)が関数\(f_0(x)\)のゼロ点であるとします。この\(a_0\)を初期値とし、\(\lambda\)を十分小さい値にとり
      \(
      f_0(x)+\lambda_1 g(x)=0
      \)

      を解きます。ここで、\(\lambda_m\)は\(0\lt\lambda_1\lt\lambda_2\cdots \lt\lambda_{M-1} \lt \lambda_M=1\)を満たすとします。
      この方程式の解\(a_1\)が得られたら、この解を初期値とし、
      \(
      f_0(x)+\lambda_2 g(x)=0
      \)

      を解きに行き、これを\(\lambda_M=1\)まで続けます。
      こうすることにより、最後の答え\(a_M\)がもともと解きたかった式の答えとなっているのです。

      関数の一回の評価に時間がかかってしまう場合、この解き方の方が早く解けます。
      しかし、関数\(g(x)\)によってゼロ点の位置が増えたり、消えたりしてしまう時にこの方法は使うことが出来ず、連続的に解が動いてくれないと使うことが出来ない点に注意してください。

  • 導関数の近似


    導関数を求める簡単な方法は、差分を使うことです。すなわち、導関数\(f'(x)\)を差分に置き換えます。
    \(\displaystyle
    f'(x)\approx \frac{f(x+h)-f(x)}{h} + O(h)
    \)

    すると、漸化式は
    \(\displaystyle
    x_{n+1}=x_{n}-h\frac{f(x_{n})}{f(x_n+h)-f(x_n)},~~n=0,1,2\cdots,
    \)

    と得られます。
    この方法は導関数を必要とせず、実装が簡単であるという意味で有用な式です。

    導関数を求める際の前進差分の刻み幅\(h\)は、計算機の扱える丸め誤差\(\epsilon\)(単精度ならば約\(\epsilon=10^{-8}\)、倍精度ならば約\(\epsilon=10^{-16}\))の\(1/2\)乗の2倍、つまり
    \(
    \displaystyle h\sim 2\sqrt{\epsilon}
    \)

    とすると良いです。

    この刻み幅の導出は折りたたんでおきますが、以下のように導出することが出来ます。

参考文献[1][3][5]

1次元ニュートン法のプログラム


以下のプログラムは、1次元のニュートン法のプログラムで、
方程式
\(
x^2-4 = 0
\)

を初期値\(x_0=3\)から始めて解くプログラムです。

このコードを動かすと、ステップを繰り返すごとに急激に解に近づいていっていることが分かります。
(デフォルトではコメントアウトしてあります。)

$ gfortran main.f90
$ ./a.out
    1       2.1666666616021075
    2       2.0064102565311974
    3       2.0000102400381690
    4       2.0000000000262466
    5       2.0000000000000000
   2.0000000000000000       ---result
$

上には出力していませんが、初期値の評価もしています。
なので、ニュートン法は6回呼び出されたことになります。
ニュートン法1回当たり関数は2回呼び出されますので、12回関数が呼び出されたことになります。

例えば比較対象として、二分法を考えますと、1回目に関数は3回,そのほかで1回ずつ呼び出されます。
二分法の解は1ステップ当たり区間の半分になりますので、16桁一致するまでには50回ほど繰り返す必要があります。そのため、52回位関数が呼び出されます。

ニュートン法の方が圧倒的に早いことが分かります。

複素関数の1次元ニュートン法のプログラム
$ gfortran main.f90
$ ./a.out
                    1 ( -1.2933333304320405     , 0.72000000272341036     )
                    2 (-0.78207788228079722     , 0.60930726266833468     )
                    3 (-0.43844290526861285     , 0.73503785880633421     )
                    4 (-0.50851135934742864     , 0.89043164725052448     )
                    5 (-0.50009896196621662     , 0.86666395460253720     )
                    6 (-0.49999991066297655     , 0.86602581130981315     )
                    7 (-0.49999999999985084     , 0.86602540378453374     )
                    8 (-0.50000000000000000     , 0.86602540378443860     )
 (-0.50000000000000000     , 0.86602540378443860     )  ----result
$

ソースコードはこちら。

2次元のニュートン法


続いて方程式
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f(x,y)&=0,\\
g(x,y)&=0
\end{aligned}
\right.
\end{eqnarray}
\)

を満たす\(x=a,y=b\)を見つけることを考えます。

1次元の時と同様に、任意の点\(x’,y’\)周りにおけるテーラー展開は、
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f(x,y) &= f(x’,y’)+f_x(x’,y’)\cdot(x-x’)+f_y(x’,y’)\cdot(y-y’)+ O(\varepsilon^2)\\
g(x,y) &= g(x’,y’)+g_y(x’,y’)\cdot(x-x’)+g_y(x’,y’)\cdot(y-y’)+ O(\varepsilon^2)
\end{aligned}
\right.
\end{eqnarray}
\)

ここで、
\(\displaystyle \left. f_x(x’,y’)=\frac{\partial f(x,y)}{\partial x}\right|_{x=x’,y=y’}\)
\(\displaystyle \left. f_y(x’,y’)=\frac{\partial f(x,y)}{\partial y}\right|_{x=x’,y=y’}\)
\(\displaystyle \Delta x=|x-x’|=O(\varepsilon^1),~~\Delta y=|y-y’|=O(\varepsilon^1)\)
と書きました。

1次元の時と同様に初期値\((x_0,y_0)\)が解\((a,b)\)に近ければ、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
0 &= f(x_0,y_0)+f_x(x_0,y_0)\cdot(a-x_0)+f_y(x_0,y_0)\cdot(b-y_0) +O(\varepsilon^2)\\
0 &= g(x_0,y_0)+g_x(x_0,y_0)\cdot(a-x_0)+g_y(x_0,y_0)\cdot(b-y_0) +O(\varepsilon^2)
\end{aligned}
\right.
\end{eqnarray}
\)

を同時に満たす\((a,b)\)を探せばよい、ということになります。

行列形式で書けば
\(
\begin{align}
\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)
\left(
\begin{array}{c}
a-x_0 \\
b-y_0
\end{array}
\right)
=
-\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)+O(\varepsilon^2)
\end{align}
\)

と書けます。\(a,b\)について解けば
\(
\begin{align}
\left(
\begin{array}{c}
a \\
b
\end{array}
\right)
=
\left(
\begin{array}{c}
x_0 \\
y_0
\end{array}
\right)
-\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)^{-1}
\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)+O(\varepsilon^2)
\end{align}
\)

となります。丁度一次元のニュートン法と似た形になりました。

右辺の
\(
\begin{align}
\mathbf{J}_0=\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)
\end{align}
\)

は\(x=x_0,y=y_0\)のヤコビアンです。また、
\(
\mathbf{a}=\left(
\begin{array}{c}
a \\
b
\end{array}
\right),~~
\begin{align}
\mathbf{x}_0=\left(
\begin{array}{c}
x_0 \\
y_0
\end{array}
\right),~~
\mathbf{f}_0=\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)
\end{align}
\)

と表記すれば、2次元のニュートン法は
\(
\mathbf{a}=\mathbf{x}_0 – \mathbf{J}_0^{-1}\mathbf{f}_0+O(\varepsilon^2)
\)

と書くことが出来ます。
近似として、\(O(\varepsilon^2)\)の項を無視すれば
\(
\mathbf{x_1}=\mathbf{x}_0 – \mathbf{J}_0^{-1}\mathbf{f}_0
\)

と書けます。繰り返せば、漸化式
\(
\mathbf{x}_{n+1}=\mathbf{x}_{n} – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

を解いていけば良い、ということになります。

ニュートン法を実行するにあたり、ヤコビアンの逆行列を解くことが一番の問題となります。

良く知られているように、2次元の逆行列は、計算コストもさほどかからずに解けてしまい、
\(
\begin{align}
\mathbf{J}_n^{-1}=
\frac{1}{f_x(x_n,y_n) g_y(x_n,y_n)-f_y(x_n,y_n)g_x(x_n,y_n)}
\left(
\begin{array}{cc}
g_y(x_n,y_n) & -f_y(x_n,y_n) \\
-g_x(x_n,y_n) & f_x(x_n,y_n)
\end{array}
\right)
\end{align}
\)

として求める事が出来ます。

参考文献[2],[6]

スポンサーリンク

2次元ニュートン法のプログラム


以下のプログラムは、2次元のニュートン法のプログラムで、
方程式
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f_1(x,y)&=x^2+y^2-1 = 0 \\
f_2(x,y)&=y-x^3 = 0
\end{aligned}
\right.
\end{eqnarray}
\)

を初期値\(x_0=2,~y_0=1\)から始めて解くプログラムです。

このコードを動かすと、ステップを繰り返すごとに解に近づいていっていることが分かります。

$ gfortran main.f90
$ ./a.out
    1       1.3571428532359113       0.2857142813732352
    2       0.9844126813492445       0.4401111876910193
    3       0.8485699953676349       0.5590406123395349
    4       0.8265084709901402       0.5633730814110012
    5       0.8260315915076198       0.5636240770681390
    6       0.8260313576542442       0.5636241621612319
    7       0.8260313576541870       0.5636241621612585
  0.82603135765418700       0.56362416216125855       ---result
$

この問題の答えは\(f_1=0,~f_2=0\)の陰関数の交点です。
図示すれば、下図の紫色の経路を辿って解に収束していくことが分かります。

ひとつ、注目しておかなければならない点として、解に近いからと言って、その解に収束するとは限らないということに注意してください。
図中の水色の線は初期値こそ右上の解に近いですが、実際に計算してみるともう一つの解に収束していることが分かります。

適当な初期値を用意し、ニュートン法が収束する点を調べるという問題はニュートン写像と呼ばれる分野です。
力学系とかそういう言葉が出てきます。
私も過去に複素関数でやったことがありますので、リンクだけ載せておきます。
ニュートン写像に現れる綺麗な画像 -シキノート

2次元の複素関数のニュートン法

置いておきます。

多次元のニュートン法

さて、この形まで持ってくると多次元への拡張が簡単にできます。
多次元のニュートン法は
\(
\mathbf{x}_{n+1}=\mathbf{x}_n – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

の形がそのまま使えます。
ここで、
\(
\begin{align}
\mathbf{x}_n=
\left(
\begin{array}{c}
x_1\\
x_2\\
\vdots \\
x_N
\end{array}
\right)_n
,~~
\mathbf{f}_n=
\left(
\begin{array}{c}
f_1(\mathbf{x}_n)\\
f_2(\mathbf{x}_n)\\
\vdots \\
f_N(\mathbf{x}_n)
\end{array}
\right)
,~~
\mathbf{J}_n=
\left.\left(
\begin{array}{cccc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_N} \\
\vdots&\vdots&\ddots&\vdots \\
\frac{\partial f_N}{\partial x_1} & \frac{\partial f_N}{\partial x_2} & \cdots &\frac{\partial f_N}{\partial x_1}
\end{array}
\right)\right|_{\mathbf{x}=\mathbf{x}_n}
\end{align}
\)

です。

問題となるのは(ヤコビアンの逆行列)×(f)の計算です。
この計算を数値計算的に何とかするべく改良が重ねられてきたのが準ニュートン法(本稿では書きません)と呼ばれるアルゴリズムです。

逆行列を数値的に計算することは宜しくないですので、
\(
\begin{align}
\mathbf{r}&=\mathbf{J}_n^{-1}\mathbf{f}_n \\
\to \mathbf{J}_n \mathbf{r} = \mathbf{f}_n
\end{align}
\)

という方程式を解く問題に帰着させます。これが解ければ、ニュートン法の次のステップを
\(
\mathbf{x}_{n+1}=\mathbf{x}_n – \mathbf{r}
\)

として求める事が出来るのです。

高次のニュートン法?


1次元のニュートン法は
\(\displaystyle
a = x_0-\frac{f(x_0)}{f'(x_0)}+O(\Delta^2)
\)

の形で求めていきますが、二階微分が分かっていればもっと収束が早くなりそうです。
実際、
\(\displaystyle
a = x_0-\frac{f^{\prime}(x_0)}{f^{\prime\prime}(x_0)}\pm\sqrt{{f’}^2(x_0)-2f(x_0)f^{\prime\prime}(x_0)}+O(\Delta^3)
\)

と解けます。しかし、2次関数ですので
解が見付からないかもしれない
という可能性があります。1次のニュートン法であれば、導関数がゼロでない限り、必ず\(y=0\)と交わる点が見付かります。しかし、2次で解に近づけると交点が見つからないことがあるかもしれません。

なので、1次のニュートン法を繰り返して使う方が(2次と比べたら)安定な解法なのです。

参考ページ

[1]1 Newton 法:壱変数の場合
https://www.astr.tohoku.ac.jp/~chinone/Planck/Planck-node6.html
[2]2 Newton 法:弐変数の場合
https://www.astr.tohoku.ac.jp/~chinone/Planck/Planck-node7.html
[3]Newton法
http://www.misojiro.t.u-tokyo.ac.jp/~murota/lect-suchi/newton130805.pdf
[4]アルゴリズムによる誤差
http://www.aoni.waseda.jp/ykagawa/chap2html/node9.html
[5]Newton 法による方程式の近似解法
http://www.math.u-ryukyu.ac.jp/~suga/C/2004/7/node9.html
[6]多次元のニュートン・ラフソン法
http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%C2%BF%BC%A1%B8%B5%A4%CE%A5%CB%A5%E5%A1%BC%A5%C8%A5%F3%A1%A6%A5%E9%A5%D5%A5%BD%A5%F3%CB%A1
[7]山本 有作, 行列計算における高速アルゴリズム
http://www.cms-initiative.jp/ja/events/0620yamamoto.pdf
[8]中田 和秀, 大規模線形方程式を解くためのクリロフ部分空間法の前処理
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1288-6.pdf

微分方程式を解く際の計算量の工夫

運動方程式を解いて、ある特定の時刻の値が欲しいとします。
特定の時刻がいつなのか、あらかじめ分かっている場合、与えられた初期条件から特定の時刻まで計算し、データを貯めておいて後で使えば良いです。

特定の時刻が分からない時に計算量の問題が生じます。
例えば、時刻\(t=0\)の初期条件が与えられた時、時刻\(t_1=100,t_2=101, t_3=102\)の時刻が知りたいとします。ここで、時刻\(t_2\)の値は実は\(t_1\)の関数の値を評価してから初めて決まるとします。

この場合の計算方法としては、

  1. 毎回、時刻\(t_0\)から初めて\(t_1,t_2,t_3\)まで計算
  2. あらかじめ微分方程式を一度細かく計算しておき、欲しい時刻\(t=t’\)を補間で求める
  3. あらかじめ微分方程式を一度細かく計算しておき、欲しい時刻\(t=t’\)に近い点から計算

だと思います。最も簡単に実装できるのは1の方法で、欲しい時刻が数点であれば計算時間は問題にならないほど少なくなるでしょう。2の方法は精度を気にする場合には使用できません。ここでは3番目の方法を考えます。

計算方法を図示したものが以下のものです。

初期条件からいちいち出発して\(t=t’\)に到達するのではなく、あらかじめ計算していた参照点から出発して計算する、ということです。
ここで問題なのが、もし横軸が時刻だった場合、時間を逆向きに遡って良いか?という問題です。

これは物理の問題で言えば、解こうとしている問題は時間可逆か?ということで、言い換えれば系のエネルギーは保存しているか?ということです。もしも抵抗力など、エネルギーが散逸する場合は、いくら参照点が近いからといって逆向きに進んではいけないということです。

ここでは、微分方程式を解く際に戻ることはせず、必ず求めたい時刻よりも前の参照点からスタートする、というプログラムを作ります。

以下のように実装できます。