ナルミヤの備忘録(仮)

ナルミヤが学んだことなどを書き記していくブログ(方向性模索中。)

【領域図示】2変数陰関数不等式の領域を図示する【Python + matplotlib】

はじめに

数値解析のレポートで、修正オイラー法や、4次のルンゲクッタ法の安定領域の図示にかなり時間食われたので、今後そんなことが無いようにここに書き残して置くことにしました。

やりたいこと

今回自分がやりたかったことは、修正オイラー法の安定条件である

$$ \left|1+ \lambda h+ \frac{(\lambda h)^2}{2}\right| < 1 \hspace{0.4cm} \lambda h \in \mathbb{C} \tag{1} $$ や4次のルンゲクッタ法の安定条件である $$ \left|1+ \lambda h+ \frac{(\lambda h)^2}{2} + \frac{(\lambda h)^3}{6} + \frac{(\lambda h)^4}{24}\right| < 1 \hspace{0.4cm} \lambda h \in \mathbb{C} \tag{2} $$

の領域を図示することでした。
ちなみに、以下では$z= \lambda h=x+iy \, \, x,y \in \mathbb{R}$として表して行きます。

できたこと

式(1)や式(2)から下の図が出来上がりました。
(1)が左で、(2)が右です。

f:id:buddasls54:20190105044332p:plainf:id:buddasls54:20190105044429p:plain
図1.安定領域のプロット:修正オイラー法(ホイン法)(左)とルンゲクッタ法(右)

手順

式(1)や(2)の領域を図示するには以下のような方法で行きます。

  1. 必要なライブラリのインポート
  2. 式(1)や(2)式の両辺を2乗して整理した式を$x,y$の条件式$f(x,y) < 0$として求める
  3. 手順2で求まった、陰関数型の条件$f(x,y) < 0$に対して図をプロットする

では、わかりやすいように、まず、式(1)から見て行きましょう。

実行環境はJupyterです。

手順1. 必要なライブラリのインポート

インポートしましょう

import matplotlib.pyplot as plt
import numpy as np
import matplotlib
matplotlib.use('Agg')
%matplotlib inline

from sympy import *

手順2. 両辺を2乗して整理した式を$x,y$の条件式$f(x,y) < 0$として求める

Sympyのシンボリックと、expand機能を使います。
細かいとことや、他のことが知りたい人はここをみるといいと思います。
実際のコードは以下のようになります。

x = Symbol('x', real=True) # xを実数のシンボリックと定義
y = Symbol('y', real=True) # yを実数のシンボリックと定義
z = Symbol('z') # zをシンボリックと定義

## z = x+iyと定義する ##
z = x + y*I

expr = 1 + z + (z**2)/2 # 修正オイラー法の安定条件の絶対値の中の式を変数定義
expr2 = expand(expr_z.conjugate()*expr_z) -1 # 修正オイラー法の安定条件の(左辺)^2-(右辺)^2を計算して展開
print(expr2)

これの出力は以下のようになりました。

x**4/4 + x**3 + x**2*y**2/2 + 2*x**2 + x*y**2 + 2*x + y**4/4

これで、

$$ f(x,y) = \frac{x^4}{4} + \frac{x^2 y^2}{2}+ \frac{y^4}{4} + x^3 + xy^2 + 2x^2 + 2x $$

とすることで、両辺を2乗して整理した式を$x,y$の条件式$f(x,y) < 0$として求めることができました。

手順3.手順eで求まった、陰関数型の条件$f(x,y) < 0$に対して図をプロットする

さて手順2で$f(x,y) < 0$という形までできたので、この陰関数表示されている領域をプロットします。
実際のコードは以下のようです。細かい説明はあとでします。

## プロットする領域を定める ##
max_range = 3
delta = 0.025 # 刻み幅
xrange = np.arange(-max_range, max_range, delta)
yrange = np.arange(-max_range, max_range, delta)
X, Y = np.meshgrid(xrange,yrange)

## Figureを作成 ##
fig = plt.figure(figsize=(5, 5))
# FigureにAxes(サブプロット)を追加
ax = fig.add_subplot(111)

## 軸の設定 ##
plt.axis([-max_range, max_range, -max_range, max_range])
plt.gca().set_aspect('equal', adjustable='box')


## 軸を中心の方に持っていく ##
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data',0))

## 軸にラベルをつける ##
ax.set_ylabel("$Im$",y = 1)
ax.set_xlabel("$Re$", x=1)

## 描画 ##
Z=(X**4)/4+(X**2)*(Y**2)/2 + (Y**4)/4 + X**3 + X*Y**2 + 2*X**2 + 2*X
plt.contour(X, Y, Z, [0], colors='black', corner_mask=True) # 境界線をプロット
plt.plot(X[Z < 0], Y[Z<0], color='gray',alpha=0.5) # 内部をgrayで塗る

## グラフを保存 ##
plt.savefig("./figure/heun_stb.png")  # figureディレクトリがないとエラーになるのでその時はfigureディレクトリを作ってあげること
plt.show()

これで次の図ができます。

f:id:buddasls54:20190105044332p:plain
図2.生成された図

それでは各行の説明をして行きます。

プロットする領域を定める

max_range = 3
delta = 0.025 # 刻み幅
xrange = np.arange(-max_range, max_range, delta)
yrange = np.arange(-max_range, max_range, delta)
X, Y = np.meshgrid(xrange,yrange)

この5行は、描画する範囲を決めて、その範囲内で、xrangeyrangeというdelta刻みの配列を作って、それを格子状に貼った二個の変数X,Ynumpy.meshgridによって生成しています。

軸の設定

## 軸の設定 ##
plt.axis([-max_range, max_range, -max_range, max_range])
plt.gca().set_aspect('equal', adjustable='box')


## 軸を中心の方に持っていく ##
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data',0))

## 軸にラベルをつける ##
ax.set_ylabel("$Im$",y = 1)
ax.set_xlabel("$Re$", x=1)

コメントの通り、軸の表示範囲と、軸の表示場所、そして、軸のラベルの設定などをしています。細かいことは、公式リファレンスへGO!()

境界線と境界線内部領域の描画

さて、ここが本題です。

#描画
Z=(X**4)/4+(X**2)*(Y**2)/2 + (Y**4)/4 + X**3 + X*Y**2 + 2*X**2 + 2*X
plt.contour(X, Y, Z, [0], colors='black', corner_mask=True)
plt.plot(X[Z < 0], Y[Z<0], color='gray',alpha=0.5)

まず、手順2で求められた$f(x,y)$の$x$を$X$に$y$を$Y$にしたものを$Z$に代入しておきます。それが1行目のコードです。
次に、どちらが先でもいいんですが、境界線の描画からして行きます。もちろんこれがなくてもいいんですが、あったほうが「それっぽい」。陰関数の境界線は、matplotlib.pyplotcontourメソッドを用いて書きます。これはもともと、標高線をかくものみたいです。
例えば、

plt.contour(X, Y, Z, [0], colors='black', corner_mask=True)

とすると、Zの値が0となる境界線を引いてくれるらしいです。ちなみに、

plt.contour(X, Y, Z, [0,1,-1], colors='black', corner_mask=True)

とかすると、Zの値が0と-1と1の境界線を引いてくれる。colors='black'を指定しなければ、標高せんの色はカラフルになります。corner_maskはなんでTrueにしたのか覚えてないです!w

そして、境界線内部をプロットして行きます。3行目ですね。

plt.plot(X[Z < 0], Y[Z<0], color='gray',alpha=0.5)

これには普通に、Z< 0という条件のもと生成されたnumpyのboolean配列を用いて作られたX,Yの部分配列をgrayでプロットしていくだけです。
下の軸が見えるようにalpha=0.5とすることで透明度をあげてます。

最後に保存する

plt.savefig("./figure/heun_stb.png")  
plt.show()

保存します。
なお、このパス指定だと、figureディレクトリがないとエラー吐かれますので、その時は大人しく、figureディレクトリを作りましょう。もちろん、保存先先のパスを変えてもいいです。

以上!

ルンゲクッタ法の安定領域にも同じことをしてみる

手順2.$f(x,y)<0$の形にする

$f(x,y)<0$の形にするために、$f(x,y)$を求めます。

from sympy import *

