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

位相速度、群速度とは!?調べてみました!!!!

近頃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

exp(-ikx)/x の無限区間に渡る積分を数値的に計算したい。

問題:以下の積分を数値的に計算せよ。
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)


1/xのフーリエ変換である
\( \displaystyle
\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)

数値的に計算します。

この積分値は一意には決まりません。
複素関数論によると、この積分値は\(x=0\)周りに存在する特異点を上周りに迂回して積分するか、下周りに迂回して積分するか、で答えが異なります。

特異点を上周りに回る場合、答えは
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
-2\pi i~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
であり、特異点を下周りに回る場合、
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x-i\varepsilon}dx =
\left\{
\begin{aligned}
0~~~(k\gt 0)\\
2\pi i~~~(k \lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
となります。
私が疑問に思うのは、積分すると本当にステップ状の関数が現れるのか?ということです。にわかには信じられません。

ここでは、上を回った時の計算にのみ注目します。
答えに出てくる係数を減らしたいため、被積分関数
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx
\)

を考えます。

数値計算方法


さて、問題を整理しましょう。今取り組む問題は、

問題:以下の積分を数値的に計算せよ。
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)

です。読み取れる情報は

・始点は、実軸上のマイナス無限大
・終点は、実軸上のプラス無限大
・\(x=0\)に特異点が存在する

という情報です。これらは理由なく勝手に変えてはいけません。
被積分関数が収束するならば、xがどんなに大きくても無限遠まで振動していきますが、フレネル積分のように積分結果は収束することを期待しましょう。

さて、\(x=0\)に特異点があるために、実軸上のマイナス無限からプラス無限まで積分していこうとすると\(x=0\)を通過しようとする時に無限大になるため、オーバーフローします。
ということで、積分経路を特異点周りで迂回します。

始点と終点は決められていますが、積分経路の指定はありませんので何も問題ありません。
積分経路は以下の図のように取ります。

このようにとると、積分経路上に特異点はないですので計算することが出来ます。
始点と終点は十分大きい値、ここでは(±2000,0)という値にします。

横軸をkとして、積分を計算すると以下のようになります。

複素関数論から、特異点を複素平面上方を迂回した場合に得られる解は
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} -\frac{1}{2\pi i}\int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
1~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
です。数値計算によって得られた解は、\(k\lt 0\)で値0, \(k\gt 0\)で値1であるので、数値的にも正しそうです。
しかし、厳密にそれぞれの値になっているわけではないことに気が付きます。なぜならば、\(k=0\)周りでは値は収束していないように見えます。
この原因は始点と終点としてとった、大きい値(±2000,0)が十分ではない事を意味しています。kが小さいと、被積分関数の\(\exp(-ikx)\)による振動はゆっくりになるため、収束が遅くなっている、と考えられます。

ジョルダンの補題の利用


この振動を無くすためには被積分関数が漸近領域で収束するように積分経路を変えればいいのです。

始点と終点は変わりません。図は\(k\)が負の時の経路です。

経路の取り方は無数にあり、上の経路は一つ前の計算で選んだ経路を出来る限り変更しないように経路をとっているだけにすぎません。

ジョルダンの補題を利用すれば、被積分関数が収束する領域であれば積分の寄与はゼロなので、数値的に積分する場合の視点、終点を複素平面上へ移すことが出来ます。
この考えの下実際に計算すると積分結果はこうなります。

綺麗なエッジが作られました。これならば解析解と数値解が一致した、と言えるでしょう。複素平面上の数値計算は一筋縄ではいかないことが実感できました。

数値計算プログラム


プログラムは以下のものを用いています。

数値積分はquadpackのガウス=ルジャンドル求積法の自動積分コードを複素関数用に変更して用いています。
実際に使用しているquadpackコードは以下のもので、
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag.f90
これと、下のはmainプログラムです。

program main
  implicit none
  integer::ier,Nt,i,j,sn
  double precision::eps
  complex(kind(0d0))::a,b,s,s0
  complex(kind(0d0))::xa(1:5),xb(1:5)
  complex(kind(0d0)),external::g
     
  eps=1d-10
  ier=1

  Nt=1000
  do i=1,Nt
     t=i*(1d0-(-1d0))/dble(Nt)-1d0
     s=dcmplx(0d0,0d0)
     
     sn=1
     if(t.gt.0d0)sn=-1
     if(t.eq.0)sn=0
     
     xa(1)=dcmplx(-2000d0,sn*2000d0)
     xb(1)=dcmplx(-1d0,0d0)

     xa(2)=dcmplx(-1d0,0d0)
     xb(2)=dcmplx(-1d0,1d0)

     xa(3)=dcmplx(-1d0,1d0)
     xb(3)=dcmplx( 1d0,1d0)

     xa(4)=dcmplx( 1d0,1d0)
     xb(4)=dcmplx( 1d0,0d0)

     xa(5)=dcmplx( 1d0,0d0)
     xb(5)=dcmplx( 2000d0,sn*2000d0)  

     do j=1,5
        s0=dcmplx(0d0,0d0)
        call cqag_sk(g,xa(j),xb(j),eps,s0,3,ier)
        s=s+s0
     enddo
     write(10,*)t,dble(s),dimag(s)
  enddo  
     
  stop
