「数学」カテゴリーアーカイブ

古典的な一次元波動方程式のグリーン関数

本稿のpdfはこちらをどうぞ
https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/08/Green_1dclassical5.pdf

問題


古典的な一次元波動方程式で、非斉次項を持つ問題

を解くことを考えます。この形は、有限範囲内しか相互作用を引き起こさない物に十分遠方から波が入射してきた場合を考える際に現れたりします。

この形の方程式が出てくる状況を考えてみましょう。

の問題を考えます。移項して

ですので、グリーン関数を用いると

となります。ここで、

を満たす関数(一般解)であり、グリーン関数を
で与えられます。ここを満たす関数(特殊解)としました。
式(4)の右辺に解\(f(x,t)\)が含まれていますが、このままの形でとどめておきます。この形にしておくと、ほかの式に\(f(x,t)\)をダイレクトに代入できたり、近似を使うときに便利な形です。

(式(7)はありません)

本稿の目的は、\(\hat{D}\)が\(\hat{D}=\frac{\partial^2}{\partial t^2}-a\frac{\partial^2}{\partial x^2}\)で与えられている場合に、\(G(x,t; x’,t’)\)の具体的な形を導出することです。

導出


今、解きたい問題は、

です。デルタ関数をそのまま扱う場合、積分が絡んでこないと扱うのが難しいです。そうでなければ、フーリエ変換を用いて波数・周波数空間で解いていくのが良いでしょう。グリーン関数や、デルタ関数のフーリエ変換を考えると

となります。ここで、基底関数が\(e^{ikx}, e^{-i\omega t}\)と符号が異なることに注意しましょう。
また、積分区間はいつも負の無限大から正の無限大であり、いちいち書くのは面倒なので省略しています。
フーリエ変換として考えるのではなく、基底関数へ射影した時の係数と考えましょう。

実際に式(8)に代入すると

となります。

ここで\(e^{i[k(x-x’)-\omega(t-t’)]}\)はゼロにならないので、その前の項[ ]内がゼロにならなければなりません。よって

が導けます。よって式(9)に代入して

が計算できれば位置と時間のグリーン関数が得られます。

さて、この積分を計算してみましょう。
右辺の\(e^{i\omega t}/\omega\)の形を持つ積分は複素関数論において頻出する積分です。この積分の特徴として、\(\omega\)に渡る積分を実軸上で実行する際に、たまたま\(\omega=\pm \sqrt{a}k\)に等しい点を通ってしまうため、実軸上では被積分関数が発散する、という特徴があります。この発散する点は極と呼ばれており、積分を行う際に積分経路をうまく選んで極を迂回しなければなりません。しかし、厄介なことに極を上に回るか下に回るかで積分結果が変わってしまいます。

もし、天才でしたら当たり前のように適切な積分経路をすぐに思いつくでしょうが、我々凡人には具体的に試してみるくらいしか思いつきません。ですのでちょっとズルしながら因果律を満たす解を探しましょう。

今、我々はグリーン関数について考えています。ということは、計算した結果の境界条件は、因果律を満たすような結果が欲しいわけです。
つまり、グリーン関数が得られた場合、$\theta(t-t’)$の形が出てくるであろうことは予想できます。そのため、これを考慮して極を回る方向を考えてグリーン関数を求めていきましょう。

式(14b)の第一項

について考えます。ヘヴィサイド関数を積分表記で表すと、

と書けることを使います。ここで\(i\)は虚数単位です。もう嬉しいですね。かなり似通った形であることが分かるでしょう。
式(16)の形に持っていくためには変数変換を行えば良さそうです。つまり、極は上に回れば因果律を満たすグリーン関数が得られそうです。

では具体的に計算していきましょう。変数変換\(y=\omega+\sqrt{a}k\)を使って、

となります。分母に\(+i\varepsilon\)を加えるようにして

と求めることができます。

続いて式(14b)の第二項

についても同様に計算していきましょう。変数変換\(y=-\omega+\sqrt{a}k\)を使って、

ヘヴィサイド関数を考えると、

の形が使えそうです(式(16)の\(\omega\to-\omega\)の変数変換ですね)。ですので、

と求められます。

式(18c)と式(22c)より、式(14b)に代入すると

となります。

さて、\(\frac{\sin x}{x}\)の形はsinc関数と呼ばれているほど面白い性質を持つ関数です。そのフーリエ変換は矩形になることが知られています。その計算は実際に
\begin{align}
\int_{-\infty}^{\infty}\frac{\sin(k)}{k}e^{ikx}\frac{dk}{2\pi}=\frac{1}{2}\Bigl[\theta(x+1)-\theta(x-1)\Bigr]
\end{align}
となります。

立ち戻り、変数変換\(\kappa=\sqrt{a}k(t-t’)\)を考えれば、

と、グリーン関数が求められました。グリーン関数は\(x-x’, t-t’\)の形をしていることが分かります。そのため、

と書いても良いことが分かるでしょう。ここで、

です。

まとめと補足


まとめ


まとめますと、偏微分方程式

の因果律を満たす解は、

となります。以上から、波動方程式

の因果律を満たす一般解は

となります。

補足


補足します。因果律を満たすグリーン関数である、ということが暗に意味される場合、積分(29a)の時間に渡る積分はグリーン関数に含まれるヘヴィサイド関数を先に利用して、

と積分範囲を組み込んで表記する場合があります。実際の例が式(29c)であり、その時間積分の積分範囲のようになることを見越している、ということです。

少しグリーン関数の物理的な意味を考えてみましょう。
グリーン関数は、ある特定の時刻、ある特定の位置でデルタ関数という非常に特殊な衝撃が生じた結果を記述します。
言い換えれば、偏微分方程式で表されている系が、デルタ関数によってどのような応答を起こすか?という関数です。

式(27b)を見てみると、1つのヘヴィサイド関数\(\theta(t-t’)\)と1つの矩形関数\(\text{rect}\left(\frac{x-x’}{2\sqrt{a}(t-t’)}\right)\)が使われています。
\(\theta(t-t’)\)の部分は、まさに因果律を示しており、衝撃が加わる前の時刻\((t-t’ \lt 0)\)においては系の応答は無いよ、ということを表しています。至極全うでしょう。

矩形関数の部分は、衝撃が生じた後にその衝撃は波及する範囲を表しています。具体的に、時刻\(t=t’\), 位置\(x=x’\)において衝撃が生じた場合、矩形関数の中身が値を持つ範囲は\([x’-\sqrt{a}(t-t’), x’+\sqrt{a}(t-t’)]\)の範囲です。つまり、衝撃が時々刻々と広がっていることを示しています。グリーン関数を図示すると、このような感じです。

衝撃が加わる\(t=0\)を境にグリーン関数が値を持ち始め、それが時間とともに広がっていくのが分かると思います。

ついでにグリーン関数の広がる速度を計算してみましょう。