x = Symbol('x', real=True)
y = Symbol('y', real=True)
z = Symbol('z')

z = x + y*I

expr_rk4 = 1 + z + (z**2)/2 + (z**3)/6  + (z**4)/24
expr_rk4_power2 = sy.expand(expr_rk4.conjugate()*expr_rk4) -1
print(expr_rk4_power2)

出力は

x**8/576 + x**7/72 + x**6*y**2/144 + 5*x**6/72 + x**5*y**2/24 + x**5/4 + x**4*y**4/96 + x**4*y**2/8 + 2*x**4/3 + x**3*y**4/24 + x**3*y**2/6 + 4*x**3/3 + x**2*y**6/144 + x**2*y**4/24 + 2*x**2 + x*y**6/72 - x*y**4/12 + 2*x + y**8/576 - y**6/72
-31/36

でした。

$$ \begin{align} f(x,y) = \frac{x^8}{576} &+ \frac{x^7}{72} + \frac{x^6 y^2}{144} + \frac{5x^6}{72} + \frac{x^5 y^2}{24} + \frac{x^5}{4} + \frac{x^4 y^4}{96} + \frac{x^4 y^2}{8} + \frac{2x^4}{3} +\frac{x^3 y^4}{24} \\ & + \frac{x^3 y^2}{6} + \frac{4^3}{3} + \frac{x^2 y^6}{144} + \frac{x^2 y^4}{24} + 2x^2 + \frac{xy^6}{72} - \frac{xy^4}{12} + 2x + \frac{y^8}{576} - \frac{y^6}{72} < 0 \end{align} $$

ですね。これは手計算で求めたくねえ。。。

手順3.図示

max_range = 3
delta = 0.025
xrange = np.arange(-max_range, max_range, delta)
yrange = np.arange(-max_range, max_range, delta)
X, Y = np.meshgrid(xrange,yrange)

# Figureを作成
fig = plt.figure(figsize=(5, 5))
# FigureにAxes(サブプロット)を追加
ax = fig.add_subplot(111)

#軸の設定
plt.axis([-max_range, max_range, -max_range, max_range])
plt.gca().set_aspect('equal', adjustable='box')

# 軸を中心の方に持っていく
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data',0))

ax.set_ylabel("$Im$",y = 1)
ax.set_xlabel("$Re$", x=1)

#描画
Z2 = (X**8)/576 + (X**7)/72 + (X**6)*(Y**2)/144 + 5*(X**6)/72 + (X**5)*(Y**2)/24 + (X**5)/4 + (X**4)*(Y**4)/96 + (X**4)*(Y**2)/8 + 2*(X**4)/3 + (X**3)*(Y**4)/24 + (X**3)*(Y**2)/6 + 4*(X**3)/3 + (X**2)*(Y**6)/144 + (X**2)*(Y**4)/24 + 2*(X**2) + X*(Y**6)/72 - X*(Y**4)/12 + 2*X + (Y**8)/576 - (Y**6)/72

plt.contour(X, Y, Z2, [0], colors='black', corner_mask=True)
plt.plot(X[Z2 < 0], Y[Z2<0], color='gray' ,alpha=0.5)
plt.savefig("./figure/rk4_stb.png")
plt.show()

出力は

f:id:buddasls54:20190105044429p:plain
図3.出力されたルンゲクッタ法の安定領域

ちゃんとできました。
まあ、この方法、$f(x,y)<0$にした後に、Zに直接、$f(x,y)$を入力してるのすごい脳筋な方法なんですけどね。
でも変に考えずに、こっちの方が早かったので...。アホなので。

最後に

珍しく連投となりましたが、これは早く片付けたかったので、連投も構わず、連投しました連投。
さて、数値解析のレポートも終わったし、試験勉強と、もう一つのレポートしなきゃ。まずい...!!!
俺ガイルも読みたいのじゃ〜。

では今日はここまで。

参考

あけましておめでとうごさいます

新年の挨拶

あけましておめでとうございます。なるみやです。
再履修している数値解析のレポートやってて完全に遅くなりましたが、新年あけて2019年ということで去年の振り返りと、今年の抱負などを上げていきたいと思います。備忘録として。

去年の振り返り

1月から3月

去年、2018年は色々なことがありました。本当に色々なことが。
まず、初っ端の1月、前年に引き続き、やる気も出ずにレポートも書かずになろう系小説「デスマーチからはじまる異世界狂想曲」を無気力に読んでたら、東京に観測史上初の大雪が降った日に振ったと思ったら、当時付き合ってた彼女に振られ、その日のうちにこのブログにその時思ったことを書き残し3月に公開しました。(現在、その記事?はアーカイブ化して非公開になってます。)
しかし、そのおかげか、1月終わりに「このままじゃダメだ!とにかく試験を頑張らなければ!!」という心理状況から試験勉強を頑張ったのですが、大学入ってから初めて単位を落としたと思ったら、ボロボロとそこから10単位以上を落としました。でも、それでも20単位ほど取れたのは、あの時、試験勉強を頑張ったからだと思います。(いや、まじであれなかったら留年確定してた...)
その後、2月に1年間やっていたサークルの運営の引き継ぎをしました。代替わりはしましたが、3月の追い出しコンパまでは幹事を担当し、そこで食べもしないで卓を周るごとに飲んでいたら、無事粗相しました。初めての粗相。
そうそう、2月終盤には人生で今後も語れるほどのことがありました。5日間で火山全店舗スタンプラリーしたことです。...あ、「火山」とは、「石焼ラーメン火山」のことで、当時、北は仙台、南は湘南までの全21店舗で経営してる「石焼ラーメン」というジャンルを売り出してるラーメン屋さんです。普通に美味しい。そこで全店舗巡ってかつ、9種類の石焼ラーメンを全て食べたらJTB旅行券5万円分がもらえるというキャンペーンをしてました。いつも仲良くバカしてる仲間の一人が、いきなり「このスタンプラリー行こうぜ」と言い始めたので、バカな僕ともう一人はじゃあ行くかとなって栃木、福島、新潟、宮城、群馬、埼玉、神奈川の7県を5日で回るというアホな計画を立てました。移動は全てレンタカー。寝泊まりは我らが故郷快活クラブ(ネカフェ)。ネカフェのシャワーや近くの温泉とか行きました。新潟とか滞在時間2時間ぐらいです。45分が石焼ラーメン火山で、あとは移動()。なんだこれ。でも、何一つ気を使うメンツではなかったのでめちゃくちゃ楽しかったです。やばい語り始めると止まらない。いや、本当に楽しかった。また、語る気力があればこのブログでも別に語りましょうかね。備忘録として。

f:id:buddasls54:20190105030806j:plainf:id:buddasls54:20190105030814j:plain
スタンプラリーの様子(左)と、火山の石焼ラーメン(右)

あ、そうそう、3月には高校同期と草津温泉にも行きました。安い宿を見つけたし、いいお店で飲んだこともあり、こちらもめちゃくちゃいい旅行となりました。また行きたいですねえ。草津
そんな感じで、サークル行事と旅行と帰省で春休みは終わりました。学業面での進捗ゼロやん...やば。

4月から7月

4月になり、3年生となりました。3年になったら思ったより学校が大変そうで、もう最低でも3年生の間は恋愛はできないな...と思った矢先に彼女ができました。(ここに書いたら怒られるかな?)
カミングアウトした人々に、「お前、恋愛しないいうとったやんけ」とお叱りを受けました。すまん、彼女が可愛かったんや。
そんなプライベートな出来事とは関係なく学業イベントは進んで行きました。初めての連続でした。制御論や回路学などの理論から、現代制御論数値計算だったり、データ構造やアルゴリズムのためのプログラミング、電子工作や3DCADや3Dプリンターを使った実装など、まさに工学部の始まりといった感じでした。楽しかったです。6月になると、プロジェクト演習と呼ばれる、同期3~4人が各研究室に配属されてのプロダクト作成もしました。楽しかったですねえ。そうこうしてるうちに7月になり期末テストもぼちぼちに学校は夏休みに入りました。

8月から9月