end program main

function g(z)
  implicit none
  complex(kind(0d0))::g
  complex(kind(0d0)),intent(in)::z

  double precision::pi=dacos(-1d0)
 
  g=exp(-dcmplx(0d0,1d0)*z*t)/(z)
  g=g*(-1d0/(2d0*pi*dcmplx(0d0,1d0)))

  return
end function g

オリジナルのquadpackコードは
http://www.netlib.org/quadpack/(パブリックドメイン)です。

1/(x^n)のフーリエ変換

その他フーリエ変換
ヘヴィサイド関数のフーリエ変換
1のフーリエ変換
\(x\)のフーリエ変換
\(1/x\)のフーリエ変換
\(1/(x^n)\)のフーリエ変換 ← 今ここ

\(\frac{1}{x}\)のフーリエ変換から\(\frac{1}{x^n}\)のフーリエ変換を導けます。
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx =-\pi i ~\text{sgn}(k)
\end{align}
\)

(導出は\(\frac{1}{x}\)のフーリエ変換へ)

部分積分を用いると\(\frac{1}{x^2}\)のフーリエ変換
\(
\begin{align}
&\int_{-\infty}^{\infty}\frac{1}{x}\left(\frac{1}{-ik}e^{-ikx}\right)’dx \\
&=\left[\frac{1}{x}\cdot\frac{1}{-ik}e^{-ikx}\right]_{-\infty}^{\infty} -\frac{1}{ik}\int_{-\infty}^{\infty}\frac{1}{x^2}e^{-ikx} dx
\end{align}
\)
第一項はゼロになります。第二項と残りの項と\(1/x\)のフーリエ変換の右辺が等しいので、
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x^2}e^{-ikx} dx =-\pi k ~\text{sgn}(k)
\end{align}
\)
が得られます。
これを繰り返すことで\(\frac{1}{x^n}\)のフーリエ変換
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x^n}e^{-ikx}dx &=-\pi i \frac{(-ik)^{n-1}}{(n-1)!}~\text{sgn}(k)
\end{align}
\)

が導かれます。

1/xのフーリエ変換

\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx =-\pi i ~\text{sgn}(k)
\end{align}
※ただし、\(x=0\)での値は定義しない。また、\(\text{sgn}(x)\)は符号関数

数値計算による上式の確かめは1/xのフーリエ変換をご覧ください。

その他フーリエ変換
ヘヴィサイド関数のフーリエ変換
1のフーリエ変換
\(x\)のフーリエ変換
\(1/x\)のフーリエ変換 ← 今ここ
\(1/(x^n)\)のフーリエ変換


実軸上の積分は複素平面の上から近づいた値と下から近づいた値の平均をとって定めるとします。すなわち、
\(
\begin{align}
&\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx \nonumber \\
&=\lim_{\varepsilon\to +0} \frac{1}{2}\left[\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx +\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}e^{-ikx}dx \right]
\end{align}
\)
で表されると仮定します。
そして、第1項、第2項をそれぞれ計算していきます。

\(
\displaystyle \lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x\pm i\varepsilon}e^{-ikx}dx
\)

を計算するにあたり、複素x平面上の漸近で収束する半径が\(k\)の符号によって異なるので場合分けをしてます。

\(k\gt 0\)の時

まず、\(\varepsilon\)の符号が正である場合を考えます。
この実軸上の積分は実軸への近づき方を指定して
\(
\begin{align}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx =
\left[\int_{-\infty}^{-\varepsilon’}+\int_{\Gamma_+}+\int_{\varepsilon’}^{\infty}\right]\frac{e^{ikx}}{x}dx
\end{align}
\)

です。ここで、\(\int_{\Gamma_+}\)は特異点周りを複素平面上半面を通る経路です。
この積分を計算するために閉路積分を考えます。すなわち、以下の図の左側のように回った積分を考えます。

すると、この閉経路に対して、
\(
\begin{align}
\oint \frac{e^{-ikx}}{x}dx&=
\lim_{R\to\infty}\left[\int_{-R}^{-\varepsilon’}+\int_{\Gamma_+}+\int_{\varepsilon’}^{R}+\int_{\Gamma’}\right]
\frac{e^{ikx}}{x}dx\\
&=\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx+\lim_{R\to\infty}\int_{\Gamma’}\frac{e^{-ikx}}{x}dx
\end{align}
\)
が成立します。今、閉経路内に特異点が一つ存在し、負の方向に回っているので、左辺\(\oint \frac{e^{-ikx}}{x}dx=-2\pi i\)です。また、大きく外側を回る経路\(\Gamma’\)の積分はゼロです。よって、
\(
\begin{align}
\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx =-2\pi i.
\end{align}
\)

