最終更新: 2017/3/8

kvライブラリの応用例

柏木 雅英
本ページでは、kvライブラリを少し複雑な実問題に応用した例をいくつか紹介する。

2-transistor circuit equation which has 5 solutions

Yusuke Nakaya, Tetsuo Nishi, Shin'ichi Oishi, and Martin Claus: "Numerical Existence Proof of Five Solutions for Certain Two-Transistor Circuit Equations", Japan J. Indust. Appl. Math. Volume 26, Number 2-3 (2009), 327-336.

の論文にある例。 2つのトランジスタを含む回路の動作点(解)の数は従来最大3つであると 考えられていたが、この論文で5つの解を持つ例が示された。 その方程式は以下の通り。未知数は\(V_1, V_2, V_3, V_4\)。 \[ \begin{pmatrix} 1 & -\alpha_r & 0 & 0 \\ -\alpha_f & 1 & 0 & 0 \\ 0 & 0 & 1 & -\alpha_r \\ 0 & 0 & -\alpha_f & 1 \end{pmatrix} \begin{pmatrix} \frac{I_s}{\alpha_f}(e^{\frac{V_1}{V_T}}-1) \\ \frac{I_s}{\alpha_r}(e^{\frac{V_2}{V_T}}-1) \\ \frac{I_s}{\alpha_f}(e^{\frac{V_3}{V_T}}-1) \\ \frac{I_s}{\alpha_r}(e^{\frac{V_4}{V_T}}-1) \end{pmatrix} + \begin{pmatrix} 2G_b + G_c & -(G_b + G_c) & -2G_b & G_b \\ -(G_b + G_c) & G_b + G_c & G_b & 0 \\ -2G_b & G_b & 2G_b + G_c & -(G_b+G_c) \\ G_b & 0 & -(G_b + G_c) & G_b + G_c \end{pmatrix} \begin{pmatrix} V_1 \\ V_2 \\ V_3 \\ V_4 \end{pmatrix} + \begin{pmatrix} G_c V_{cc} \\ G_b V_s - G_c V_{cc} \\ G_c V_{cc} \\ G_b V_s - G_c V_{cc} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \] \begin{align} \mathrm{where} \\ G_b &= 1/R_b \\ G_c &= 1/R_c \\ R_b &= 10^4 \\ R_c &= 5 \times 10^3 \\ V_{cc} &= -5 \\ \alpha_f &= 0.99 \\ \alpha_r &= 0.5 \\ I_s &= 10^{-9} \\ V_T &= 0.053 \\ V_s &= -0.64 \end{align} この方程式の5つの解を非線形方程式の全解探索プログラムで探索してみる。

このプログラムを実行すると、 
[4]([0.70356598164369243,0.70363160049497453],[-0.72265625,-0.72143622077293478],[0.74230000554159702,0.74233295475972683],[0.6198074492962462,0.61996680696900542])(ex)
[4]([0.70358963169344701,0.70358963169362677],[-0.72180712343566756,-0.72180712342886677],[0.74231647296775893,0.74231647296781967],[0.61988728360925959,0.61988728360955603])(ex:improved)
[4]([0.71987945704442757,0.71991981499962843],[-0.0034523581433403496,-0.00244140625],[0.7355023219369261,0.73552255604547645],[0.57546705397045905,0.57563759442400664])(ex)
[4]([0.71990164129087852,0.71990164129099532],[-0.003335867014043125,-0.0033358670081318586],[0.73551253992991871,0.73551253992997967],[0.57555293496204418,0.57555293496258942])(ex:improved)
[4]([0.74230139160766606,0.7423315678063569],[0.61981417508245506,0.61996014351137508],[0.70356604997670513,0.70362832229358996],[-0.72265625,-0.72149841028493844])(ex)
[4]([0.74231647296775959,0.74231647296781911],[0.61988728360926281,0.61988728360955337],[0.70358963169344879,0.70358963169362488],[-0.72180712343560172,-0.7218071234289427])(ex:improved)
[4]([0.7355024825637294,0.73552237646543628],[0.57546830409154858,0.57563623001366749],[0.71987991176740129,0.71991980140167178],[-0.003440782931533068,-0.00244140625])(ex)
[4]([0.73551253992991894,0.73551253992997923],[0.57555293496204606,0.57555293496258542],[0.7199016412908793,0.7199016412909951],[-0.0033358670140009599,-0.0033358670081407915])(ex:improved)
[4]([0.72925176532520441,0.72932695555387029],[0.47024015360242243,0.47265286297307325],[0.72925747404052254,0.72932241572766921],[0.47042364565773697,0.47250638247688682])(ex)
[4]([0.72928963256368528,0.72928963256369895],[0.47145516180306784,0.47145516180350878],[0.72928963256368528,0.72928963256369895],[0.47145516180306701,0.47145516180350833])(ex:improved)
ne_test: 295668, ex_test: 145259, ne: 147816, ex: 5, giveup: 0
のようになった。 すなわち、V1-V4を-10から10の範囲として探索したところ、 295668回の非存在テストと145259回の存在テストを行い、147816個の非存在領域と 5つの存在領域の特定に成功し、5つの解は