矩形関数の中身が値を持つ範囲\(x\)は、時間\(\Delta t=t-t’\)の間に、距離\(x-x’=\pm \sqrt{a}(t-t’)\)だけ波及しますから、衝撃の結果が速度\(\pm\sqrt{a}\)で広がっていることが分かります。
これは、古典的な波動方程式でよく知られているように、系の速度が

の系の速度は\(\pm\sqrt{a}\)である、という事実と一致します。
つまり、衝撃が加わった時刻以外では、系そのものの性質しか存在しませんので、衝撃が波及していく速度が一致することは何ら不思議ではありません。むしろ、一致しなければなりません。

デルタ関数、ヘヴィサイド関数に関係する数式

本稿のpdfはこちら
https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/08/Green_related_functions1.pdf

定義


ヘヴィサイド関数

符号関数

矩形関数

デルタ関数

デルタ関数\(\delta(x)\)は、

を満たす関数です。この定義式を満たす場合、デルタ関数は以下の性質を持つことが分かります。

フーリエ関数

性質、関係式


積分表示

定積分が関係する場合


デルタ関数

ヘヴィサイド関数


sinc関数に関係する式

畳み込み積分


デジタル信号処理と偏微分方程式の関係

問題提起


デジタル信号で良く以下の図のように入力と出力が説明されます。

ここで疑問が生まれました。それは

  • 時間だけの関数であるならば、通過前後の位置はどこに記述されている?
  • 出力が、入力とシステムの畳み込み積分である根拠は?

です。通過位置については、どこにも位置が現れていないのです。また、畳み込みは確かに正しそうですが、どこから畳み込みが現れたのか根拠が分かりません。
これの解決を目指します。

本稿の結論は、散乱問題の考え方と同じで良い、です。つまり、

  • 位置と時間は、偏微分方程式の分散関係によって結び付いているので、時間又は位置について分かればもう片方は自動的に来まる。
  • 初期状態から終状態が生まれる際にグリーン関数を用いて書くことができるので、それが畳み込みの形となり、出力となる。

ということです。

序論


デジタル信号処理に現れる入力・出力を表す関数は何か?


波は偏微分方程式として記述されることが多く、基本的に波は位置と時間の関数として書かれます。
しかし、デジタル信号処理の分野では位置に関する表記を消去して時間のみの関数(時間信号)について解くだけで良し、とされます。
波はそもそも偏微分方程式ですので、何らかの仮定をおいて位置は考えなくて良し、としているはずですが、その仮定とは何なのでしょうか。
調べても『当たり前』とされている部分のようで、数式で示されている資料を見つけることができませんでした。

この疑問を詳しく説明します。
デジタル信号処理では、波と相互作用する何かを”システム”という言葉で表現し、その”システム”を特徴づける”インパルス応答”やそのフーリエ変換である”周波数特性”として扱います。
そして、”入力”と”システム”の相互作用を畳み込み積分した結果を”出力”として表します。
ここで私が疑問に思ったのは、

  • 「インパルス応答を用いて畳み込み、は妥当そうだが、それはなぜokなのか?」
  • 「現実には”システム”を通過する前後の位置という情報があるはずなのに、それはどこに行った?」

という2点です。

畳み込み積分が出てくる数学上の話はグリーン関数が想像されますが、やはり検索しても見当たらないので自分で導出を試みました。
いろいろ悩みましたが、以降の導出は私の導き、納得できた結論です。考え方は散乱問題と同じです。
つまり、畳み込みが出てくるということは非斉次方程式が出てきて、グリーン関数を用いて解が表される方程式があるということです。

導出


関係性の導出


波動方程式

を解くことを考えます。ここで、\(\hat{D}\)は位置\(x\)と時間\(t\)に関する微分演算子です。具体的に例えば
\begin{align}
\hat{D}=\frac{\partial^2}{\partial t^2}-a\frac{\partial^2}{\partial x^2}
\end{align}

みたいな感じです。相互作用を表現する\(s(x,t)\)の持つ性質として
\begin{align}\displaystyle
s(x,t)\Bigr|_{x\gt |a|}=0
\end{align}
を仮定します。

この仮定の下で、初期状態は\(x\lt -|a|\)の領域のみで波形が与えられると考えます。つまり、初期状態\(f_0(x,t)\)は

を満たす状態として与えられます。この条件下で、式(1)の解を

の形で探します。\(g(x,t)\)は未知の、これから求めたい関数です。
式(1)に代入すると\(g(x,t)\)について

を満たす関数として書けます。式(4)の解\(g(x,t)\)は

とグリーン関数を用いて書くことができます。ここで、\(g_0(x,t)\)は
\begin{align}
\Bigl[\hat{D}(x,t)+s(x,t)\Bigr]g_0(x,t)=0
\end{align}
を満たす解であり、\(G(x,t;x’,t’)\)は

を満たす解です。補足しますが、右辺を見ると\(x-x’, t-t’\)の形でしか存在しません。これはつまり、左辺のグリーン関数も\(x-x’, t-t’\)の関数になっていることを意味します。つまり、
\begin{align}
G(x,t;x’,t’)=G(x-x’,t-t’)
\end{align}
という形になっており、グリーン関数は4変数関数ではなく、2変数関数であることが分かります。どちらの表記でも構わないので適宜用いることにします。
もし初期状態が存在しない、つまり\(f_0(x,t)=0\)ならば\(g(x,t)=0\)が想定されるため、

でなければなりません。よって、

と書けます。ここで\(H(x,t;x’,t’)\)は

を満たす関数で、俗にいうインパルス応答を記述します。ここでもGと同様に
\begin{align}
H(x,t;x’,t’)=H(x-x’,t-t’)
\end{align}
と書くことができます。以上から、式(1)の解は

と書けます。ここで解の意味を考えます。右辺第一項は初期状態(入力)であり、右辺第二項は、初期状態が\(s(x,t)\)と相互作用した結果生じる相互作用項(出力)を意味します。

仮に初期状態が

とかけていた場合を考えます。式(10)の右辺第二項を計算すると

とまとめられます。ここで、

と定義しました。以上から、解は

となります。良く知られているように、波数空間(もしくは、初期状態を角周波数\(\omega\)の変数とするならば角周波数空間)で相互作用の結果は、相互作用(インパルス応答)のフーリエ変換 \(H(k,\omega(k))\)と初期状態のフーリエ変換との積となります。

位相速度、群速度とは!?調べてみました!!(パロディ)

近頃SNSやインターネットでよく見る言葉、「位相速度」「群速度」

元ネタや原理を分かった気になって使っている方、多いのではないでしょうか。
せっかくなので調べてみました!!!

  1. そもそも「位相速度」「群速度」って何?
  2. 学会騒然。ある波数近傍だけ値を持つ場合にテーラー展開!?
  3. 位相速度はどうやって!?驚きの仮定
  4. 簡単だと思いました?注意です!
  5. いろいろなところに群速度に関する議論が!
  6. 最後に
  7. Note

そもそも「位相速度」「群速度」って何?


なにはともあれwikipedia!!
関連するキーワードをササッと知るには強い味方です。
「位相速度」「群速度」を調べるとこのように説明されていました。

