第 5 章 結論
D. CONTIN プログラムについて
92
93
リーマン・ルベーグの補助定理を表す式(D.1)と式(D.3)より、
F s d
F s A dy b
a k
b a k
k ( ) ( ) lim ( ) ( ) sin( ) … (D.5) が成立することは明らかである。
したがって、式(D.3)および式(D.4)において s(λ)を求めようとした場合、任意の大きさの振幅 A とある程度大きな角振動数ωを持ったs(
)Asin(
)も解の一つとして考えられうるというこ とを意味する。したがって、解として考えうる要素はある意味無限大であり、解を一意に定める ことは困難を極める。以上のように、数値逆ラプラス変換解析の解は一意に定まらず、高度な逆問題(inverse problem) かつ非適切な問題(ill-posed problem)となっていることがわかる。また、以上の問題は測定誤差を 含む積分方程式を解くということ自体に付き纏う問題である。積分方程式を解く上でも考えうる 解を絞るために境界条件を与えるように、本研究における数値逆ラプラス変換解析においても解 の存在範囲を狭めるための妥当な条件を制約として課す必要がある。
本研究に用いるプログラムである CONTIN は、このような制約を容易に課すことができ、かつ 安定解と適当なモデルとの間の妥協点を自動選択することが可能である。数値逆ラプラス変換解 析を達成する上で CONTIN は重要なツールであり、本研究でも CONTIN を用いて研究を行った。
なお、CONTIN 以外にも、数値逆ラプラス変換解析が可能なプログラムとして FTIKREG[D1]や FROG がある。次節では、数値逆ラプラス変換解析において重要な Tikhonov regularization につ いて説明を行う。
D.2 Tikhonov regularization
[24]観測値の行列をy、誤差に対するそれぞれの測定点に関する重みの行列をω、作用素行列をA、 未知入力行列をxとするとき、最小二乗評価式V(0) は、
minimum )
( ) 0
( yAx 2
V
… (D.6)で表される。ここで、先験的知識を以下の形で考慮することにより、式(D.6)の解の存在範囲を限 定することができる。
minimum )
( )
( )
( yAx 2 S x
V
… (D.7)な お、S(x)は解を 安定化 させるた めの x の関数 、αは 正の定 数であ り適切 化パ ラメー タ (Regularization Parameter)と呼ばれる。
例えば、xがxaに近い値をもつことが事前にわかっている場合には、S(x)として∥x－xa∥2を選 べば、式(D.7)は、
minimum )
( )
( yAx 2 xxa 2
V
… (D.8)となり、式(D.8)は式(D.7)にx = xaという制約を重み付きで課した形となる。
以上のように、先験的知識により考えうる適切化項を重み付きで最小二乗項に付加し、解の存 在範囲を狭める手法を Tikhonov regularization と呼ぶ。
94
D.3 CONTIN の基本事項
本節では、CONTIN の基本事項を説明する。なお、詳細は参考文献[20]-[22]を参照されたい。
CONTIN による数値逆ラプラス変換解析は、前節で述べた Tikhonov regularization によって達成 される。式(D.3)において、任意のNL個のLkiの重みのオフセットβiを付加すると、
k N
i i k i b
a k k
L
L d
s F
y
1
) ( )
( … (D.9)
となる。CONTIN は上記のような第一種のフレドホルム積分方程式やヴォルテラ積分方程式を解 くことができ、式(D.9)の数値積分を以下の線形代数方程式に変換する。
k N
i i k i N
m
m m k m k
g L
L s
F c
y
1 1
) ( )
( … (D.10)
ここで cmは求積法の重み関数である。数値逆ラプラス変換解析を行う際には、NLを 1 として単一 のオフセットβを仮定し、Fk(λm)を exp(-λmtk)に置き換えて解析する。式は以下の式(D.11)のよう に表される。
1
( ) exp( ) ( ) , 1
Ng
k k m m k m k t
m
y f t c t s B k N
… (D.11) また、式(D.11)は、以下の行列式(D.12)で表すことができる。1 1 1 1 1 1
1
1 1
1 1
exp( ) exp( ) exp( ) 1 ( )
( )
exp( ) exp( ) exp( ) 1 (
( )
( ) exp( ) exp( ) exp( ) 1
g g
g g
t t t g g t
m m N N
k m m k N N k
k
N N m m N N N N
c t c t c t s
f t
c t c t c t s
f t
f t c t c t c t
1
)
( )
g
t
m
k N
N
s B
… (D.12)
式(D.12)中の各測定誤差εkは未知である。式(D.12)の行列式をy = Ax +εとおくと、CONTIN は Tikhonov regularization により式(D.13)の評価関数V(α)を最小化する最適解を求める。
1 2
2 2
( ) 2( ) minimum
V M yAx rRx … (D.13)
CONTIN では、式(D.13)の重み付き最小二乗項に付加された適切化項として、以下の式(D.14)を適 用できる。
2 b 2
a[ ( )]s d
r Rx … (D.14)
式(D.14)の適切化項は、制約化条件として滑らかさの指標を付加することを意味する。よって、
適切化パラメータαが大きくなればなるほど、得られる s(λ)の関数は滑らかになる。CONTIN は 制約化項の重みαを変更しながら最小化を実施し、最適な重みαおよび最適解を決定する。
以上のように、CONTIN は安定解と適当なモデルとの間の妥協点を自動選択する。つまり、こ のような原理の難しさはある程度排除し、ユーザとしては容易にプログラムを扱うことが可能と なっている。次節では、CONTIN で数値逆ラプラス変換解析を実施する際の具体的な制御パラメ ータの設定値について説明する。
95
D.4 CONTIN による逆ラプラス変換アルゴリズムの実現
CONTIN を扱うためには、CONTIN の制御パラメータを設定する必要がある。解析に必要とな る設定値を以下に列挙する。詳細は参考文献[20]-[22]を参照されたい。
--- yk : QTS の V-t 特性の V の各要素
tk : QTS の V-t 特性の t の各要素
NINTT 1. … 等間隔に位置する tkの集合数 NSTEND NT TSTART TEND … t の要素数(tkの数)、開始点、終了点
cm : 求積法の重み関数
IQUAD 3. … Simpson’s rule を指定 Ng : 式(D.10)の第一項の Ng
NG … グリッド数 NL : 式(D.10)の第二項のオフセットの数
NLINF 1. … NL=1 により単一のオフセットを仮定 Lki : 式(D.10)の第二項の重み
λm : グリッドの指定
IGRID 1. … GMNMX(1)と GMNMX(2)の間を等間隔に設定 (Ng点生成, IGRID 2 の場合は対数間隔に設定)
Fk : 積分方程式のカーネル関数
ラプラス変換の作用素となるように以下の値を設定する。
) exp(
) ,
( m tk fm mR23 R21tk mR22
F
… (D.15)上記の式は、USERK により評価されるカーネル関数の一般形である。ここで、R23=0、R21=1、 R22=1、
そして fm=1 のとき、式(D.15)は以下のカーネル関数となる。
) exp(
) ,
(
mt
kt
k mF
… (D.16)式(D.16)はラプラス変換のカーネル関数に対応する。よって、上記の条件を満たす値を設定する ために以下の設定が必要となる。
LUSER 3 -1(.FALSE.) … fm=1
RUSER 16 0. … R16=0, R21=1.0
IUSER 10 2. … R23=0, R22=1.0, R21=R202
なお、論理制御変数は、TRUE として 1、FALSE として－1 のみ許される。
数値逆ラプラス変換解析においては、解は非負であるので以下の設定を行う。
NONNEG 1. … 非負
また、ラプラス変換における積分区間をGMNMX(1)≦GMNMX(2)により定義する。
GMNMX 1 … 下限値 GMNMX 2 … 上限値
---