8月!夏休みが始まりました。しかし、僕は駒場に通い詰めの8月を過ごしました。前述した、プロジェクト演習で作ったプロダクトが中途半端な出来だったので、延長戦としてもう一度作りたいと思っていたら、そのお許しが出たためです。初めてのハードとソフト全体を通したコンテンツ作成でした。不完全なものでしたが、これにかけていた1ヶ月間はとても楽しかったです。また、やりたいなぁ。
9月はバイトで終わりました。バイトで実装の連続。そして帰省。以上。語ることはありません。8月との差異よ。。。彼女との思い出はたくさんありますがね。それはここでは語れることではありません。

10月から12月

10月になり、夏休みが終わり学校が始まりました。新学期はじめ特有の意識の高さとストレスや課題のなさから学校に出ていましたが、徐々にそのバランスは消え、学校に行けなくなりました。寒くなってからは完全にお布団とお友達です。毎年この季節は本当に活動割合が減るんですよね。自分、まじ野生生物。
この3ヶ月は怠惰でした。触覚講習会とかいったり、自分でも触覚勉強したりとかはしましたが、実装もろくにできずに、レポートと受けてない授業の消化に追われ、それに終わる日々。実装しろや本当に。あの8月が完全に無駄になってますね本当。
あ、そうそう、でも11月に話すことがありました。私、2浪して東大に入って、一応、「多浪」ということになっているため、東大多浪交流会という団体に所属しているのですが、そこが駒場祭(東大の学園祭)で出した、多浪経験談が詰まった冊子「東大多浪」に寄稿させていただいたものが世に出回りました。ありがたいことに、これに対して好意的な意見ばかりでいた。しかし、ここで書いたことも4000字という字数制限のもとで書いたので書き足りませんでした。もともとそうこうした時は2万字だったし。これについても、このブログでおいおい書いて行けたらいいなと思います。
ちなみに、これ駒場祭限定販売で1冊300円で売ってるんですが、メルカリで2300円とかで売られてるのを見て、近々通販開始が決定しました。メルカリでやると思います。

f:id:buddasls54:20190105034051j:plainf:id:buddasls54:20190105034325j:plain
東大多浪交流会冊子「東大多浪」とそれがメルカリで転売されてる様子

あ、あと、10月に学科のイベントでSonyの本社を見学させてもらう機会がありました。すごく刺激的でした。特にR&Dのやつ。感動した。「こういう動きのアクチュエーターが欲しい→ない!→なら作ろう!」という発想とかを聞いて感動してました。でも惜しむらくは、オープンソースとして公開する気がないこと。あそこ、コンテンツづくりがそこまで上手くないから、ユーザーにコンテンツ作りしてもらえばいいのに。。。なまじなんでも作れるせいで、すごいもの作って終わりになってる気がする。Sony、本当勿体無い。まあ、俺が口出すようなことじゃないかもだけど。
12月にも語るイベント、1つありました。xR学生大忘年会です。あれは本当に面白かった。xRを志す学生が集まってお話したり、イベント?をこなしていくの、楽しかった...。そのあとの某inm先生を囲う会も楽しかった。あそこにいた方々の足元にでも及ぶようになるように頑張りたいと思いました。本当に。

utvirtual.tech

総括

さて、そんな感じでザーッと振り返ってきた2019年でしたが、僕にとっては非常に惜しい使い方をした一年間でした。確実に成長はしました。人脈も広げました。すごい人たちにたくさんお会いできました。でも、スキル面の成長はまだ全然足りず、また、精神的な面での成長も恋愛面以外をのぞいて見られず、実装したり創作するという面では本当に何もできませんでした。
しかし、よくよく考えて見たら、「実装できなかったー。創作できなかったー。」って言い出したのも2018年からです。そう考えてみると、一つ具体的な未来への階段は前に進めたのかな?そんな1年でした。それでも悔いは残る1年でした。特に創作面。
今年、自分にとってよかった点があるとすれば恋愛面ですね。とにかく、この一年間、恋愛面では、大変なこともありましたが、それ以上に楽しい思い出や大切な感情をもらいました。ありがたいですね。

2019年の抱負

はじめに

先日、年を越す前のこと、「青春ブタ野郎シリーズ」のアニメ12話とその次の最終話を見てました。その作中で、あるキャラクターがこういうセリフをいってました。

「私はね、咲太君。人生って優しくなる為にあるんだと思っています。昨日の私よりも今日の私が、ちょっとだけ…優しい人間であればいいなと思いながら生きています」

はっとしました。ずっと忘れていた、でも昔僕が大切にしていた価値観でした。ここ最近、学業面での至らなさから無力感しか感じておらず、すごく疲れた日々を送ってました。しまいには涙を流す日もありました。ここずっと、私は、学業面での成長以外の、大切にしていた価値観すら忘れて自分を追い込んでいました。
でも、現状、僕は追い込んだところで精神がすり減るだけ。それなら、元に戻ってこの価値観を大切にしようと思いました。でも、やはり、僕の人生をもう一段先に進めるために、「創作と実装」は必ず必要なことです。これを成し遂げなければ、ずっとこれに劣等感?を抱いたまま生きていくことになる。そんなの嫌です。(そういや東大に固執した時もこんな精神状況だったな。。。)
あと、周りにいる大切な人たちはより一層大切にして行きたいです。

抱負

ということで2019年の今のところの抱負は以下の2つです。

  1. 今日、昨日の自分より優しくある。2018年の僕より優しくある。
  2. 創作して実装する。そして、それを公開する。
  3. 僕にとって大切な人をより大切にする。

この3つは必ず成し遂げます。
ついでに、院試もしっかりやって行きたい研究室に行ってやるんです!思う存分研究してやるんです!

最後に

というわけで、2019年の抱負でした。
年始まって1週間以内に書けてよかった。
さて、この記事を2019年終わりに、満足げな顔で読めるように努力しますかね。

では、今日はこの辺で。

【Python】多項式代数方程式の解を一気に見つけるアルゴリズムの実装【数値計算】

はじめに

学校のレポートをやるために、$ z^{ n } + a_{ n - 1 } z^{ n - 1 } + \cdots + a_0 = 0 \, ( z \in \mathbb{ C } ) $という多項式代数方程式の全ての根を数値的に求める必要が出たのでPythonで実装しました。
Newton法とか、1個だけ求めるアルゴリズムは知っていはいましたが、全て一気に解ける方法はないものかと調べていたら...、あったあった。先日買った「数値計算の常識」という本に解き方が載っていました。この本いいですね。気になる人の為にAmazonのリンク貼っておきます。(「数値計算の常識」のAmazonリンク)

なお、Python数値計算・数値シミュレーションの優秀なライブラリであるScipyにも代数方程式をとくメソッドがありましたが、今回はNumpyのみを使って自作してみました。
調べてたら、Numpyにもありました。多項式代数方程式を解く方法。Python優秀すぎる。まじでライブラリ使わないで実装する意味のなさ...。
という訳で、Scipy、Numpyのライブラリを用いて実装する方法については、おまけに記しました。

解きたい方程式

解きたい方程式を確認しておきます。次の式(1)です

$$ p(z) = z^{ n } + a_{ n - 1 } z^{ n - 1 } + \cdots + a_0 = 0 \hspace{0.3cm} ( z \in \mathbb{ C } ) \tag{1} $$

連立方程式を解くアルゴリズム

今から記すアルゴリズム連立法と言われるものらしいです。
では、連立法の流れをみていきましょう。

連立法の流れ

連立法は、基本、反復法のアルゴリズムまんまです。