群速度(ぐんそくど、英: group velocity)とは、複数の波を重ね合わせた時にその全体(波束)が移動する速度のことである。

群速度 -wikipedia

位相速度(いそうそくど、英語: phase velocity)は、位相、すなわち波の山や谷の特定の位置が移動する速度のことである。

位相速度 -wikipedia

どうやら、波束波の速度を表現する際に都合が良い言葉なようです。

波束とは何なのでしょう…?また変な言葉が出てきました。wikipediaの出番です!!

波束(はそく、英: wave packet, wave train)は、局所的に存在する波うち/波動であり、移動する1個の波動の塊のようにふるまう。

波束 -wikipedia

うーん、つまり一回だけ波打って、ある程度まとまった波みたいなものですかね。
水面に水が一滴だけ落ちたときに発生する波みたいなものでしょう。

具体的にどういう意味なのでしょうか。気になりますよね!

皆さんも群速度、位相速度の導出が気になったら「tweet」、お願いします!

学会騒然。ある波数近傍だけ値を持つ場合にテーラー展開!?


角周波数が波数の関数として書けている場合、分散関係があるといいます。
つまり、角周波数\(\omega\)が波数\(k\)に依存して\(\omega=\omega(k)\)と書けている場合、波形のほとんどは

とように書けます!!
本来、位置\(x\)と時間\(t\)で書ける関数は、\(k\)と\(w\)の二重積分で表現されなければなりませんが、分散関係のおかげで\(\omega\)が強制的に決まってしまうので、一重積分で良くなるのですね!
うーん、とても嬉しい!!!!

…でも、あまりに一般的過ぎて分かりませんね。
しかし!!
応用で重要なのは、ゆっくり振動する包絡線早く振動する波掛け算で表されていることが多いのです。

つまり、これから考える波が持つ波数は、\(k= k_0\)の周りしか値を持たないのです。ここで、早く振動する波の波数を\(k_0\)と置きました。

\(k_0=5\)の場合、これから考えていく\(f(k)\)はこんな感じの特徴を持っています!ババン!

こんな特徴を持つのならば、\(k=k_0\)以外ではどうせ\(f(k)\to 0\)となるので、考えてもしょうがないっていう近似が使えそうですよね。

そ こ で !

波数\(k=k_0\)周りでいろいろテーラー展開して使用していきます。つまり、\(k=k_0\)の近傍さえ合っていれば、\(k_0\)以外でどんな振る舞いになろうとも\(f(k)\)がゼロになるので大丈夫でしょう!ということです。

早速、分散関係についてテーラー展開を行うと

となります。式を簡単に書くために\(\omega'(k_0)\equiv \frac{d\omega}{dk}\bigr|_{k=k_0}\)と置きました。

計算を進めると、

より、

のように書けます。つまり、式(4)が言っていることは、
時刻\(t\)の波形とは、初期状態\(t = 0\)の波形 \(f(x,0)\)を\(\omega'(k_0)t\)だけ平行移動させた波に、時間だけに依存する位相に関する変化分 \(e^{−i[\omega(k_0)−k_0 \omega'(k_0)]t}\)が掛け合わされたもの、と表すことができる。と言っています。

以上から、初期状態の波が時刻\(t\)に至るまでの速度は\(\omega'(k_0)\)で決まることが分かり、これは
群速度
と呼ばれます!!

位相速度はどうやって!?驚きの仮定


さて、ここまでで群速度が求まりましたが、位相速度が求められていません。
そこで、もう少し具体的な波の形を定義しましょう。
波が、

と書けている場合を想定します。ここで、振幅\(A(x,t)\)は実数値関数で波の包絡線を表しています。また\(A(x,t)\)の変化は \(k_0\) に比べて非常にゆっくり変化しているものとします。
つまり、

今まで考えてきた仮定が全部使える

んですね!

式(5)に代入しますと

となります。つまり、振幅の変化はいつも初期状態の単なる平行移動として表現することができ、位相の変化も簡単に表せられます。

そこで、振幅が平行移動する速度を群速度\(v_g\)、位相が平行移動していく速度を位相速度 \(v_p\) として定義すると、 式(6c) よりそれぞれ

と定義してしまうのが良さそうです!

簡単だと思いました?注意です!


注意しなければならないのは、群速度と位相速度の概念は、

振幅と位相を記述する部分がはっきりと区別できるような状況でなければ、定義できない

ということです!振幅部が早く振動していたら、これまでの議論が使えなくなることに注意しましょう。
その場合、群速度と位相速度の概念が変わってしまうので、意味をなさなくなってしまいます。

いろいろなところに群速度に関する議論が!


様々なところにリンクがあります。やはり皆、群速度の概念を理解するのに苦労しているようですね。

群速度の謎 -平林 浩一, (C) 2006

「yam@広島大」物理化学Monographシリーズ ”3. 物体の速度と物質波の速度 -E=hνの本質の理解-”

最後に


いかがでしたか?気に入ったなら是非、記事の下にある「like!」ボタンや「tweet」して皆さんに共有して下さいね!
それでは皆さん、良い物理ライフを!

Note


本記事の大枠の構成はいかがでしたか? -ニコニコ大百科を参照しています。

陰関数の等高線を数値計算で求める

まとめ


関数

を求める問題を考えます。
\(f(x_0, y_0)=0\)を満たす\((x_0, y_0)\)が与えられている場合、\((x_1\equiv x_0+\Delta x, y_1\equiv y_0+\Delta y)\)を求めることを考えていきます。
まず、

の変換を行い、初期値\((x_0, y_0)\)を\((x’_0, y’_0)\)に変換します。
ここで、\(\theta\)は定数で

です。\((x’_0, y’_0)\)を新たな初期値として微分方程式

を解き、\(x_0′\)から微小距離\(h\)だけ発展させた\(x_1’\equiv x_0′ + h\)とその値\(y_1’\equiv y_0′ + \Delta y\)を求めます。
その後、逆変換

を経ることで\(f(x_1,y_1)=0\)を満たす\((x_1,y_1)\)の値を得ることができます。

導出


陰関数

を満たす\((x,y)\)を求めます。ある点\((x_0, y_0)\)において\(f(x_0, y_0)=0\)を満たすことが分かっている場合、
陰関数の全微分より、\((x_0, y_0)\)の近傍で

が満たされています。\(f(x_0,y_0)=0\)であり、\(\Delta x\approx \Delta y \to 0\)ならば、\(O(\Delta^2)\)の項を無視して

が成立する変化のみが許されます。両辺を\(\Delta x\)で割って、\(\Delta x\to 0\)の極限をとれば、

を得ます。すなわち、式(9)を何らかの手段、例えば数値計算で解けば等高線を得ることができます。

発散への対処


さて数値計算によって微分方程式(9)を解いて等高線を描く場合、\(y\)が\(x\)の一価関数になることはほとんどありません。
そのため、式(9)の右辺の分母がゼロになることがあり、微分が発散するため計算を行うことができなくなることがあります。言い換えれば、発散しない範囲しか計算できなくなるということです。
これを回避するためには様々な方法があると思いますが、例えば変数変換を適宜行いながら積分を実行していくことでこの問題を回避することができます。

発散を回避するためのアイデアは以下の通りです。
\(y\)を\(x\)の関数と見たとき、微分が発散するということは、分母がゼロになることを意味しています。
しかし、その点において90°回転させて\(y\)を\(x\)の関数としてみれば微分はゼロに等しいため、問題なく計算を進めることができます。
図示すれば以下の通りになります。

左図は座標を変化せずにそのままの関数ですが、右図のようにある点を計算したいときにその点近くの傾きと等しいように回転した座標から見るようにします。こうすることで発散を回避することができます。

これから変数変換について考えていきます。新たな変数\(x’, y’\)を用いて\(x=x(x’, y’),~y=y(x’, y’)\)(またはその逆関数)と書けているならば、

と書いても良いです。すなわち式(9)を解く代わりに\(x’\)と\(y’\)に関する方程式

を解いても良いのです。具体的に変数変換


を考えます。これは座標を\(\theta\)回転させた座標系への変換です。もちろん、\(\theta=0\)のとき、\(x=x’, y=y’\)が成立し、恒等変換となります。微分方程式は

であり、\((x’, y’)\)の空間で解いていきます。

さて、回転角\(\theta\)は任意ですが、式(15)の右辺が発散しないような\(\theta\)を選ぶ必要があります。
また、数値計算を行う上では、傾きがほとんど0に近いほうが有利です。そのため、方針としては現在の傾きを計算し、その傾き分だけ逆方向に座標回転させることでこれを実現させます。
すなわち、

分だけ回転させます。現在の傾きと座標軸の傾きを一致させることで傾きがほとんどゼロの座標軸に回転させるのです。

プログラム

以下にプログラムを示します。
必要な情報は\(f(x,y)\)の\(x, y\)方向の偏微分\(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}\)の解析解です。関数名”fx”, ”fy”にそれぞれを入力してください。
計算は刻み幅制御の陽的ルンゲクッタ法を用いています。