が示されます。

また、特異点を下に回る積分(図の左の経路)では、閉経路内に特異点は存在しないので、\(\oint \frac{e^{-ikx}}{x}dx=0\)です。なので、
\(
\begin{align}
\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}e^{-ikx}dx =0.
\end{align}
\)

と導けます。

\(k\lt 0\)の時

まず、\(\varepsilon\)の符号が正である場合を考えます。

kが正の時と異なる点は、漸近で収束する領域が複素x平面の上半面であるという点です。

図の左のような閉経路を考えます。

同様に計算を行うと
\(
\begin{align}
\oint \frac{e^{-ikx}}{x}dx&=
\lim_{R\to\infty}\left[\int_{-R}^{-\varepsilon’}+\int_{\Gamma_+}+\int_{\varepsilon’}^{R}+\int_{\Gamma’}\right]
\frac{e^{ikx}}{x}dx\\
&=\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx+\lim_{R\to \infty}\int_{\Gamma’}\frac{e^{-ikx}}{x}dx
\end{align}
\)
今、閉経路内に特異点を含まないので、左辺\(\oint \frac{e^{-ikx}}{x}dx=0\)、また遠方ではゼロに近づくので、\(\int_{\Gamma’}\frac{e^{ikx}}{x}dx=0\)です。よって、
\(
\begin{align}
\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx =0.
\end{align}
\)

また、特異点周りを正方向に回る場合、閉経路内に特異点を含むので、
\(
\begin{align}
\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}e^{-ikx}dx =2\pi i.
\end{align}
\)

となります。

\(k= 0\)の時

\(k=0\)では被積分関数は漸近で収束しません。なので漸近で大きく回った時の値がゼロになりません。

ここからは誤解を恐れずに計算します。
今求めたい積分は
\(
\int_{-\infty}^{\infty}\frac{1}{x}e^{-i0x}dx =\int_{-\infty}^{\infty}\frac{1}{x}dx
\)

です。実軸上の積分がこのページで書いた
\(
\begin{align}
&\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx \nonumber \\
&=\lim_{\varepsilon\to +0} \frac{1}{2}\left[\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx +\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}e^{-ikx}dx \right]
\end{align}
\)
で定義されているのだとすれば、
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}dx =\lim_{\varepsilon\to +0} \frac{1}{2}\left[\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}dx +\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}dx \right]
\end{align}
\)

です。簡単のため、特異点を含まない閉経路を取れば第1項は
\(
\begin{align}
\lim_{R\to \infty}
\left[\int_{-R}^{-\varepsilon’}+\int_{\Gamma_+}
+\int_{\varepsilon’}^{R}+\int_{\Gamma’} \right]\frac{1}{x}dx =\oint \frac{1}{x}dx=0
\end{align}
\)

なので、
\(
\begin{align}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}dx +\int_{\Gamma’} \frac{1}{x}dx =0
\end{align}
\)

となります。外側を回る経路\(\Gamma’\)を通った積分は
\(
\begin{align}
\int_{\Gamma’} \frac{1}{x}dx &=\int_0^\pi \frac{1}{re^{i\theta}}ire^{i\theta} d\theta \\
&=\pi i
\end{align}
\)

と計算できてしまうので、代入すれば
\(
\begin{align}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty} \frac{1}{x+i\varepsilon} dx=-\pi i
\end{align}
\)

となります。また、同様にして
\(
\begin{align}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty} \frac{1}{x-i\varepsilon} dx=\pi i
\end{align}
\)

を得ます。
よって、実軸上の積分の式に代入すれば
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}dx =\lim_{\varepsilon\to +0} \frac{1}{2}\left[(-\pi i) +\pi i \right] =0
\end{align}
\)

を得ます。すなわち、フーリエ変換後の\(k=0\)の値はゼロということです。

計算結果のまとめ


以上をまとめると、3つの関係式
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
-2\pi i~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)

\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x-i\varepsilon}dx =
\left\{
\begin{aligned}
0~~~(k\gt 0)\\
2\pi i~~~(k \lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)

\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}e^{-i0x}dx=0
\end{align}
\)

を得ることが出来ました。
これを実軸上の積分の式に代入すると、
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx =-\pi i ~\text{sgn}(k)
\end{align}
\)

となります。ここで、\(\text{sgn}(x)\)は符号関数で、
\(
\begin{eqnarray}
\text{sgn}(x) =
\left\{
\begin{aligned}
-1~~~(x\lt 0)\\
0~~~(x=0) \\
1~~~(x\gt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
で定義される関数です。