V1 V2 V3 V4
(1) [0.70358963169344701,0.70358963169362677] [-0.72180712343566756,-0.72180712342886677] [0.74231647296775893,0.74231647296781967] [0.61988728360925959,0.61988728360955603]
(2) [0.71990164129087852,0.71990164129099532] [-0.003335867014043125,-0.0033358670081318586] [0.73551253992991871,0.73551253992997967] [0.57555293496204418,0.57555293496258942]
(3) [0.74231647296775959,0.74231647296781911] [0.61988728360926281,0.61988728360955337] [0.70358963169344879,0.70358963169362488] [-0.72180712343560172,-0.7218071234289427]
(4) [0.73551253992991894,0.73551253992997923] [0.57555293496204606,0.57555293496258542] [0.7199016412908793,0.7199016412909951] [-0.0033358670140009599,-0.0033358670081407915]
(5) [0.72928963256368528,0.72928963256369895] [0.47145516180306784,0.47145516180350878] [0.72928963256368528,0.72928963256369895] [0.47145516180306701,0.47145516180350833]

のように精度保証付きで得られた。

Sturm-Liouville型固有値問題

以下のようなSturm-Liouville型固有値問題の 固有値\(\lambda\)と固有関数\(y(x)\)を精度保証付きで 計算する問題を考える。 (E. Haiere, G. Wanner, S. P. Norset: "Solving Ordinary Differential Equations I Nonstiff Problems", Springer Series in Computational Mathematics, Volume 8, 1993 のp.109から引用) \begin{align} & ((1 - 0.8\sin^2x)y')' -(x-\lambda)y = 0 \\ & y(0) = y(\pi) = 0 \end{align} \(y'=z\)とおいて変形すると、 \begin{align} y' &= z \\ z' &= \frac{(1.6\sin x \cos x)z + (x-\lambda)y}{1 - 0.8\sin^2 x} \\ y(0) &= y(\pi) = 0 \end{align} となる。shooting methodを行なうために、\(\lambda\)を未知関数とし、 その微分が\(0\) (従って定数関数となる)とする。 \begin{align} y' &= z \\ z' &= \frac{(1.6\sin x \cos x)z + (x-\lambda)y}{1 - 0.8\sin^2 x} \\ \lambda' &= 0 \\ y(0) &= y(\pi) = 0 \end{align} このままだと境界条件が1つ足りないが、\(y(x)\)は定数倍の自由度があるため、 \(y'(0) = 1\)と固定してしまう。 \begin{align} y' &= z \\ z' &= \frac{(1.6\sin x \cos x)z + (x-\lambda)y}{1 - 0.8\sin^2 x} \\ \lambda' &= 0 \\ y(0) &= 0 \\ y(\pi) &= 0 \\ z(0) &= 1 \end{align} こうすれば、\(y(0), z(0)\)固定で、\(\lambda(0)\)を変化させながら 初期値問題を計算し、\(y(\pi)\)を\(0\)にするようなshooting methodとして 定式化できる。 これを解くために作成したプログラムを示す。 これは、固有値の近似値 \(2.1224\) を与え、 その解を軽く5回ほどNewton法で改善した後、 Krawczyk法で解の存在を検証するものである。 更に、解の存在検証に成功したなら、改めてその(区間)固有値を初期値として 初期値問題を計算する。それが精度保証された固有関数である。 固有関数の値は\(2^{-7}\)間隔で密出力する。

実行結果は以下の通り。

すなわち、固有値は、 \[ [2.122442998078283,2.1224429980782956] \] に存在し、固有関数を描画すると

のようになる。

更に、固有値の近似として \begin{align} & 3.6078 \\ & 6.0016 \\ & 9.3773 \\ & 13.7298 \\ & 19.053 \end{align} を用いて同様の計算を行なうと、固有値は \begin{align} & [3.6078101471559307,3.607810147155961] \\ & [6.0016132242501019,6.0016132242501695] \\ & [9.3772846770430859,9.3772846770432033] \\ & [13.729809464569756,13.729809464569945] \\ & [19.053238318275418,19.053238318275711] \end{align} のように精度保証され、それらの固有関数を先の固有関数と併せて描画すると、

のようになった。

更に、プログラムをわずかに変更して固有値の全解探索をさせてみる。 探索範囲を[0,20]としてそこに含まれる全ての固有値を計算するプログラムは、

の通り。実行結果は、 
[1]([2.0986626091980888,2.161736165269347])(ex)
[1]([2.122442998078283,2.1224429980782951])(ex:improved)
[1]([3.6062683995352315,3.609976767107979])(ex)
[1]([3.6078101471559311,3.6078101471559605])(ex:improved)
[1]([6.0004425545043877,6.0028424238613436])(ex)
[1]([6.0016132242501019,6.0016132242501686])(ex:improved)
[1]([9.3771351052845197,9.3775358438725825])(ex)
[1]([9.3772846770430859,9.3772846770432033])(ex:improved)
[1]([13.72974816853938,13.729896530042041])(ex)
[1]([13.729809464569756,13.729809464569942])(ex:improved)
[1]([19.053166776104859,19.053276214137953])(ex)
[1]([19.053238318275418,19.053238318275711])(ex:improved)
ne_test: 128, ex_test: 6, ne: 56, ex: 6, giveup: 0
のようになった。これが[0,20]に含まれる全ての固有値であり、すなわち 先に近似値を元に計算した固有値が第1-第6固有値であることを示している。