具体例


具体的に関数\(f(x,y)=x^3-2xy+y^3\)の等高線を考えていきましょう。
関数の値がそれぞれ、\(0.01,1,2,5\)に等しいときの等高線を考えます。
具体的な初期値は、
\(\displaystyle
f(-2,1.3657)\approx 0.01,~f(-2,1.4646)\approx 1,~f(-2,1.5568)\approx 2,~f(-2,1.7977)\approx 5
\)
とします。
具体的に解いたのが左の図です。

”正解”としたgnuplotの等高線表示の結果と同じ結果が得られました。座標回転がうまく働いているようですね。微分の発散を抑えるだけでなく、数値計算量も抑えられているようです。

刻み幅\(h\)が正である場合と負である場合で進む方向が変わります。
また計算の途中で分岐がある場合、どうなるかは保証できません。

また、閉じた等高線が複数ある場合、その全てを計算することは本稿のプログラムだけではできません。初期値を求める別のプログラムを作成する必要があります。
更に、等高線がある値に等しい量を計算することはできません。初期値の値に等しい等高線を引きます。

追記


調べ物をしていたら

照井章, 佐々木建昭著『27.代数関数の陰関数描画について』
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0941-27.pdf
という素晴らしい資料がありました。

この文献では2つの方法が紹介されており、
1.\(f(x,y_0)=0\)を解いて\(x=x_1,x_2,\cdots,x_n\)を解く
2.陰関数定理から求める
という方法です。資料では、1.の方法を採用していましたが、本稿の方法は2.の方法でした。
1.はもれなく等高線を見つけることができますので、考えられる結果を得たい場合によい方法です。
一方、2.は全てのゼロ点を見つけることができませんが、計算が早くできます。別の方法などで注目したい等高線が分かっている場合に良い方法でしょう。

効率的に解くためには、場合によって使い分けることが一番です。

負×負=正はなぜ?

負は、正の世界と負の世界を入れ替える因子である。

Q氏とR博士の会話1


Q「なぜ負に負を掛けると正になるのですか?」

R「まず、正の世界に正の数を掛けたら正になる。」

Q「納得します。」

R「では、正の世界に負の数を掛けたらどうなるかな?」

Q「負でしょう。」

R「そうだ。正の世界から見たら、負の数を掛けることで世界が反転するのだ。」
R「では負の世界から、世界を反転させる負の数を掛けたらどうなるかな?」

Q「負の反転は正ですから、正です。」

R「その通り。なので負に負を掛けると正なのだ。」

Q氏とR博士の会話2


Q「正の数に、正を掛けたら正になるのですよね。」

R「そうだ。」

Q「負の数に、負を掛けたら正になるのですよね。」

R「そうだ。」

Q「自分の世界の数同士をかけたらどちらも正はおかしくないですか?」
Q「正×正=正に、負×負=負になるのが、道理でしょう?」

R「その通りだ。君の疑問は全く正しい。だが、人類の都合で決められてしまったのだ。」
R「正の数しか無かった頃、人類は”正×正=正”を基準にした。その後に”数”を負に拡張したのだ。」
R「負の数に広げても”正×正=正”を保ちつつ、”正×負=負”または”負×正=負”になるように拡張したいと考えた。」
R「そう制約したから、”負×負=正”となってしまったのだ。負は正の世界と負の世界を入れ替える因子である、と。」

Q「なるほど。もし仮に”負×負=負”を基準に選んでいたら、正の数が正の世界と負の世界を入れ替える因子となるわけですか。」
Q「”負×負=負”を基準とするならば、”正×負=正”または”負×正=正”であって、”正×正=負”になるように。こういった世界を作ることも可能だったんですね。」

R「その通りだ。どちらの考え方も正しく、優劣はないよ。」

Q「ならば、無理やり”正×正=正”と”負×負=負”を満たすように強引に数学を組み立てたら、”正×負”はどうなるのでしょう?」

R「未定だ。定義できない。正×負の単位を数えるための、新たな数の拡張が必要になるだろうね。」

秋間補間

秋間補間と呼ばれる補間方法があります。

こんな補間結果となります(akimaとmodified akimaの線が秋間補間の結果です)。

秋間補間の特徴は、

  • 与えられた点を必ず通る
  • オーバーシュートが発生しにくい
  • 補間したい点から遠い補間点の影響を受けない
  • 導関数の連続性は保証されない
  • 元のデータ点を大きく逸脱するような補間は苦手(線形的な補間に近い)
  • 極値におけるデータ点など、元の関数の特徴的な値が分かっているときに良い補間結果を与える

ということです。オーバーシュートとは、元のデータが急峻に変化した場合、その前後の補間結果が振動してしまい大きく外れてしまうことを意味します([5]にオーバーシュートの様子を見ることができます)。
ただし、3次スプライン補間のような微分の連続性が保証されません。

それでは計算方法を追っていきます。

導出