すなわち、$n$個の$\nu$次の近似解$z_1^{(\nu)}, ..., z_n^{(\nu)}$に対して、ある$n$次関数$\varphi_1(z_1,...,z_n),...,\varphi_n(z_1,...,z_n)$を用いて$\nu+1$次の近似解 $$ z_i^{(\nu +1)} = z_i^{(\nu)} + \varphi_i(z_1^{(\nu)}, ..., z_n^{(\nu)}) \tag{2} $$
で計算して、ある正数$\varepsilon$について、ある$i$が$|p(z_i^{(\nu})| < \varepsilon$となったら反復を停止するものです。

これには、以下の手順が含まれます。

  1. 初期値$(z_1^{(0)}, ..., z_n^{(0)})$を定める
  2. 反復関数$\varphi_i(z_1, ..., z_n)$の定めて、$\nu$近似解$(z_1^{(\nu)}, ..., z_n^{(\nu)})$から$\nu + 1$近似解$(z_1^{(\nu+1)}, ..., z_n^{(\nu+1)})$を作る。
  3. 作られた$\nu + 1$近似解に対して誤差を見積もり反復を停止するか否かを判定する

初期値$(z_1^{(0)}, ..., z_n^{(0)})$を定める

では、まずどう初期値を求めるのがいいでしょうか。
初期値は求める値、すなわち、方程式の根から近い方がいいです。そのため、根が存在する範囲をまず調べます。
そのために、根の重心$b$を求め、その重心から全ての根を含む最大半径$r$を得る必要があります。

根の重心は、解と係数の関係から、解を $ \alpha_{ 1 }, \alpha_{ 2 }, ..., \alpha_{ n }$ とすると、 $$ b = - \frac{ a_{ n - 1 } }{ n } = - \frac{ \alpha_{ 1 } + \alpha_{ 1 } + \cdots + \alpha_{ n } }{ n } \tag{3} $$ で求められます。

半径$r$の求めかたはテクニカルです。

すなわち、根の重心周りで$p(z) $を展開した式 $$ p(\zeta + b) = \zeta^n + c_{n-1} \zeta^{n-1} + \cdots + c_1 \zeta + c_0 \tag{4.1} $$ の係数を組立除法(あとで説明します。)(もしくはテイラー展開?)で求めて、実係数の方程式 $$ x^n - |c_{n-2}| x^{n-2} - \cdots - |c_1|x - |c_0| = 0 \tag{4.2} $$ を考えると、式(4)を$x^n$でくくると、 $$ x^n \{ 1 - ( \frac{|c_{n-2}|}{x^2} +\cdots + \frac{|c_1|}{x^{n-1}} + \frac{|c_0|}{x^n} ) \} = 0 $$ となり、 $$ f(x) = \frac{|c_{n-2}|}{x^2} +\cdots + \frac{|c_1|}{x^{n-1}} + \frac{|c_0|}{x^n} $$ は$x>0$で$f(x) > 0$でかつ単調減少であり、$x<0$で$f(x) < 0$でかつ単調増加ゆえ、$g(x) = 1$と$x = r > 0$で実数の範囲内でただ一つの共有点を持つことがわかります。ゆえに、式(4)にはただ一つの正の解があり、それを$r$とします。
このとき、式(1)の複素数根は根の重心$b$から半径$r$の円内に含まれます。数値計算的には、この$r$をNewton法で少し大きめの初期値から始めて得ることとなります。

そして、根の存在範囲を得られたら、次の初期値から出発するのがいいらしいです。(なんでかは知りません。調べてないです。すいません。)

$$ z_i^{0} := b + r \exp [j (\frac{2\pi (i-1)}{n} + \frac{3}{2n})] \tag{5} $$ ただし、$j$は虚数単位です。

反復関数の定め方

次に反復関数$\varphi_i(z_1, ..., z_n)$をいかに定めるかです。これこそ連立法のキモです。
今回は、連立法でも2次法と言われるものを用います。これは、Newton法の反復法の2次収束の考え方を背景に持ちます。
すなわち、Newton法での反復関数 $$ \varphi_i(z_1, ..., z_n) = - \frac{p(z_i^{(\nu)})}{p'(z_i^{(\nu)})} $$ において、$p(z) = \prod_{k =1 }^n (z- \alpha_k )$とすると、 $$ p'(z) = \sum_{k=1}^n \prod_{j \neq k} (z- \alpha_j ) $$ となるため、$z_i^{(\nu)}$が$\alpha_i$に十分近いとして、$z_i^{(\nu)}- \alpha_i \simeq 0$を用いると、 $$ p'(z) = \prod_{j \neq i} (z_i^{(\nu)}- \alpha_i) $$ となり、他の$j \neq i$についても$z_j^{(\nu)} \simeq \alpha_j$を用いると、 $$ p'(z) = \prod_{j \neq i} (z_i^{(\nu)}- z_j^{(\nu)}) $$ となるため、結局、反復関数を以下のようにおきます。 $$ \varphi_i(z_1, ..., z_n) := - \frac{p(z_i)}{\prod_{k = 1, \, k \neq i}^n (z_i - z_k)} \tag{6} $$

こうすると、反復列は2次収束します。

停止条件 : 誤差を見積もる(Smithの方法)

そうして作った反復列に対し、停止条件を考えるために、誤差の見積もりをします。

式(1)の解の近似値$z_1,...,z_n$があるとき、真の存在範囲を厳密に求める方法があります。この方法の根拠は長くなるので割愛しますが、大まかに次の方法で求められます。
まず、全ての$i$について、$|p(z_i)| < \delta_i$となるような$\delta_i$を計算して、それを用いて$\varphi_i$の上界$\tilde{\varphi_i}$を $$ \varphi_i(z_1,...,z_n) \le \tilde{\varphi_i}(z_1,...,z_n) := \frac{|p(z_i)| + \delta_i}{\prod_{j \neq i} |z_i - z_j|} $$ と算出しておきます。
そうして、$n$個の円、 $$ C_i := \{ z | |z_i + \varphi_i (z_1,...,z_n) - z| \le (n-1) \tilde{\varphi_i}(z_1,...,z_n) \} $$ を作ると、全ての根は$\cup_{i=1}^n C_i$の中にあります。
この$C_i$の半径$(n-1) \tilde{\varphi_i}(z_1,...,z_n)$が十分に小さい正数$\varepsilon$を用いて、 $$ (n-1) \tilde{\varphi_i}(z_1,...,z_n) < \varepsilon $$ となった時に、反復列は収束したとみなします。
ごちゃごちゃ言いましたが、実際には、$|p(z_i^{(\mu)})|$の計算をして、それが、十分に小さい正数$\varepsilon$を用いて、 $$ \forall i \in \{1,2,..., n\} \hspace{0.3cm} |p(z_i^{(\mu)})| < \varepsilon $$ となる時に反復を停止する方法でもいいです。

以上が連立法の流れとなります。

(補足) 組み立て除法

ここで式(2)を得るために用いる組立除法について軽く触れておきます。
組み立て除法とは、テーラー展開した時の係数を求めるアルゴリズムといっても過言ではなさそうです。

すなわち、 $$ \begin{array}{ccl} p(x) &=& p(\alpha) + p'(\alpha)(x-\alpha) + \frac{p''(\alpha)}{2!}(x-\alpha)^2 + \cdots tag{7} \\\ &=& c_0 + c_1 \zeta + c_2 \zeta^2 + \cdots \end{array} $$ に対して、係数$c_k = \frac{p^{k}(\alpha)}{k!}$を求めるのに、わざわざ微分してべき乗計算するのは計算コストが高いので、工夫をしましょう。式(7)より $$ \begin{array}{ccl} p(x) &=& p(\alpha) + (x-\alpha) ( p'(\alpha) +\frac{p''(\alpha)}{2!}(x-\alpha) + \cdots ) \\\ &=& p(\alpha) + q^{(1)}(x, \alpha) \end{array} $$ となるため、$p(\alpha)$は$p(x)$を$x-\alpha$で割った時の余りとなることがわかります。そのため、$p(x)$を$x-\alpha$で割って余りを求めることを考えましょう。この時、次の表を用いて求められます。(証明は割り算の筆算を書いてみるとわかります)
係数 $a_0 $ $a_1$ $a_3$ $\ldots$ $a_{n-1}$ $a_{n}$
$\times \alpha$ $a_0 \alpha$ $b_1^{(1)} \alpha$ $\ldots$ $b_{n-2}^{(1)} \alpha$ $b_{n-1}^{(1)} \alpha$
$a_0 $ $b_1^{(1)}$ $b_2^{(1)}$ $\ldots$ $b_{n-1}^{(1)} $ $b_{n}^{(1)} = p(\alpha)$
これにより、$p(\alpha)$は求められました。
また、$q^{(1)}(x, \alpha) = a_0 + b_1^{(1)}x^{n-2} + \cdots + b_{n-1}^{(1)}$となるため、この$q^{(1)}(x, \alpha)$に対して、再び、 $$ q^{(1)}(x, \alpha) = p'(\alpha) + (x-\alpha) q^{(2)}(x, \alpha) $$ とすると、$p'(\alpha)$は$q^{(1)}(x, \alpha)$を$(x-\alpha)$で割った時の剰余で、$q^{(2)}(x, \alpha)$は商。先ほどと同じようにして、それぞれ求められる。これを繰り返して行くことで、 $$ \frac{p^{k}(\alpha)}{k!} = b_{n-k}^{(k+1)} $$ として求められます。

実装

という訳で、2次の連立法で多項式方程式を求める場合、実装を次のような流れで進めます。

  1. 係数配列$a = [a_0, a_1,...,a_{n-1}]$で多項式を定義する
  2. Newton法の実装する
  3. 根の存在範囲を求める。

    • 根の重心を求める
    • 組み立て除法により、式(4.1)の係数$c_i$を求める
    • 式(4.2)の解をNewton法により求める
  4. 初期値を求めて反復法を開始して、根の配列$z = [z_0, z_1,..., z_{ n -1 }]$を求める。

このように進めて行きます。

配列coefficientsを係数にもつ多項式関数を作る

def getGfunc(coefficients):
    """N次式を作る
    coefficients = [a_0, a_1,...., a_{N-1}]として、a_0 + a_1 x^1 + ... + a_{N-1} x^{N-1}を得る。
    --Input---
    coefficients : ndarray. N個の要素をもつ
    
    -- return --
    g : coefficientsを係数とするzに関するN次式。
    """
    
    def g(x):
        func_val = 0
        for index in range(len(coefficients)):
            func_val += coefficients[index] * (x**(index))
        return func_val
    
    return g

複素数版Newton法の実装

Newton法で解を見つけるときは以下のコードでできます。複素数版Newton法です。

def newtonMethod(func, init_val, max_step = 50):
    delta_val = complex(0.0001, 0.0001) # 微小区間の幅
    
    # 以下Newton法
    for i in range(max_step): # 回数制限をかける
        f = func(init_val)
        df = (func(init_val + delta_val) - func(init_val - delta_val)) / (2*delta_val) # 微分値の取得
        init_val = init_val - f/df # 解の更新

    ans_x = init_val # 近似解
    ans_y = func(ans_x) # そのときのyの値(確認用)
    print('x : '+str(ans_x))
    print('y : '+str(ans_y))
    return ans_x


def func_y(x_dat):
    func = 2*(x_dat**2) - 3*x_dat - 4
    return func

if __name__ == '__main__':
    init_val = complex(4, 0) # 初期値
    print(newtonMethod(func_y, init_val))

最大半径を求める

組立除法で、式(2)の係数$c_i$を求め、さらに、式(2)の正の根$r$をNewton法で求めることを目指す。
def getRadiusBykumitateZyoho(coefficients):
    """組立除法とニュートン法を用いて最大半径を求める
     --Input---
    coefficients : ndarray. N個の要素をもつ係数配列.
    
    -- return --
    g : float. coefficientsを係数とするzに関するN次式の解に対して、
        解の重心からの全ての解が含まれる最大半径rを求める。
    """
    
    dim_n = len(coefficients) - 1
    
    # 根の重心
    root_centroid = coefficients[-2] / dim_n
    
    # rを求める際にできる関数の係数を格納するリスト
    new_coefficients = []
    
    b = 0
    for k in range(dim_n):
        for i in range(dim_n - k):
            b = coefficients[i] + b * root_centroid
        new_coefficients.append(-abs(b))
    
    new_coefficients.append(1)
    # 新しい関数を作る
    new_func = getGfunc(new_coefficients)
    
    radius = newtonMethod(new_func, 10)
    return abs(radius)

連立法

# zについてのN次の代数方程式をとく
def solveEquationNdim(coefficients, max_step=10**2):
    """連立法を用いて、代数方程式をとく"""
    
    dim_n = len(coefficients) - 1
    root_centroid = - coefficients[-2]/dim_n
    
    # 求める関数
    objectFunc = getGfunc(coefficients)
    
    # 最大半径を求める
    radius = getRadiusBykumitateZyoho(coefficients)
    
    # 初期値を設定する
    z_list = []
    for i in range(dim_n):
        z_list.append(root_centroid + radius*np.exp(1j *(2*np.pi*((i+1) - 1)/dim_n + 3/(2*dim_n))))
    
    
    def timesFunc(val_list, index):
        # 分母
        denominator = 1
        val_i = val_list[index]
        for j in range(len(val_list)):
            if j == index:
                continue
            denominator *= (val_i - val_list[j])
        return denominator
    
    # 求まった根の個数
    ans_num = 0
    
    # 更新された根を格納するリスト
    new_z_list = [0]*dim_n
    
    # 関数値を格納するリスト
    func_val_list = [0]*dim_n
    
    for i in range(max_step):
        for i in range(len(z_list)):
            fun_val = objectFunc(z_list[i])
            func_val_list[i] = fun_val
            
            denominator = timesFunc(z_list, i) # 分母を計算
            psi = - fun_val/denominator # 更新する値を計算
            new_z = z_list[i] + psi # 更新
            
            new_z_list[i] = new_z
            new_radius = new_z + abs(fun_val)/abs(denominator)
        z_list = new_z_list
        
        if (np.max(abs(np.array(func_val_list))) < 1e-8):
            print(i+1,"回で停止")
            return np.sort(np.array(z_list))
        
    print("収束しなかった")
    return np.sort(np.array(z_list))

おまけ

Scipy、Numpyのライブラリを用いて実装する方法について軽くみていきます。

Numpyのライブラリで解く方法

Numpyには多項式を扱うためのライブラリがあります。それを使うことで一瞬で解が求まるみたいです。すげえ。
Numpyには多項式を扱う際に、二つの方法があるみたいです。

1個目の方法 : poly1dメソッド

1個は、poly1dメソッドを用いる方法。以下のコードでできます。
以下のコードでは$x^{ 3 } - 2 x^{ 2 } + x - 2 = 0$を解いてます。

# Define a polynomial x^3 - 2 x^2 + x - 2.
p = np.poly1d([1, -2, 1, -2])

# print polynomial
print(p)

# Solve f(x) == 0.
print("roots: ", p.roots)

こうすると、

   3     2
1 x - 2 x + 1 x - 2
roots:  [ 2.00000000e+00+0.j -4.16333634e-16+1.j -4.16333634e-16-1.j]

と出力されます。

ちなみに、p(0)とかすると、$x=0$を代入した値を得られるみたい。

2個目の方法 : polynomialメソッドを使う

polynomialメソッドもべき乗多項式を扱うものです。実装例は以下のよう。
なお、先ほどと同じく、$x^{ 3 } - 2 x^{ 2 } + x - 2 = 0$を解いてます。

from numpy.polynomial.polynomial import Polynomial

# Define a polynomial x^3 - 2 x^2 + x - 2.
testEqF = Polynomial([-2, 1, -2, 1])
print("polynomial: ", testEqF)

# Solve f(x) == 0.
print("roots: ", testEqF.roots())

注意すべきは、先ほどと係数の指定する順番が真逆だということ。
出力は、

polynomial:  poly([-2.  1. -2.  1.])
roots:  [-1.83880688e-16-1.j -1.83880688e-16+1.j  2.00000000e+00+0.j]

となりました。根の表示は先ほどとは少し順番が異なりますね。

Scipyのライブラリでとく方法

こちらはScipyのシンボリック機能とsolveメソッドを使う方法。
お次は、$2x^{ 3 } + 3x - 5 = 0$を解いてみます。以下のコードででできます。

from sympy import *
x=Symbol('x')                  # 文字'x'を変数xとして定義

"""
方程式: solve
例: 2 * x **3 +3*x-5=0の根を求める
"""
sol=solve(2 * x **3 +3*x-5, x)
print(sol)

出力は、

[1, -1/2 - 3*I/2, -1/2 + 3*I/2]

Iというのは虚数を表すのでしょう。こっちの方が直感的。

最後に

今回の記事、思ったより長かった...。疲れた...。
まあ、でもいつかまた必要になると思うからね。。。。
でも、まじでライブラリ使わないならPython使う意味ほぼないし、C++あたりで実装し直そうかな...。

全く関係ないけど、この記事を書きながらスマブラやってました。楽しい。
キングクルール使いやすくて良いね。
ということで、今回はここまでです。バイバーイ。

参考文献

連立法

Numpyライブラリで解く

Scipyで代数方程式を解く

【Python】複素線積分の数値的解法 【数値計算】

はじめに

脳磁逆問題を解く際に複素線積分する必要があっただが、Pythonで複素線積分をしているものが見つからなかったので、記事を書こうと思いました。
まあ、複素積分って実質は線積分に過ぎないので、線積分が数値的に求められたら、それを複素数にしたら終わりです。

構成としては、

  1. 数値積分の基礎理論を触れる
  2. 閉曲線に沿った数値積分の理論について触れる
  3. 実装

って感じでいきます。

理論

数値積分の基礎

数値積分には、代表的に台形則や、中点則、シンプトン則があるが、それぞれ、積分を「求積問題」として、求める面積を細かく区間で分けて、どのように近似していくかという概念からきてるものであった。 それぞれの近似の概念としては以下のような感じ

  • 台形則...積分を狭い台形の足し合わせとして近似
  • 中点則...積分を微小区間の中点を高さの代表値とする幅の狭い長方形の足し合わせとして近似
  • シンプソン則...微小区間において関数$f(x)$を2次曲線で近似し、微小区間の端点とその中点とでの近似2次曲線上の点を用いて台形を作って積分を近似。(台形則の精度をあげている)

基本的に、下にいくほうが精度が高い(はず)

以下の図はシンプソン則の概念をわかりやすく図に起こしたもの。(重田 出 : 数値積分と数値微分(基礎)から引用)

f:id:buddasls54:20181212220615p:plain
図:シンプソ則について

閉曲線に沿った数値積分の理論

閉じた曲線$C$に沿う複素積分 $$ I := \oint_C f(z) \mathrm{d}z \tag{1} $$ について考えた時、まず、閉じた曲線を円に限る。(任意の閉曲線での線積分って円周での線積分で代用できなかったっけ?)
この時、文献1によると円に沿った線積分は上端も下端もなく、シンプソン則をはじめとした積分の上端と下端からの効果を補償したものだから、そんなことを考えない台形則で十分らしい。
ってな訳で、台形則だけで数値積分を行う。

という訳で、式(1)は

$$ I \simeq \sum_{k=0}^{N -1 } ( f( z_k ) + f(z_{ k -1 }) ) \frac{z_{k +1 } - z_k}{2} \tag{2} $$

で数値積分を行うが、この時、$z = z_0 + r\mathrm{e}^{i\theta}$は円周上と限られているので、極座標での数値積分を行うといい。
という訳で、結局考える積分

$$ I = \int_0^{2\pi} f(z_0 + r\mathrm{e}^{i\theta}) ri\mathrm{e}^{i\theta} \mathrm{d}\theta \tag{3} $$

となる。これに台形則を適用すると、

$$ I \simeq \sum_{k=0}^{N -1 } f(z_0 + r\mathrm{e}^{i\theta_k}) ri\mathrm{e}^{i\theta_k} \Delta \theta \tag{4} $$

となる。ただし、

$$ \Delta \theta = \frac{2\pi}{N}, \hspace{0.5cm} \theta_k = k\Delta \theta \tag{5} $$

とした。

理論は以上。

数値複素積分を実際にやってみる

さて、理論も済んだし、複素積分をやっていきましょ。題名の通りPythonを使います。理由は僕が一番使い慣れてるから。
という訳で実際のコードは以下のよう。

import numpy as np

# 複素積分
_int_step_num = 10**3
_int_center_point = 0 + 0j
_int_r = 1.1

def numerical_int( func, z0=_int_center_point, r=_int_r, Ni=_int_step_num ):
    """"関数fを複素積分する。
    z0は積分半径の中心。rは積分経路の半径。Niは刻み数。
    数値積分の値を返す。"""
    
    I = 0 # 積分値
    delta_theta = 2*np.pi/Ni
    
    for k in range(Ni):
        theta_k = k*delta_theta
        if np.isnan(func(z0 + r* np.exp(complex(0, theta_k )))):
            continue
        I += func(z0 + r* np.exp(complex(0, theta_k )))*   r*1j * np.exp(complex(0, theta_k )) *delta_theta
        
    return I

刻み数は実験的にやったら、刻み数が$10^{3}$ぐらいでも精度良さげという適当な理由で決めた。

実験

$1/z$で実験してみた。厳密解は$2\pi i$。

f:id:buddasls54:20181212224508p:plain
図:実験結果

良さげ。

以上!!!

参考文献

最後に

これはアドベントカレンダーのProny法の実装に繋がる準備としてブログ化しました。アンド、将来僕がまた複素積分するときに困らないように。 今月は、これと、

の更新が目標。(できるとは言ってない。) 今月は5回の更新を目指すぞー!おー!

では、今日はここまで!!バイバイ!!!

【逆問題】脳磁逆問題をProny法をといてPythonで実装【理論編】

はじめに

初めましての人は初めまして。そうでない人はお久しぶりです、なるみやです。
自分は、東京大学計数学科というところに所属しているのですが、これはその計数工学科と物理工学科からなる応用物理学門(略称AP)で企画されたAPアドベントカレンダー2018/12/08分のものです。

Prony法とは、大雑把に言えば、逆問題の一種に対する解法です。得られている結果から要因を求めてやろうというやつの磁気バージョン。例えば、脳磁場のポテンシャルエネルギーが測定できてる時に、その脳磁の源となる磁気双極子の個数と、その場所、そして磁気モーメントを求めてやろうというときに用いられる一つの解法です。

これができたら、頭のどこに双極子の源(「ソース」と言います)があるかわかるため、人の頭の中の活動が、頭の中を侵さずに計測できてすごくハッピーというわけです!!!すごい!!!!

構成

この記事は以下のように進めて行きます。

  1. 準備編:磁場の「複素ポテンシャル」という概念の導入と、それを用いた双極子の表し方。
  2. 脳磁逆問題の設定
  3. 逆問題の解析(多重極子展開と双極子パラメタの導入)
  4. Prony法を用いた解法(変数変換による解法)
  5. 3の数値計算(N次方程式の解の求め方)
  6. 考察(笑)

今回の記事では、1-4まで書いて行きます。長くなりそうなので。5-6は、次の記事で書いて行きます。(記事が出来次第リンクを貼ります。)
そうなんです!!今回プログラムは一切出てきません!!!タイトルにPythonとかあるのに!!!!
つまるところ、この記事は理論編と実装編との2部構成となり、これは、その理論編です。タイトルにも書いてあるしね。そのまんま。

導入が長くなりましたね、さて初めて行きましょー。

1. 準備:磁場の「複素ポテンシャル」という概念の導入

磁場のポテンシャル

最初に断っておきますと、当方、電磁気が苦手なので、変なところがあったらすいません。反省します。
話は、まず、Maxwell方程式の第一式、ガウスの法則から出発します。
一様空間中で原点に点電荷$q$がある時は、ガウスの法則より、 $$ \boldsymbol{\nabla} \cdot \boldsymbol{E} = \frac{q\delta(\boldsymbol{r})}{\varepsilon_0} $$ にしたがって電場が生じます。ただし、$\delta(*)$はディラックデルタ関数です。この時、ポテンシャル$V$、つまり、$\boldsymbol{E} =- \boldsymbol{\nabla} V$を考えると、 $$ - \boldsymbol{\nabla}^2 V = \frac{q\delta(\boldsymbol{r})}{\varepsilon_0} $$ となります。(これは、ポアソン方程式。)
このポアソン方程式の解は3次元だと $$ V = \frac{q}{4\pi \varepsilon_0 r} $$ 2次元だと、 $$ V = - \frac{1}{2\pi \varepsilon_0} \log r $$ と表されます。ここでは、この2次元平面でのポテンシャルについて話して行きます。

複素ポテンシャルについて

さて、先ほどのポテンシャル$V$を複素ポテンシャルに拡張して行きましょう。まず、点電荷の位置を$z=x+iy$とします。なお、ここから先、$ \varepsilon_0=1$と仮定して進めて行きます。いちいちこれを書いていくのも面倒だし。

このとき先ほど求めたポテンシャルは $$ V = -\frac{1}{2\pi} \log |z| $$ となります。このポテンシャルの実部として、複素数として拡張します。
この時に、虚部としてちょうどいいのってなんですかね。....そう、共役調和な関数(コーシーリーマンの関係式を満たす関数)だね!!そうすれば、拡張されたポテンシャルが正則なものになって都合がいい!!
というわけで、 $$ f(z) = u + iv $$ で、係数はひとまず置いておくと$u(z) = \log |z|$とし、 $$ \frac{\partial u}{\partial x} = \frac{\partial v}{\partial x} , \, \frac{\partial u}{\partial y} = - \frac{\partial u}{\partial x} $$ となる$v(z)$を求めると、それは$v = \arg z$なので、求める複素ポテンシャルは $$ f(z) = - \frac{q}{2\pi} ( \log |z| + \arg z ) = - \frac{q}{2\pi} \log z $$ となります。

複素ポテンシャルの物理的な意味

複素ポテンシャルの物理的意味を考えると、

  • $\log |z|$=一定...等電位線
  • $\arg |z|$=一定...電気力線

と考えられます。すなわち、複素ポテンシャルは、等電位線と電気力線の2つを1で表してるものと言えます。

双極子の複素ポテンシャル

次に、原点に$+q$を, $z=-\epsilon \in \mathbb{C}$に$-q$の電荷をおくと、この時に生成される複素ポテンシャルは、重ね合わせの法則から、

$$ \begin{array}{cl} f(z) &= - \frac{q}{2\pi} \log z + \frac{q}{2\pi} \log (z+\epsilon) \\ &= - \frac{q}{2\pi} \log z + \frac{q}{2\pi} (\log z + \log (1 + \frac{\epsilon}{z}))\\ &= \frac{q}{2\pi} \log (1 + \frac{\epsilon}{z}) \\ &= \frac{q}{2\pi} \left( \frac{\epsilon}{z} - \left( \frac{\epsilon}{z} \right)^2 + ... \right) \end{array} $$

ただし、最後の等式は$\log (1+x)$のテイラー展開を用いました。これを$\epsilon \to 0$とする、すなわち、二つの単位電荷を無限に近づけると、

$$ f(z) \to \frac{m}{2\pi z} $$

となります。これが、双極子の複素ポテンシャルです。ただし、$ q \epsilon =m $とした。$m$は物理的にはモーメントを意味します。
同様にして、
逆向きの双極子を無限に近づけると四重極子$1/z^2$が生じ、 逆向きの四重極子を無限に近づけると八重双極子$1/z^4$が生じ...となっていきます。
一般に、2n重極の複素ポテンシャルは$\frac{1}{z^n}$と書けます。

原点が$z_0$の場所に移ったときは、$z \to z-z_0$とすればいいです。

2. 脳磁逆問題の設定

さて、準備が終わったところで問題設定です。これは、最初に言ったのと同じです。

逆問題:
観測場所$z$として、二次元領域(人の頭とかを考えてください)$\Omega$中の$z_1,...,z_N$の中にモーメント$m_1,...,m_N$の双極子がある。 このとき、$ \Omega $の境界線、$ \partial \Omega $ 上で、複素ポテンシャル $$ f(z) = \sum_{k=1}^{N} \frac{m_k}{2\pi (z-z_k)} $$ を観測した。このとき、$z_1,...,z_N, m_1,...,m_N,N$を推定せよ。

というものです。

さて、この解き方の発想としては、いくつかあります。
今回はその中でもProny法をご紹介します。
まずは、この逆問題が解きやすくなるように少し工夫をして行きましょー。

3. 逆問題の解析(多重極子展開と双極子パラメタの導入)

本当に、$f(z)$だけで、$N$個のソースの位置とかモーメントとか求められるのお??とお思いになるでしょう?
でもできるんです!(ジョン・カビラ風)
問題が難しいのなら、解きやすい方向に持っていってあげればいいんです!!
では、そのやり方を見て行きましょう。

ローラン展開と複素ポテンシャル

さて、この逆問題を解いていくために、まず、$f(z)$のローラン展開をします。なんでローラン展開を考えるかの理由は、すぐ後にわかります。
ローラン展開とは、簡単に言っちゃえば、テイラー展開の拡張版です。(ローラン展開についてわからない人は、この記事(EMANの物理学)をどうぞ。)

さて、$f(z)$を原点周りでローラン級数展開をすると、

$$ \begin{array}{cl} f(z) &= \frac{1}{2\pi} \sum_{k=1}^{N} \frac{m_k}{z(1-z_k / z)} \\\ &= \frac{1}{2\pi} \sum_{k=1}^{N} \frac{m_k}{z} \left( 1 + \frac{z_k}{z} + (\frac{z_k}{z})^2 +\dots \right) \\\ &= \frac{1}{2\pi} \left( \frac{\sum_{k=1}^{N} m_k}{z} + \frac{\sum_{k=1}^{N} m_k z_k}{z^2} + \dots \right) \end{array} $$

2個目の等式には、$|z_k /z| <1$を使いました。$z$はセンサの観測位置で、$z_k$はソースの位置です。
さて、最後にできてた式を見ると、結局、右辺は、

(原点においた双極子) + (原点においた四重極子) + ...

となる。つまり、ローラン展開は、多重極展開であると言えます

観測した$f(z)$から双極子パラメタを定める

ソースの位置$z_1,...,z_N$やモーメント$m_1,..,m_N$を直接求めるのは難しいので、$f(z)$を観測して双極子パラメタを決めて、そのパラメタから、求めて行きたいと思います。
まずは、$f(z)$を$\partial \Omega$上でそのまま積分してみることとします。今回、$\partial \Omega$を単位円と仮定します。
このとき、

$$ \int_{\partial \Omega}\frac{dz}{z} = \int_0^{2\pi} \frac{i \mathrm{e}^{i\theta}}{\mathrm{e}^{i\theta}}d\theta = \int_0^{2\pi} d\theta = 2\pi i \\\ \int_{\partial \Omega}\frac{dz}{z^2} = \int_0^{2\pi} \frac{i \mathrm{e}^{i\theta}}{\mathrm{e}^{i2\theta}}d\theta =\int_0^{2\pi} \frac{i }{\mathrm{e}^{i\theta}}d\theta = 0 $$ 以下、同様にすると、 $$ \int_{\partial \Omega}\frac{dz}{z^2} = \int_{\partial \Omega}\frac{dz}{z^3} = \dots = 0 $$

となるので、結局、

$$ \begin{array}{cl} \int_{ \partial \Omega} f(z) dz &= \frac{1}{2\pi} \int_{\partial \Omega} \left( \frac{\sum_{k=1}^{N} m_k}{z} + \frac{\sum_{k=1}^{N} m_k z_k}{z^2} + \dots \right) dz\\\ &= i \sum_{k=1}^{N} m_k \end{array} $$

となります。積分しただけで、モーメントの和の$i$倍が出ちゃいました!!!
それに、この計算過程をみると、要は、$1/z$の項が残り、それ以外は0になるっぽい!!!
というわけで、お次は、$f(z)$に$z$をかけたものを積分にしてみると...

$$ \begin{array}{cl} \int_{ \partial \Omega} f(z) z dz &= \frac{1}{2\pi} \int_{\partial \Omega} \left( \sum_{k=1}^{N} m_k z + \frac{\sum_{k=1}^{N} m_k z_k}{z} + \frac{\sum_{k=1}^{N} m_k z_k}{z^2} \dots \right) dz\\\ &= i \sum_{k=1}^{N} m_k z_k \end{array} $$

次は、モーメントにソースの位置をかけたものの和が出てきました!!!いい感じですね!!!
これを繰り返していくと、結局

$$ c_n := \frac{1}{i} \int_{ \partial \Omega} f(z) z^n dz = \sum_{k=1}^{N} m_k z_k^n $$

が得られます。この$c_n$をパラメータとして、逆問題を解いていけそうですね。

3章のまとめ

$c_n$というパラメータのおかげで逆問題は以下のように書き換えられました!

逆問題の違う表現:
$$ c_n := \frac{1}{i} \int_{ \partial \Omega} f(z) z^n dz \hspace{0.5cm} (n=0,1,...,N) \tag{1} $$ が与えられたとき、 $$ \sum_{k=1}^{N} m_k z_k^n = c_n \tag{2} $$ を満たす、$z_k,m_k \, (k=0,1,..,N)$および$N$を求めろ。

というわけで、次はこれを解くことを目指して行きます。

Prony法を用いた解法(変数変換による解法)

さて、では、書き換えられた問題を解いて行きましょー。解き方は色々あるみたいです。知ってるだけで4通りあります。(全て演習の授業で説明されましたものです。)
$N$個を厳密に推定する方法はありますが、計算コストが高いことなどがあり、$N=10$などと仮定して位置やモーメントを求めたのちに、モーメントが小さいものはカウントしないという手法を取るのが実用的みたいです。ここでも、その手法に乗っ取り、$N$は既知として進めて行きます。

まず単純に考えてみる

単純に書き換えられたものを考えてみると、それは以下の非線形連立方程式を解くこととなります。 例えば$N=2$のとき

$$ \left( \begin{array}{ccc} 1 & 1 \\\ z_1 & z_2 \\\ z_1^2 & z_2^2 \\\ z_1^3 & z_2^3 \end{array} \right) \left( \begin{array}{c} m_1 \\\ m_2 \end{array} \right) = \left( \begin{array}{c} c_0 \\\ c_1 \\\ c_2 \\\ c_3 \end{array} \right) $$

でも、この連立方程式をそのまま解くのは大変なので、ここは少し工夫を試みます。

視点を変えてみる(変数変換)

さて、ここからトリッキーな感じになってきます。
突然ですが、

$$g(x)= \prod_{i=0}^N (x- z_i)$$

とおくと、$g(x)=0$の根は$z_i$となるため、

$$ g(x)= x^N + a_{N-1} x^{N-1} + \dots + a_1 x + a_0 $$

としたときの係数$a_i$が定まれば、そのN次の方程式$g(x)=0$の根を求めることで$z_i$が求められることとなります。
これは、実質、$z_i$を$z_i$の基本対称式$a_i$に置き換えたことに相当します。(それは、解と係数の関係からわかることです。)

そして解法

次なるターゲットは$a_i$です。この$a_i$さえ求めることができれば、話はほぼ片付きます。Prony法も最後の段階です。

さて、$z_i$は$g(x)=0$の根なので、全ての$i=1,...,N$で、当然$g(z_i)=0$です。
よって、

$$ \sum_{i=1}^{N} g(z_i) m_i z_i^{k}= 0 \hspace{0.5cm} (k=0,1,...,N-1) \tag{*} $$

です。左辺を考えると、

$$ \begin{array}{cl} (\mbox{左辺}) &= \sum_{i=1}^{N} m_i ( z_i^{N+k} + a_{N-1} z_i^{N+k -1} + \dots + a_0 z_{i}^{k} ) \\\ &= \sum_{i=1}^{N} m_i z_i^{N+k} + a_{N -1 } \sum_{i=1}^{N} m_i z_i^{N+k -1 } + \dots + a_0 \sum_{i=1}^{N} m_i z_i^{k} \\\ &= c_{N+k} + a_{N-1}c_{N+k -1 } + \dots + a_0 c_{k} \end{array} $$

となります。よって式(*)は、

$$ c_{ N +k } + a_{ N-1} c_{ N+ k -1 } + \dots + a_0 c_k = 0 \hspace{0.5cm} (k=0,1,...,N-1) $$

となります。これを$k$について連立させると...

$$ \left( \begin{array}{ccccc} c_0 & c_1 & c_2 & \dots & c_{N -1 } \\\ c_1 & c_2 & & & \vdots \\\ c_2 & & \ddots & & \vdots \\\ \vdots & & & \ddots & \vdots \\\ c_{N -1 } & \ldots & \ldots & \ldots & c_{2 N -2 } \end{array} \right) \left( \begin{array}{c} a_0 \\\ a_1 \\\ \vdots \\\ a_{N -1 } \end{array} \right) = \left( \begin{array}{c} c_N \\\ c_{N +1 } \\\ \vdots \\\ c_{2 N -1 } \end{array} \right) $$

という、線形連立方程式ができます。ということで、この連立方程式を解くことで$g(x)$の係数$a_i$が求まります。
ちなみに、ここで出てきた行列、

$$ \left( \begin{array}{ccccc} c_0 & c_1 & c_2 & \dots & c_{N -1 } \\\ c_1 & c_2 & & & \vdots \\\ c_2 & & \ddots & & \vdots \\\ \vdots & & & \ddots & \vdots \\\ c_{N -1 } & \ldots & \ldots & \ldots & c_{2 N -2 } \end{array} \right) $$

はハンケル行列という名前がついており、$N \times N$の行列ですが、$2N -1$個のパラメータだけで行列が確定されます。
これによって、$g(x)$が求められたので、$g(x)=0$の方程式を解くことで$z_i$が求まります。
また、逆問題の元の式(2)より、

$$ \left( \begin{array}{cccc} 1 & 1 & \dots & 1 \\\ z_1 & z_2 & & z_N\\\ \vdots & \vdots & & \ddots & \vdots \\\ z_1^{N -1 } & z_2^{N -1 } & \ldots & z_N^{N -1 } \end{array} \right) \left( \begin{array}{c} m_1 \\\ m_2 \\\ \vdots \\\ m_{N} \end{array} \right) = \left( \begin{array}{c} c_0 \\\ c_1 \\\ \vdots \\\ c_{N -1 } \end{array} \right) $$

の関係式がわかっていて、先ほど、$z_i$が求められたので、この線形連立方程式を解いてしまば$m_i$が求められるということになります。
以上がProny法の理論です。(多分)

最後に

というわけで、今回はここまでです。 お疲れ様でした!!!(僕も疲れた!)
この記事としては、珍しくデスマス調にした気がしますが、なれませんね。このブログ、あくまでも備忘録として使ってるので。

元々は、我らが学科のコース長で、逆問題の権威でいらっしゃる奈良先生が、授業でお話して宿題にもなった「Prony法を数値計算でとけ」というものに対して、「東大多浪」という冊子の執筆をしていたため、満足にできなかったことから、そのリベンジの意味で書いてます。あと、ネットで調べた時にあまりにも日本語の記事が少なかったため、残して置いたら、あとあとProny法を実装する人の力になれるかなと思ったのもあります。(ほんとか?)
なお、ここでは理論よりも、実装(数値計算)の法に重きを置きますので、正直、Prony法については大雑把に説明して終わります。どうしても、Prony法について気になる方は「数学セミナー」の5-7月の数理のクロスロードを読むことをお勧めします。

さて、今回の記事はどうだったでしょうか?僕ら計数工学科は普段こんな勉強しています!

...ごめんなさい、嘘つきました。全てが全てこんな内容やってるわけじゃないです。これは授業で取り扱った中でもトップレベルに面白いことでした。しかし、それでも、計数工学科の本質とする学問はこういうものにあるところにあると思います。
そもそも、計数工学科の「計数」とは、計測と数理を足し合わせてできた言葉です。言うなれば、計測数理工学科です。
だから、成り立ちは、計測対象に数理モデルを立ててエンジニアリングをしていく学問を扱うところ、それが計数工学科の期限です。
今では、それに加えてアクチュエータなど、計測だけでなく制御なども行なっていますが、それの本質も、働きかける対象を計測して数理モデルを立てることです。
僕も最近、全ての本質はセンシングにありと思うようになりました。というのも、僕の標榜とする、「人工現実感」は人に現実ではないものを現実と感じさせることがゴールだと思っています。ということは、人がどのようにして「現実」を感じ、その人が「現実感」を得る感覚をいじれるようになれば、 「人工現実感」は達成されると思うからです。すなわち、人のセンシング機構を知り、それを制御できれば現実は作れる!!!...はずです。ですので、この僕の考えと、計数工学科の理念はすごくマッチしていて、ここ最近は、特に、計数にきてよかったなあ...と思っています。いやぁ、過去の俺、グッジョブ!!!

というわけで、なんか、計数についての所感を語ってしまいました。隙あれば自分語り。恥ずかしいですね。
長さも1ブログ記事としては、かなり長くなってきたので、ここら辺で引き上げたいと思います!! あくまでも、この記事、自分としては次の記事となる数値計算がメインとなるのですがそこまでちゃんと書けるかな?() それでは、ここで退散します!!ではでは、またどこかで!!!

参考