秋間補間は1970年に秋間浩によって提案された補間方法です[1]。地図の等高線を表示するために生成されたようです。その後修正された秋間補間が紹介されました[2]。

具体的な計算方法を見ていきます。
補間は補間したい点を中心に前後のデータ点を用いるので、両端では特別な扱いが必要です。それを区別するために領域を3つに分けて考えます。具体的には下の図のような感じです。

領域IIにおける補間


まずは領域IIを考えます。簡単に考えるために、以下の図のように補間に用いるための6点\((x_i, y_i), i=0,5\)があった時、中間の領域である\(x_2\le x\le x_3\)の中の値を補間して求めたいと考えます。

補間は3次関数で行います。すなわち\(x_2\le x\le x_3\)において

の形で補間を行います。係数\(a_0, a_1, a_2, a_3\)は補間に用いる6点のデータ点から計算できて、

の通り求められます。ここで、\(q_1, q_2, m_1,\cdots, m_5\)は


の通り求められます。\(w_1, w_2\)は使用する秋間補間の種類によって変わります。通常の秋間補間では、

の通り計算し[1,3,4]、修正秋間補間においては

で計算します[2,3,4]。

領域I, IIIにおける補間


領域I, IIIについては以下の関係を満たすように、本来は存在しない補間用のデータ点を予測します。
領域I:最も左端のデータ点を\(x_3\)、その右隣のデータ点を\(x_4, x_5\)と置きます。すると、\(x_3\)より左のデータ点\(x_2, x_1\)を関係式(7)より予測します。
領域III:最も右端のデータ点を\(x_3\)、その左隣のデータ点を\(x_2, x_1\)と置きます。すると、\(x_3\)より右のデータ点\(x_4, x_5\)を関係式(7)より予測します。

ここで、\(m_{ij}\)は

を意味します。すると、秋間補間で用いる(本来は存在しない)端のデータ点は、

領域I)

領域III)

と求められます。

プログラム例


上の式をプログラムすれば、このようになるかと思います。
注意してプログラムを作らないとゼロ割が発生します。\(x_i\ne x_j\)を仮定していれば、式(5),(6)の\(w_1, w_2\)を計算する際にこれが発生します。なので、この回避を取り入れてプログラムを作成しています。

計算例


参考文献[4]に秋間補間と修正秋間補間の振る舞いの例があります。それを本稿のプログラムを用いて計算すると同じ振る舞いを示すことが分かります。

それでは、3次スプライン補間と秋間補間、修正された秋間補間を比べて、秋間補間にどのような特徴があるかみましょう。

オーバーシュートの抑圧


本当の答えを定義することは難しいですが、変化する点(\(x=2.8\sim 3.2\))近傍で直線とすれば傾きは5になります。一階導関数はどの補間方法も連続(右から極限と左から極限が一致の意味で)ですが、秋間補間は二階導関数は連続ではありません。3次スプラインは二階導関数まで連続であるという仮定の下での補間なので、一応連続的な変化をしています。

秋間補間では3次スプライン補間よりはオーバーシュートが抑えられているものの、近傍 \((x=2.5, 3.5)\)において外れた値となっています。
修正された秋間補間では、その点が改善されオーバーシュートが発生していないことが分かります。

振動関数の補間

振動関数を補間する例として、元の関数が\(y=\sin(x)\)だったとします。

sin関数

この場合は元の関数が十分滑らかであり、高階微分も連続であることが分かっています。与えられたデータ点も誤差がないので、3次スプラインが有効でしょう。
秋間補間では本来極値をとるはずの\(x=n\pi/2,~(nは整数)\)あたりで値があまり大きくならないことから、かなり局所的な補間である特徴が分かります。すなわち、大きくオーバーシュートしない反面、予測は得意ではないということです。これを回避するためには、元の関数の極値の値など特徴的なデータが必要であることが分かります。

二次元平面の補間は本稿には載せていませんが、[4,6]にどのような違いが表れるか図があるので、興味がある方は参照してみてください。
確かに、本来の目的である地図の等高線作成には大きな効果が期待できそうなことが分かります。おそらく、取得するデータ点も標高など、地形の特徴的な点におけるデータを集めることが主となるので、秋間補間との相性が良いのでしょう。

参考文献


[1]Akima H. A New Method of Interpolation and Smooth Curve Fitting Based on Local
Procedures. Journal of the ACM, 17(4), 589 ― 602. (1970).
doi:10.1145/321607.321609

[2]Akima H. A method of univariate interpolation that has the accuracy of a third-degree
polynomial. ACM Transactions on Mathematical Software, 17(3), 341 ― 366. (1991),
doi:10.1145/114697.116810

[3]
makima 修正Akima区分的3次エルミート内挿 -MathWorks
interp1 1次データ内挿(テーブルルックアップ) -MathWorks

[4]Cleve Moler, Makima Piecewise Cubic Interpolation

[5]スプラインによる補間

[6]格子点に内挿して等値線をひく

複素関数の微分と積分

複素関数論のメモです。

  1. 微分
  2. 積分
  3. 参考文献

微分


コーシー・リーマンの関係式(まとめ)


複素関数がコーシー・リーマンの関係式を満たすと、その複素関数は通常の変数通りに微分することが出来ます。複素数\(z=x+iy\)を引数に持つ複素関数\(f(z)=f(x+iy)\)が、実数値関数\(u(x,y), v(x,y)\)を用いて

と書けているとき、コーシー・リーマンの関係式

を満たすならば、\(f(z)\)は正則関数と呼ばれます。ここで、
\(f_x(x,y)=\frac{\partial f}{\partial x}, f_y(x,y)=\frac{\partial f}{\partial y}\)を意味します。式(2)を満たす場合、複素関数\(f(z)\)の複素数\(z\)による微分は

と書くことができ、式(3), (4)は同じ値となります。

コーシー・リーマンの関係式(導出)


複素関数の引数の数は1つであり、それは複素数です。すなわち1変数関数です。故に、複素関数の微分を定義しようとすれば、
\begin{align}
\frac{df(z)}{dz}
\end{align}
という量は1つに定まるはずです。しかし実数の空間で生活している我々には、1つの複素数を実部と虚部に分けて考えてしまう癖のようなものがあります。そのせいで、複素関数は平面に記述される関数であり、実部方向と虚部方向をそれぞれ独立に考えることが出来る2変数関数だと勘違いしてしまいますが、そうではありません。

すなわち、2変数の実数値関数と1変数の複素関数は全く異なる関数であり、同じ様に扱ってはいけません。しかし、我々は1次元上に複素数を表現する方法を持っておらず、本来、1変数の複素関数を2変数の実数値関数のように表現しなければイメージすることが出来ません。この場合、複素関数の表現に何らかの制約がかかることがイメージできるでしょう。その制約が正則である、という条件であり、コーシー・リーマンの関係式なのです。

では実際に本来1変数関数の複素関数を実軸と虚軸という、あたかも2変数関数のように拡張して考えて、微分を考えてみましょう。

前提として、複素関数の微分
\begin{align}
\frac{df(z)}{dz}
\end{align}
は一意に定まるとします。この下で2変数関数として\(f(z)=f(x,y)\)を拡張した場合、微分したい点への近づき方が二通りあることに気が付きます。それは、\(x\)方向と\(y\)方向ということです。2変数関数であれば、偏微分の事です。どんな方向から近づいても微分値は一意に定まると考えているので、両結果は等しくなければなりません。
微分は実数値関数と同様に差分の極限として定義して、
\begin{align}
\frac{df(z)}{dz}\equiv \lim_{\Delta z\to 0}\frac{f(z+\Delta z)-f(z)}{\Delta z}
\end{align}
と書けるとします。ここで、\(\Delta z = \Delta x+i\Delta y\)を意味しており、\(\Delta x, \Delta y\)は実数です。

点\(z=(x, y)\)に実数軸に沿って近づく場合、\(\Delta y=0\)を考えればよいので、微分は差分の極限として
\begin{align}
\frac{df(z)}{dz}=\lim_{\Delta x\to 0}\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}
\end{align}
と書かれます。ここで\(h\)は実数です。また、虚数軸に沿って近づく場合\(\Delta x=0\)を考えればよいので、
\begin{align}
\frac{df(z)}{dz}=\lim_{\Delta y\to 0}\frac{f(x,y+\Delta y)-f(x,y)}{i\Delta y}
\end{align}
と書かれるはずです。複素関数の微分は一意に決まることから、二つの近づき方を考えたとしても一致していなければなりません。複素関数が\(f(z)=u(x,y)+iv(x,y)\)と実部と虚部があらわに分離されて書かれていたとするならば、
\begin{align}
&\lim_{\Delta x\to 0}\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}\
&=\lim_{\Delta x\to 0}\frac{[u(x+\Delta,y)+iv(x+\Delta,y)]-[u(x,y)+iv(x,y)]}{\Delta x} \
&=\lim_{\Delta x\to 0}\frac{[u(x+\Delta,y)-u(x,y)]+i[v(x+\Delta,y)-v(x,y)]}{\Delta x} \
&=\frac{\partial u}{\partial x}+i\frac{\partial v}{\partial x} \
\end{align}
と書けます。また、虚軸方向は
\begin{align}
&\lim_{\Delta y\to 0}\frac{f(x,y+\Delta y)-f(x,y)}{i\Delta y}\
&=\lim_{\Delta y\to 0}\frac{[u(x,y+\Delta y)+iv(x,y+\Delta y)]-[u(x,y)+iv(x,y)]}{i\Delta y} \
&=\lim_{\Delta y\to 0}\frac{-i[u(x,y+\Delta y)-u(x,y)]+[v(x,y+\Delta y)-v(x,y)]}{\Delta y} \
&=-i\frac{\partial u}{\partial y}+\frac{\partial v}{\partial y}
\end{align}
と書くことが出来ます。
両者はもともと同じものであるとしてきましたから、実部と虚部で比較します。
すると、

という関係式が導けるという訳です。これがコーシー・リーマンの関係式と呼ばれるもので、1変数関数の複素関数を実部、虚部という2変数関数に無理矢理拡張した場合に現れる条件です。
複素関数の微分が定義できるとき、複素関数の実部虚部は、この関係式を満たしていなければなりません。ある範囲においてすべての点で複素関数の微分が定義できることを、正則であると言ったり、滑らかである、解析的であると言います。
”複素関数の滑らか”は通常の”滑らか”と意味が異なることに注意しましょう。通常は不連続点が無いことを指しますが、複素関数の場合はそれ以上の制約(コーシー・リーマンの関係式)を含んでいます。

また、ここでは述べませんが、複素関数は1度微分できると何回でも微分することが出来ます。
さらに微分が定義できる正則関数であれば定義域を超えた領域へ関数が拡張できる場合があります。これが解析接続と呼ばれるテクニックです。コーシー・リーマンの関係式を制約として今まで捉えてきたのですが、この制約があるからこそまだ見ぬ定義域外をこの関係式を手掛かりに関数を拡張できるのです。

例として、究極的に滑らかそうに見えるけど微分できない関数を上げましょう。
\begin{align}
f(z)=z=x+iy
\end{align}
はコーシー・リーマンの関係式を満たすため、複素関数として微分が定義でき、正則です。しかし、その複素共役を取った関数
\begin{align}
f(z)=z^*=x-iy
\end{align}
はコーシー・リーマンの関係式を満たさないため、複素関数として微分が定義できず、正則ではありません。近づき方によって微分の値が変わってしまうのです。

通常の微分が可能という意味を明確にするため、例を挙げます。\(f(z)\)が\(z\)の多項式で書かれている場合、微分は

として求める事が出来ます。通常と同じですね。これは、複素数の成分を持たない場合には実軸上で考えたものと同等になるので、実数軸上の微分を自然に拡張したようになっています。

積分


特異点周りの積分


端的にまとめます。
点\(a\)を囲うように複素平面上の閉じた積分経路の方向を\(C\)で指定するとします。\(z\)の多項式で表現される関数\(f(z)=(z-a)^n\)の積分は、

となります。点\(a\)の周りを反時計回りに1周回る場合、\(C\)はパラメータ\(t\)を用いて\(z(t)\)と書くと、
\begin{align}
z(t)=a+re^{it},~(0\le t \le 2\pi)
\end{align}
と書けます。\(r\)は積分経路の半径を意味します。式(6)の計算をすれば、

となります。\(n=-1\)の時、右辺の積分は\(2\pi\)に等しい事は明らかです。それ以外の\(n\ne -1\)ならば、

となります。よって、

が導かれます。

留数定理


一つ前の、ある点周りの積分をもう少し一般的に考えましょう。\(z=a\)を中心とした閉曲線に沿った複素関数\(f(z)\)の複素積分\(\oint_C f(z)dz\)は

と書けます。ここで\(\text{Res}[f,a]\)は留数と呼ばれる量です。これから留数について述べていきます。

\(\text{Res}[f,a]\)は、\(f(z)\)が点\(a\)周りで

とローラン展開出来るとします。この時、\(a_{-1}\)を留数と呼びます。すなわち、

を意味します。

では、留数を実際に求めていきましょう。\(f(z)\)が点\(a\)周りでローラン展開可能な場合、\(n\)の下限が重要になります。下限を極の位数\(k\)と表現します。\(k\)は、

を満たすときの値として定義できます。

もし、関数\(f(z)\)が\(z=a\)で\(k\)位の極であれば、\(z=a\)周りで

と展開できます。
留数は\(a_{-1}\)なので、\(a_{-1}\)について式変形を行っていけば、

と計算できます。

参考文献


原惟行、松永秀章著、『複素解析入門』(第3刷2010, 初版2007)

Hyper-dual numbersによる二階偏微分の計算

Hyper-dual numbersと呼ばれる、実数を拡張した考えを導入すると導関数が計算できます。

あらかじめ、引数がHyper-dual numbersである時の数々の関数の定義を実装しておけば、その関数の組み合わせで作られる関数の一階導関数、二階導関数のほぼ厳密な答えを得ることが出来ます。
この考えに従う関数の微分方法は、Forward型の自動微分と呼ばれます。

  1. Dual number
    1. Dual numberと導関数
  2. Hyper-dual numbers
    1. Hyper-dual numbersの演算
    2. Hyper-dual numbersのプログラム
    3. ヘッセ行列
  3. 参考文献

Dual number


Dual numberという考えがあります[2,3,5]。

これは実数を拡張する、という考えで複素数に似た考え方です。

良く知られている実数の拡張方法の一つは、複素数です。
通常の実数に\(i^2=-1\)を満たす\(i\)という数を付加するのが複素空間です。

拡張の方法は何も\(i\)だけではありません。
例えば、\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす\(\epsilon\)という数を付加することもできます。

この\(\epsilon\)は複素数ではありません
複素平面上でこの性質を満たす数は無いことからも、この\(\epsilon\)を追加するということは新しい方向への数の拡張です。
実際、複素平面上で\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす数があるのかを探しても
\(
\begin{align}
(a+ib)^2=a^2-b^2+i2ab=0
\end{align}
\)

であるので、これを満たすのは\(a=b=0\)しか存在せず、複素平面上の数ではないことが分かります。

そんな新しい方向\(\epsilon\)の平面で定義された数をDual number(二重数、または双対数)と呼びます。

二重数\(a\)は実数部と非実数部から構成されており、
\(
a=a_0+a_1\epsilon
\)

と書くことが出来ます。ここで、\(a_0, a_1\)は実数であり、 \(\epsilon\)は、虚数単位\(i\)に倣って二重数単位とでも名付けておきましょう。

Dual numberと導関数


二重数は関数の導関数と大きな関係がある事を示しましょう。
テーラー展開を行います。
関数\(f\)の\(x\)周りの展開は刻み幅\(\Delta\)を用いて
\(\displaystyle
f(x+\Delta)=f(x)+f'(x) \Delta+\frac{1}{2!}f^{\prime\prime}(x)\Delta^2+\cdots
\)

と記述することが出来ます。もし仮に\(\Delta\)が純非実数だとしましょう。すなわち、\(\Delta=h\epsilon\)とします。ここで、\(h\)は実数です。テーラー展開の式に代入すれば、
\(
\begin{eqnarray}
f(x+h\epsilon)&=& f(x)+f'(x) h\epsilon+\frac{1}{2!}f^{\prime\prime}(x){h}^2\epsilon^2+\cdots \\
&=& f(x)+f'(x) h\epsilon
\end{eqnarray}
\)

となるわけです。注記しますが、上の式は\(h\)の1次で打ち切っているのではなく、厳密にイコールが成り立っています。
この式が言っているのは、関数\(f(a), (a\text{は二重数})\)を計算し、その非実数部を\(h\)で割ると関数の導関数になっているということです。

二重数の非実数部を取り出す関数\(\text{Dual}\)を定義します。
すると、導関数は
\(\displaystyle
f'(x)=\frac{1}{h}\text{Dual}(f(x+h\epsilon))
\)

として得られます。
言葉で書けば、二重数空間に拡張した\(f(x+h\epsilon)\)を計算すると、その非実数部に導関数が現れる、ということです。

二重数のプログラムはJeffrey Fikeさんによる[4]にありますので、そちらをご参考ください。

Hyper-dual numbers


さて、ここまでで二重数の概念を簡単に説明し、新しい実数の拡張を行いました。
Dual Numberのままでは高階導関数は得られません。なぜなら、導関数の二次以降は二重数の性質\(\epsilon^2=0\)によってゼロになるからです。

そこで、若干工夫します。
Dual Numberを二種類用意することで二階微分を得ることが出来ます[2,3]。
すなわち、
\(\begin{gather}
a=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
\epsilon_1^2=\epsilon_2^2=(\epsilon_1\epsilon_2)^2=0,~\epsilon_1\ne 0,~~\epsilon_2\ne 0,~~\epsilon_1\epsilon_2 \ne 0
\end{gather}
\)

と実数を拡張します。この様に拡張した数\(a\)をHyper-dual numbersと呼びます[2,3]。

Hyper-dual numbersの演算


Hyper-dual numbersである
\(
\begin{align}
a&=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
b&=b_0+b_1\epsilon_1+b_2\epsilon_2+b_3\epsilon_1\epsilon_2
\end{align}
\)

を用意します。和、積、商はそれぞれ
\(
\begin{align}
a+b&=(a_0+b_0)+(a_1+b_1)\epsilon_1+(a_2+b_2)\epsilon_2+(a_3+b_3)\epsilon_1\epsilon_2\\
ab&=a_0b_0+(a_0b_1+a_1b_0)\epsilon_1+(a_0b_2+a_2b_0)\epsilon_2+(a_0b_3+a_1b_2+a_2b_1+a_3b_0)\epsilon_1\epsilon_2\\
\frac{1}{a}&=\frac{1}{a_0}-\frac{a_1}{a_0^2}\epsilon_1-\frac{a_2}{a_0^2}\epsilon_2+\left(-\frac{a_3}{a_0^2}+\frac{2a_1a_2}{a_0}\right)\epsilon_1\epsilon_2
\end{align}
\)
と定義されます[2]。任意の関数は
\(
f(x)=f(x_0)+x_1f'(x_0)\epsilon_1+x_2f'(x_0)\epsilon_2
+\left(x_3f'(x_0)+x_1x_2f^{\prime\prime}(x_0)\right)\epsilon_1\epsilon_2
\)

として計算することが出来ます。なので、この結果から\(\epsilon_1\epsilon_2\)の係数として関数の二階微分が得られます。

Hyper-dual numbersのプログラム


実際にプログラムを組みましょう。
Hyper-dual numbersの型を持つ変数はFortran90では定義できません。
なので、構造体を利用して自分で型と、演算を定義します。

Fortran90のプログラムは以下の通りになるかと思います。
基本的な四則演算、基本的な初等関数のHyper-dual numbersの演算をモジュールとして書いています。

上のモジュールと下のメインプログラムを一緒にコンパイルすることにより、関数の二階微分が得られます。

例として
2変数関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)

の偏微分
\(\begin{align}
\frac{\partial f}{\partial x},~~\frac{\partial f}{\partial y},~~\frac{\partial^2 f}{\partial x\partial y}
\end{align}
\)

を得ることを考えます。
\(
f(x+1\epsilon_1,y+1\epsilon_2)
\)

を計算すると
\(
f(x+1\epsilon_1,y+1\epsilon_2)=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2
\)

のように実数係数\(a_0, a_1, a_2, a_3\)が得られます。
すると、
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial y}=a_2 \\
\frac{\partial^2 f}{\partial x\partial y} =a_3
\end{align}
\)

として偏微分が得られます。
ちなみに、二階微分が欲しい場合は
\(
f(x+1\epsilon_1+1\epsilon_2, y)
\)

を考えると
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial x}=a_2 \\
\frac{\partial^2 f}{\partial x^2} =a_3
\end{align}
\)

として得られます。\(a_1,~a_2\)はどちらを採用しても構いません。

プログラムでは変数の型Hyperdualを持つ入力変数を\(\text{xH,yH}\),出力を\(\text{rH}\)と置いています。

program main
  use Hyperdualmod
  implicit none
  type(Hyperdual)::xH,yH,rH
 
  xH%x0 = 0.3d0 ! real part
  xH%x1 = 1d0   ! unreal part \epsilon_1
  xH%x2 = 0d0   ! unreal part \epsilon_2
  xH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  yH%x0 = 0.4d0 ! real part
  yH%x1 = 0d0   ! unreal part \epsilon_1
  yH%x2 = 1d0   ! unreal part \epsilon_2
  yH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  write(6,'(4f23.16)')xH%x0,xH%x1,xH%x2,xH%x3
  write(6,'(4f23.16)')yH%x0,yH%x1,yH%x2,yH%x3  
 
  rH = log(xH*yH**2)*exp(xH)/sqrt(sin(xH)**3+cos(yH)**3)
  !rH = asin(2d0*xH)*acos(yH)/atan(xH*yH)
  !rH = xH**yH
 
  write(6,'(4f23.16)')rH%x0,rH%x1,rH%x2,rH%x3  

  stop
end program main

ヘッセ行列


二階偏微分の計算が出来たので、ヘッセ行列が簡単に計算できます。
ルーチンを作れば、以下の通りになるかと思います。
下のプログラムは3変数関数
\(\displaystyle
f(x,y,z)=\exp(xy)\tan(z)
\)

の\(x=-2,~y=3,~z=1\)におけるヘッセ行列を計算します。

program main
  use Hyperdualmod
  implicit none
  integer::N
  double precision::x,y,z,f
  double precision,allocatable::nabla(:),Hesse(:,:)
  double precision,allocatable::w(:)
  external::func

  N=3

  allocate(nabla(1:N),Hesse(1:N,1:N))
  nabla=0d0
  Hesse=0d0
  allocate(w(1:N))
  w=0d0

  x = -2d0
  y =  3d0
  z =  1d0
  w(1) = x
  w(2) = y
  w(3) = z
  call Hessian(N,w,func,nabla,Hesse)

  f=exp(x*y)*tan(z)
  write(6,'(1e24.16)')nabla(1)
  write(6,'(1e24.16)')y*f
  write(6,'(1e24.16)')nabla(2)
  write(6,'(1e24.16)')x*f
  write(6,'(1e24.16)')nabla(3)
  write(6,'(1e24.16)')exp(x*y)/(cos(z)**2)
  write(6,'(3e24.16)')Hesse(1,1:3)
  write(6,'(3e24.16)')f*y**2, f*(1+x*y), f*y/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(2,1:3)
  write(6,'(3e24.16)')f*(1+x*y), x**2*f, f*x/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(3,1:3)
  write(6,'(3e24.16)')f*y/(cos(z)*sin(z)), f*x/(cos(z)*sin(z)), f*2d0/(cos(z)**2)

  stop
end program main

subroutine func(N,x,f)
  use Hyperdualmod
  implicit none
  integer::N
  type(Hyperdual),intent(in)::x(1:N)
  type(Hyperdual),intent(out)::f
 
  f = exp(x(1)*x(2))*tan(x(3))
 
  return
end subroutine func

subroutine Hessian(N,x,func,nabla,Hesse)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::x(1:N)
  double precision,intent(out)::nabla(1:N)
  double precision,intent(out)::Hesse(1:N,1:N)
  external::func

  integer::i,j,k
  type(Hyperdual),allocatable::xH(:)
  type(Hyperdual)::rH
 
  allocate(xH(1:N))
  do i=1,N
     xH(i)%x0 = 0d0
     xH(i)%x1 = 0d0
     xH(i)%x2 = 0d0
     xH(i)%x3 = 0d0
  enddo
 
  do i=1,N
     xH(i)%x0 = x(i)
     xH(i)%x3 = 0d0
  enddo

  do i=1,N
     do k=1,N
        xH(k)%x1 = 0d0
        xH(k)%x2 = 0d0
     enddo
     xH(i)%x1 = 1d0
     xH(i)%x2 = 1d0
     
     call func(N,xH,rH)
     nabla(i) = rH%x1
     Hesse(i,i) = rH%x3
  enddo

  do i=1,N
     do j=i+1,N
        do k=1,N
           xH(k)%x1 = 0d0
           xH(k)%x2 = 0d0
        enddo
        xH(i)%x1 = 1d0
        xH(i)%x2 = 0d0
        xH(j)%x1 = 0d0
        xH(j)%x2 = 1d0
        call func(N,xH,rH)
        Hesse(i,j) = rH%x3
        Hesse(j,i) = Hesse(i,j)
     enddo
  enddo  
 
  return
end subroutine Hessian

Hyperdual.f90にモジュールを入れ、メインプログラムをmain.f90に入れたとすると、以下の結果を得ます。

> gfortran Hyperdual.f90 main.f90  
> ./a.out
  0.1158128336233602E-01
  0.1158128336233602E-01
 -0.7720855574890680E-02
 -0.7720855574890680E-02
  0.8491012233306163E-02
  0.8491012233306162E-02
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458059E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458058E-01

実行結果の奇数行目はHyper-dual Numberによる計算結果、偶数行目は解析解を表します。
また、6行目までは一階微分、7行目以降はヘッセ行列を表します。

参考文献


[1]関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)

の偏微分\(\frac{\partial^2 f}{\partial x\partial y}\)
の計算
https://www.wolframalpha.com/input/?i=D%5Bln(x*y%5E2)*e%5E(x)%2Fsqrt(sin%5E3(x)%2Bcos%5E3(y)),+y,x%5D

[2]Jeffrey A. Fike and Juan J. Alonso,~”The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations”, 49th AIAA Aerospace Sciences Meeting(2011)
http://adl.stanford.edu/hyperdual/fike_aiaa-2011-886_slides.pdf,
J. A. Fike and J. J. Alonso. The Development of Hyper-Dual Numbers for Exact Second Derivative Calculations. AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting, January 4-7, 2011.http://adl.stanford.edu/hyperdual/Fike_AIAA-2011-886.pdf

[3]JeffreyA.Fike,~”Derivative Calculations Using Hyper-Dual Numbers”, Sandia National Laboratories (2016)
https://www.osti.gov/servlets/purl/1368722

[4]Jeffrey Fike, Aerospace Design Lab, http://adl.stanford.edu/hyperdual/

[5]松本佳彦, 新しい数をつくる, (2018) http://ymatz.net/assets/docs/20180629-jtpa-slide-mod

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

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

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

まとめ

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}
\)

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

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

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

実際に、LAPACKのLU分解を使用したプログラムの例を載せておきます。

高次のニュートン法?


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