トップ :: G 物理学 :: G06 計算;計数




【発明の名称】 連立方程式解法
【発明者】 【氏名】井上 靖秋
【目的】 解曲線を効率よく追跡し、現実の回路解析において使いやすい連続法を実現する。
【構成】 予測子修正子法において解曲線を追跡している場合に、予測子と、修正子とのなす角度が所定の基準角度より大きいか否かが各解を求める度に判断される。そして、角度が基準角度より大きい場合には、その基準となった既に求められているｊ＋１番目の解は別の軌跡に移動してしまった等、追跡を失敗してしまった可能性が大きいとしてそのｊ番目の解を放棄し、一つ前のｊ番目の解から新たなより小さなステップ幅でｊ番目の解を求め直す。
【特許請求の範囲】
【請求項１】連立方程式をＢＤＦ法を用いてホモトピー方程式の解曲線を追跡することによって解く方法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、予測子と修正子とのなす角度φj +1を算出し、この角度φj +1が所定値より大きいか否か判定する判定ステップと、前記判定ステップにおいて、前記角度φj +1が前記所定値より大きいと判断された場合には、解曲線の追跡のステップのｊ＋１番目の解を求めるステップをより小さなステップ幅によって再実行し、ｊ＋１番目の解を新たに求め直すステップと、を含むことを特徴とする連立方程式解法。
【請求項２】前記請求項１記載の連立方程式解法において、前記角度φj +1は、次式で求められることを特徴とする連立方程式解法。
φj +1＝２arcsin（‖W j +1‖／２）
W j +1＝ξ1 （ｓj +1，ｓj ）−ξ1 （ｓj ，ｓj-1 ）
ここで、‖・‖は、ユークリッドノルムを表し、ξk （ｋ＝０，１，２，３，・・・）は次のようなｋ階の差分商を表す。
【数１】

【請求項３】連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、さらに、前記ｊステップのｓ＝ｓj の近傍ｓP ＝ｓj ＋ｈP において、前記ｊステップに至る後方値を用いて局所的に解曲線を多項式で近似する近似ステップと、前記近似ステップにおいて求めた解曲線に基づいて、前記ｊステップからステップ幅ｈP 進んだ点における角度φ（ｈP ）を求める予測角度ステップと、前記予測角度ステップにおいて求められた角度φ（ｈP ）が所定の許容角度値に等しくなる場合の前記ステップ幅ｈP を算出する予測ステップ幅算出工程と、を含み、前記予測ステップ幅算出工程で求められたステップ幅を用いて、解曲線上のｊ＋１番目の解が求められることを特徴とする連立方程式解法。
【請求項４】連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法であって、追跡中の各解の算出がニュートン・ラフソン法により算出される連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、ｊ＋１番目の解の算出の際に、ニュートン・ラフソン反復によって解が非収束であるときに、非収束に対処する対処ステップを実行する非収束判断ステップと、前記ステップ幅があらかじめ定められているステップ幅の下限値より大きい場合には、ステップ幅を小さくしてから、ｊ＋１番目の解を新たに算出し直し、あらかじめ定められているステップ幅の下限値と等しい場合には、既に実行したｊ番目の解を求めるステップをより小さなステップ幅で再実行し、ｊ番目の解を求め直す前記対処ステップと、を含むことを特徴とする連立方程式解法。
【請求項５】連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法であって、追跡中の各解の算出がニュートン・ラフソン法により算出される連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、ｊ＋１番目の解の算出において、ニュートン・ラフソン反復によって解が非収束である時に、ニュートン・ラフソン反復が発散する場合は、第一対処ステップを実行し、ニュートン・ラフソン反復が非発散非収束である場合には、第二対処ステップを実行する非収束判断ステップと、ｊ＋１番目の解の算出における前記ステップ幅があらかじめ定められているステップ幅の下限値より大きい場合は、ステップ幅を小さくしてから、ｊ＋１番目の解を求め直し、前記ステップ幅の下限値と等しい場合には、既に実行したｊ番目の解を求めるステップをより小さなステップ幅で再実行し、ｊ番目の解を求め直す前記第一対処ステップと、枢軸再選択を行ってから、再びニュートン・ラフソン反復を行ってｊ＋１番目の解を求める前記第二対処ステップと、を含むことを特徴とする連立方程式解法。
【請求項６】連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法であって、追跡中の各解の算出がニュートン・ラフソン法により算出される連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、ｊ＋１番目の解の算出において、ニュートン・ラフソン反復によって解が非収束である時に、ニュートン・ラフソン反復が非発散非収束である場合に対処ステップを実行する非収束判断ステップと、枢軸再選択を行ってから、再びニュートン・ラフソン反復を行ってｊ＋１番目の解を求める前記対処ステップと、を含むことを特徴とする連立方程式解法。
【請求項７】請求項１又は２の連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、さらに、ｊ＋１番目の解を求める際における、ｊ番目の解からｊ＋１番目の解までのステップ幅が所定のステップ幅下限値以上に設定されるステップ幅設定ステップと、前記ステップ幅設定ステップで求められたステップ幅を用いてｊ＋１番目の解を求める求解ステップと、を含み、前記ステップ幅設定ステップは、固定の下限値と、ｊ番目の解を求める際に用いられたステップ幅の値と、に基づいて前記ｊ＋１番目の解を求める際の所定のステップ幅下限値を算出する下限値設定ステップ、を有することを特徴とする連立方程式解法。
【請求項８】請求項１又は２の連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、さらに、ｊ＋１番目の解を求める際における解の許容誤差を、ｊ＋２番目の解を求める際のステップ幅下限値より小さく設定する許容誤差設定ステップ、を含むことを特徴とする連立方程式解法。
【請求項９】連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法において、前記ＢＤＦ法の修正子方程式の各方程式のうちステップ幅による拘束条件を形成しているステップ幅方程式に対し、ステップ幅の値で両辺を除算するステップと、前記ステップ幅方程式に対し、所定の定数を両辺に乗ずることにより、前記修正子方程式の各方程式のうち前記ステップ幅方程式以外の方程式群に対応するヤコビ行列の絶対値の大きさと、前記ステップ幅方程式に対応するヤコビ行列の絶対値の大きさとを近づける絶対値近似ステップと、を含むことを特徴とする連立方程式解法。
【発明の詳細な説明】【０００１】
【産業上の利用分野】本発明は、回路設計の支援装置として使用される回路シミュレータや、回路解析装置等において行われている大規模な連立方程式の解法に関するものである。特に、連続法の予測子修正子法における解曲線の追跡の改良に関するものである。
【０００２】
【従来の技術】近年、回路設計の支援ツールとしてＳＰＩＣＥ系の汎用回路シミュレータが広く用いられている。この汎用回路シミュレータは例えば、Nagel L.:"SPICE2, Acomputer program to simulate semiconductor circuits", Memorandum No.ERL-M520, Univ. of California, Berklay (May 1975) （以下、文献１と呼ぶ）に記載されている。しかし、このような回路シミュレータにおいては直流動作点解析において解が一定の値に収束しない場合がしばしば生じていた。このような非収束問題の有力な解決策として、従来用いられていたニュートン・ラフソン法（以下、ＮＲ法という）より収束性の高い連続法（若しくはホモトピー法）と呼ばれる一群の方法が幅広く知られるところとなった。それに伴い、これらの方法についての研究成果も多数報告されている。例えば、小島政和：”相補性と不動点−アルゴリズムによるアプローチ”、産業図書（昭和５６年）（以下、文献２と呼ぶ）や、Garcia C. B. and Zangwill W. I. :"PATHWAYS TO SOLUTIONS, FIXEDPOINTS, AND EQUILIBRIA", Prentice-Hall (1981)（以下、文献３と呼ぶ）等にこれらに関する記載が掲載されている。これらの方法は、方程式系を連続的に変形しながら、その方程式系が描く解曲線を追跡するものであり、原理的には大域的な収束を保証できるという特徴を有する。この特徴については、上記文献２や３だけでなく、Allgowcr E.L.and Georg K.:"Numerical Continuation Methods:An introduction", Springer 一Verlag（l990）（以下、文献４と呼ぶ）や、山村清隆、堀内和夫：“非線形回路解析におけるホモトピー法の収束性について”信学論（Ａ）．Ｊ７２−Ａ，１，pp．l56 −l59 （l989一0l）（以下、文献５と呼ぶ）等がある。
【０００３】一般に連続法は解曲線の追跡の仕方によって、区分的線形化法と、予測子修正子法とに大別される。区分的線形化法は、例えば、Ohtsuki T., Fujisawa T. and Kumagai S.:"Existence theorems and a solution algorithm for piecewise一linear resistor networks", SIAM J. Math. Anal., 8. 1, pp.69 −99（Feb.l977 ）（以下、文献６と呼ぶ）や、Yamamura K. and Horiuchi K. :"A globally and quadratically convergent algorithm for solving non1inear resistive networks", IEEE Trans. Comput.−Aided Des. of Integrated Circuits & Syst., CAD9. 5, pp.487 −499 （May 1990）（以下、文献７と呼ぶ）等に記載されている。また、予測子修正子法は、例えば、Cermak I. A. :"DC Solution ofnonlinear state −space equations in circuit analysis", IEEE Trans. Circuit Theory, CT−l8. 2. pp.3l2 −3l4 （March 197l）（以下、文献８と呼ぶ）や、Ushida A. and Chua. L. O.:"Tracing solution curves of nonlinear equations with sharp turning points", Int. J. Circuit Theory & Appl., l2, pp.l−2l（l984）（以下、文献９と呼ぶ）等に記載されている。
【０００４】回路解析の分野、例えば、ＬＳＩの設計に用いられ実用化されてきたのは後者、すなわち予測子修正子法によるものが多い。このような実用化については、上記文献８だけでなく、井上靖秋：“連続法に基づく統計解析アルゴリズム”．信学論（Ａ），Ｊ７４−Ａ，１，pp．63−7l（1991‐01）（以下、文献１０と呼ぶ）や、井上靖秋：“解曲線追跡回路を用いた非線形回路の直流解析”，信学論（Ａ），Ｊ７４−Ａ，１１，pp．1647‐l655（l99l‐ll）（以下、文献１２と呼ぶ）や、Chang M.-C.,Chern J.-H. and Yang P. :"Efficient and robust path tracing a1gorithm for DC convergence problem", Proc. IEEE ISCAS'93. pp.l635 一l638（May 1993）（以下、文献１４と呼ぶ）や、Melville R.C.,TrajkovicL., Fang S.-C. and Watson L.T.:"Artificial parameter homotopy methods for the DC operating point problem", IEEE Trans. Comput.−Aided Des. of Integrated Circuits & Syst., CAD−l2, 6, pp.86l −877 （June l993 ）（以下、文献１５と呼ぶ）等に記載されている。この予測子修正子法はさらに、超平面を用いる平面法と、超球面を用いる球面法とに分けられる。このうち、球面法は、アルゴリズムが簡潔で、既存の汎用回路シミュレータに容易に実装可能であるという特徴を有する。これらについては、上記文献９や１２等の他に、He C. and Ushida A. :"A modified predictor−corrector tracing curve algorithm for solving nonlinear resistive circuits", Trans. IEICE. E74, 6, pp.l455−l462（June l991 ）（以下、文献１１と呼ぶ）や、井上靖秋：“予測子修正子法を用いた解曲線追跡アルゴリズムの幾何学的解釈”，信学論（Ａ），Ｊ７５−Ａ，１１，pp．l682−l690（l992−ll）（以下、文献１３と呼ぶ）に記載されている。さらに、関連する文献として、Yamamura K.:"Simple algoriths for tracing solution curves, IEEE Trans. Circuits & Syst., CAS−40, 8, pp.537 −54l （Aug-l993）（以下、文献１６と呼ぶ）、井上靖秋: “大規模回路の直流勤作点解析法”，信学論（Ａ），Ｊ７７−Ａ，3 ，pp．388 −398 （l994−03）（以下、文献１７と呼ぶ）、山村清隆: “ホモトピー法／Katzenelson 法ベース回路シミュレータの可能性”，第７回、回路とシステム軽井沢ワークショップ論文集，pp.13-l8（l994-04 ）（以下、文献１８と呼ぶ）や、Yamamura K.:"Improvingthe convergence of spherical algorfthms for tracing solution curves", Trans. IEICE, E77A, 6, pp.1085-1088(June 1994) （以下、文献１９と呼ぶ）等が挙げられよう。
【０００５】本願発明者らは既に、ＳＰＩＣＥ系の汎用回路シミュレータに球面法を実装して、大規模回路の直流動作点解析の実用を開始している。このことについては、上記文献１７に記載されている。しかし、実用的な観点からはまだ解決すべき課題も多い。
【０００６】一般に予測子修正子法に基づく方法は、ステップ幅を大きくとって効率を高めようとすると、疑似解に収束してしまうなどの不安定要因があり（上記文献１３参照）、ステップ幅の制御に関して敏感な特性を有している。ステップ幅を解曲線の態様に応じて変化させるステップ幅制御アルゴリズムが従来から種々提案されている。本願発明者らのこれまでの経験則によれば（上記文献１２、１３、１７参照）、大規模回路では様々な不安定現象が生じ、従来のステップ幅制御アルゴリズムでは解曲線の追跡に失敗する場合が少なくない。そのため、アルゴリズムのパラメータを少し変更して解析を再実行するという試行錯誤を要する等の問題がある。その結果、連続法を広く実用化するために、より確実で効率の良いアルゴリズムの開発が望まれている。
【０００７】
【発明が解決しようとする課題】まず、超球面を用いた解曲線追跡法を概観し、大規模回路の解曲線追跡の問題点を指摘する。
【０００８】１. 解曲線追跡法と不安定現象１. １ 超球面を用いた解曲線追跡法非線形方程式【数２】

の求解に際して、まず、この方程式の代わりに次のような形式で表されるホモトピー方程式を構成する。
【０００９】
【数３】

このホモトピー方程式（２）の解曲線を追跡して、元の方程式（１）の解に到達する方法を連続法（又はホモトピー法）という。この方法については、上記文献２、３、４などに詳しい。ここで、【数４】

であり、【数５】

をホモトピー変数という。また、式（２）は、ｔ＝０で自明解【数６】

を持ち、ｔ＝１の時に、上記式（１）と等価になるように構成されている。
【００１０】式（２）の解集合は、適当な仮定の下では、Ｒn+1 空間で分岐しない滑らかなパスとループとを描く。このことについては、上記文献２、３、４等に記載されている。このような解曲線上の任意の点を解曲線の弧長【数７】

をパラメータとして、【数８】

と表す。このとき、【数９】

は次のような微分代数方程式を解くことにより決定される（上記文献９参照）、【数１０】

式（３）は、適当な陰的数値積分公式を適用して、予測子修正子法を用いて解くことができる。ｔ＝０の自明解【数１１】

から出発して解曲線を追跡してｔ＝１の解【数１２】

に到達すればそれが式（１）の解となる。
【００１１】式（３）を数値的に解く方法として、ＢＤＦ（後退差分公式）解曲線追跡アルゴリズムが知られている（上記文献９参照）。以下、本文において述べる１次のＢＤＦ解曲線追跡アルゴリズム（以下、簡単のため単にＢＤＦ法と呼ぶ）の概要について、簡単に説明する。このＢＤＦ法は、球面法の中でも、最も簡潔で、実現容易性が高いことが知られている（上記文献１３、１６参照）。そして、既に大規模回路の直流動作点解析に対してその有効性が上記文献１７において確認されている。ＢＤＦ法では、解曲線追跡のｊ＋１ステップｓ＝ｓj+1 において、次のような１次の予測子【数１３】

を用い、また、式（３）に１次の後退差分公式（例えば、Dahlquist G. and Bjorck A.:"NUMERICAL METHODS", Prentice−Hall（l974）（以下、文献２３と呼ぶ）に記載されている）を適用して得られた次のような修正子【数１４】

を用いる。ここで、【数１５】

更に、ｈj ＝ｓj −ｓj-1 はステップ幅を表す。
【００１２】ＢＤＦ法を用いた解曲線追跡の様子は図７を用いて次のように幾何学的に説明することができる（上記文献１３参照）。まず、解曲線に沿って適当な超球面を設定する。このとき、球面の中心は直前のステップの解の点に、また、球面の半径は現在のステップのステップ幅に設定される。次に、直前の２つのステップの解を結ぶ直線と球面の交点（図７において白丸で示されている）を予測子として、修正子方程式を用いて球面と解曲線の交点（黒丸）に修正する。この操作を繰り返して解曲線を順次追跡する。一般に、予測子修正子法に基づく解曲線追跡アルゴリズムでは、解曲線の方向の変化の緩急に応してステップ幅を調節し、可能な限り大きなステップ幅で追跡するのが計算効率の面からは望ましい。
【００１３】尚、図１３には、従来の解曲線の追跡の動作の概略を表すフローチャートが示されている。
【００１４】図１３のＳ１３−１においては、ｊやｓ0 などの値の初期設定が行われる。また、Ｓ１３−２においては、予測子であるｙj+1Pが予測される。Ｓ１３−３においては、ニュートン・ラフソン法によって解が求められる。Ｓ１３−４においては、ニュートン・ラフソン反復が収束するか否かについて判定される。この判定の結果、非収束であれば、ステップ幅を小さくして再びｊ＋１番目の解が求められる。すなわち、Ｓ１３−５においてステップ幅が小さく設定された後、再びＳ１３−２に処理が移行する。Ｓ１３−５におけるステップ幅の設定は、図１３に示されているように、そのとき使用されていたステップ幅の１／２、若しくはステップ幅の下限値、のいずれか大きい方が選ばれる手法等が考えられる。
【００１５】一方、Ｓ１３−４において、収束すると判定された場合には終了か否かが判断される。終了であると判断された場合には、最終的な解が求められたことになるが、まだ解曲線の追跡の途中であると判断された場合には、ｊをインクリメントした後、所定の演算により次の解を求める際のステップ幅が設定され、再びＳ１３−２に処理が移行する。
【００１６】１. ２ 修正子方程式の反復解法解曲線追跡の各ステップで修正子方程式を解く際にＮＲ法が用いられる（上記文献９参照）。上記式（５ａ）と（５ｂ）をまとめて【数１６】

の形式で表すと、ｊ＋１ステップにおけるｍ＋１回目の反復の近似解【数１７】

を与えるＮＲ反復式は【数１８】

で与えられる。ここで、【数１９】

は修正子方程式のヤコビ行列を表す。一般に、大規模回路ではヤコビ行列【数２０】

はスパース行列となる。
【００１７】修正子方程式のＮＲ反復の収束判定基準としては、【数２１】

などを用いるのが自然である（上記文献１７参照）。ここで、【数２２】

はそれぞれＮＲ反復の許容誤差を表す。後述する図１１のごとく、上記式（９ｂ）は近似解と超球面の距離を、また、式（９ｃ）は近似解と解曲線の距離を判定するために用いられる。所定の反復回数以内で上述の判定条件が満たされないとき、ＮＲ反復は非収束であるという。
【００１８】ＮＲ反復において線形化された方程式を解く場合、回路シミュレータではＬＵ分解法が用いられる（上記文献１参照）。その際に、数値許算の丸め誤差を小さく抑え、更に、スパース行列の充填要素の発生を少なくするために、ＬＵ分解に先立って枢軸選択（最適順序付け（pivoting）ともよばれる）が行われる。このような枢軸選択の方法として、Ｍａｒｋｏｗｉｔｚの方法（Markowitz H.M.:"The elimination form of the inverse and its application to linear programming", Management Sci., 3, pp.255−269 （l957）（以下、文献２２と呼ぶ）に記載されている）が広く用いられている（上記文献１参照）。枢軸選択とは、具体的には行又は列の入れ替えによって行われる。枢軸選択は、大きく２つの目的で行われる。一つは、いわゆる充填（fill-in ）を防止するためである。一般に回路解析の際の修正子方程式のヤコビ行列はそのほとんどの要素の値が０である「スパース性」を有している。そして、このスパース性という条件のもとで非常に高速に連立方程式を解く手法が知られている。ところが、計算が進むにつれて、今まで０であった要素が次々に非ゼロ値へと変更されてしまう（充填）ことにより、計算効率が悪化してしまう恐れがある。このような場合、あらかじめ枢軸選択により、充填が生じないようにすることが一般に行われている。二つ目の目的は、対角要素（ピボット（pivot ）要素）が０に近くなってしまうこと等のいわゆる数値的な悪条件の状態を防止するためである。対角要素が０に極めて近くなっていくと、計算誤差が増大しがちになってしまう。このような２つの目的で枢軸選択は行われる。
【００１９】枢軸選択はかなりの手間を要することから、従来の回路シミュレータにおける過渡解折の場合などと同様に、予測子修正子法においても解曲線の追跡の最初のステップのＮＲ反復の開始時点で一回だけ枢軸選択を実行し、その後は、特に必要のない限りはじめに選択された枢軸を用いて計算を続行するのが自然である（上記文献１７参照）。
【００２０】１. ３ 疑似解収束現象予測子修正子法に基づく解曲線追跡法に共通した欠点として、修正子アルゴリズムが疑似解に収束して解曲線の追跡に失敗してしまうという疑似解収束現象が知られている（上記文献１３参照）。球面法におけるこのような疑似解収束の典型的な例が図８に示されている。この図８において、Γ1 は追跡すべき解曲線（パス）、Γ2 は他の解曲線（この例ではループ）を表し、更に、Ｋは超球面を表す。また、矢印は解曲線の方向を表す。この場合、球面と解曲線の交点、つまり修正子方程式を満足する解は全部で６個存在する。このうち、所望の解ａを真性解、残りのすべての解ｂ〜ｆを疑似解という（上記文献１３を参照）。解曲線の追跡のすべてのステップで、常に真正解ａに収束すれば、最終解まで到達できる。一般に、解くべき方程式の性質と用いるアルゴリズムにより、場合によっては、ｂ〜ｆのいずれの疑似解にも収束する可能性がある。このような疑似解に収束すると、解曲線の追跡に失敗して最終解が得られなくなることが多い。
【００２１】従来のアルゴリズムが疑似解に収束して、解曲線の追跡に失敗する例を図９に基づいて説明する。図９に示されている例は、解曲線Γ1 を追跡している例である。従来の球面法のステップ幅制御アルゴリズムでは、概念的に予測子と修正子の差Ｅが所定値内に収まるようにステップ幅が調整されている（例えば、上記文献９、１１〜１３、１７を参照）。すなわち、図９（ａ）に示されているように、Ｅが所定値より大きくなった場合に、その球面（破線）を棄却して、より小さな球面（実線）を用いて、予測子修正子アルゴリズムを再実行する。その結果、図９（ｂ）の破線の方向に追跡が進行するが、途中で疑似解に収束して、追跡すべき解曲線Γ1 から他の解曲線Γ2 に乗り移っていても、これを検出することができない。
【００２２】球面法におけるこのような疑似解収束現象には、以下のような様々なタイプがある。本願発明者らのこれまでの経験をもとに（上記文献１２、１３、１７参照）、実用的な回路で確認されている例が、図１０の（ａ）、（ｂ）〜（ｌ）に列挙して示されている。図１０に示されている現象は次にように分類することができる。
【００２３】
（１）帰還現象 折り返し帰還型（ａ）
半島周回帰還型（ｂ）（ｃ）
離島周回帰還型（ｄ）（ｅ）
（２）循環現象 離島・湾循環型（ｆ）
半島循環型（ｇ）
離島循環型（ｈ）（ｉ）
（３）準安定現象 半島周回型（ｊ）
半島側路型（ｋ）
離島周回型（ｌ）
このうち、準安定現象に属するものは最終解まで到達する可能性があるので、最終結果だけではそのような現象が生じたかどうかの見分けがつかない。また、帰還現象や循環現象が生じると最終解は得られない。
【００２４】このような疑似解収束現象のうち、すべての帰還現象はその現象に陥った点で修正子方程式のヤコビ行列の行列式の符号が反転することにより検出可能であるので（上記文献２〜４、１４参照）、その時点でステップ幅を適当に小さくすれば、このような現象に対し適切に対応することができると考えられる（上記文献１８参照）。但し、この際に小さくするステップ幅が十分に小さくないと、循環現象や準安定現象に移行することもある。例えば、（ｂ）は（ｇ）に、（ｄ）は（ｆ）や（ｈ）にそれぞれ移行する可能性がある。なお、ヤコビ行列の行列式の符号はＬＵ分解の結果を用いて容易に求めることができる（上記文献２３参照）。
【００２５】循環現象や準安定現象の場合でも、ヤコビ行列の行列式の符号が反転する場合は同様にして対応することができる。但し、図１０の（ｆ）、（ｇ）、（ｈ）に示すような場合に有効な対策手段は知られていない。このような疑似解収束現象は一般にステップ幅を十分に小さくすれば原理的には生じないようにすることが可能である。しかし、そのようなステップ幅を事前に予測する実用的な方法はまだ知られていない。またさらに、ステップ幅を常に小さくしたのでは計算効率が悪化してしまうという問題がある。
【００２６】１. ４ 修正子方程式の非収束現象修正子方程式のＮＲ反復が非収束の場合、適当にステップ幅を小さく取り直してＮＲ反復を再度実行する方法が用いられる（上記文献１２、１７参照）。しかし、大規模な系では、次に示すようにさまざまな非収束の状況があり、このような単純なステップ幅制御だけではＮＲ反復が収束しないことがたびたび発生する。
【００２７】（１）充填要素の多発ＬＵ分解の途中で枢軸要素が許容下限値に満たないといういわゆる枢軸エラーが発生して、ＮＲ反復が続行できなくなることがある。このような枢軸エラーの原因として、上記式（８）のヤコビ行列のｎ×ｎの主小行列がturning point で特異になること（上記文献２〜４、９参照）や、また、ＬＵ分解が数値的な悪条件（上記文献２３参照）に陥っていることなどが考えられる。
【００２８】このような場合には、枢軸選択を再実行するのが自然であるが（上記文献１７参照）、途中で充填要素が極端に多くなって枢軸選択が実用的な時間で終了しないことがある。
【００２９】（２）非発散非収束ＮＲ反復の誤差がある所定の上限値【数２３】

または、【数２４】

を超えた状態、つまり、【数２５】

のいずれかの状態を本文では「発散」と呼ぶ。また、ＮＲ反復がこのように発散するわけではないが、収束もしない状態を「非発散非収束」とよぶ。解曲線追跡の途中でいくらＮＲ反復を繰り返しても、このような非発散非収束の状態から抜け出せないことがある。このような非発散非収束現象はＬＵ分解が数値的に悪条件の状態に陥っているためと考えられる。
【００３０】この現象はステッブ幅を小さくしても解消しない。また、ＮＲ反復の許容誤差を厳しく設定するとこのような現象が生じやすくなる。解曲線追跡の途中で混入する誤差は最終解の精度に影響しないので、ＮＲ反復の収束判定基準を少し緩くすることも対策として考えられるが、その場合には次に述べる無解現象が生じやすくなってしまう。
【００３１】（３）無解予測子修正子法に基づく解曲線追跡アルゴリズムでば、ステップ幅の下限値ｈL とＮＲ反復の許容誤差【数２６】

との関係について特段の注意を要する。例えば、図１１に示されているように、直前のステップの許容誤差よりも現在のステップのステップ幅（球面の半径）が小さくなると、直前のｊステップの近似解が解曲線Γや球面Ｋj から離れてしまうことがある。その結果、次のｊ＋１ステップで解曲線Γと球面Ｋj +1とが交わらず、修正子方程式に解が存在しなくなる。このような無解現象が生ずると、ヤコビ行列の正則性は保証されなくなる。従って、ＮＲ反復の許容誤差（少なくとも）
【数２７】

は、球面の半径の下限値、つまりステップ幅の下限値ｈL （経験的にl0-4以下）に比べて十分小さな値（例えば10-5以下）に設定する必要がある。しかし、これは回路シミュレーションで普通に用いられる値（10-3程度）（上記文献１参照）に比べて非常に厳しい要求となる。このようにステップ幅の下限値に応じてＮＲ反復の収束判定基準が厳しくなると、反復回数の増大を引き起こし、前述の非発散非収束現象が生じやすくなってしまう。
【００３２】本発明は上記課題に鑑みなされたものであり、その目的は、解曲線をより確実に効率よく追跡し、現実の回路解析において使いやすい連続法を実現することである。
【００３３】
【課題を解決するための手段】本発明は球面法をもとに、解曲線をより確実に効率よく追跡するための実用的なステップ幅制御アルゴリズムを提案する。
【００３４】本発明は、上述した解曲線追跡の際に発生する間題を対策するための実用的なステップ幅制御アルゴリズムを提供するものである。本アルゴリズムは球面法に基づいており、特に上記文献１７で有効性が確認されているＢＤＦ法に適している。
【００３５】まず、第一の本発明は、上記課題を解決するために、連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、予測子と修正子とのなす角度φj +1を算出し、この角度φj +1が所定値より大きいか否か判定する判定ステップと、前記判定ステップにおいて、前記角度φj +1が前記所定値より大きいと判断された場合には、解曲線の追跡のステップのｊ＋１番目の解を求めるステップをより小さなステップ幅によって再実行し、ｊ＋１番目の解を新たに求め直すステップと、を含むことを特徴とする連立方程式解法である。
【００３６】第二の本発明は、上記課題を解決するために、上記第一の本発明の連立方程式解法において、前記角度φj +1は、次式で求められることを特徴とする連立方程式解法である。
【００３７】φj +1＝２arcsin（‖W j +1‖／２）
W j +1＝ξ1 （ｓj +1，ｓj ）−ξ1 （ｓj ，ｓj-1 ）
ここで、‖・‖は、ユークリッドノルムを表し、ξk （ｋ＝０，１，２，３，・・・）は次のようなｋ階の差分商を表す。
【００３８】
【数２８】

第三の本発明は、上記課題を解決するために、連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、さらに、ｊステップのｓ＝ｓj の近傍ｓP ＝ｓj ＋ｈP において、前記ｊステップに至る後方値を用いて局所的に解曲線を多項式で近似する近似ステップと、前記近似ステップにおいて求めた解曲線に基づいて、前記ｊステップからステップ幅ｈP 進んだ点における角度φ（ｈP ）を求める予測角度ステップと、前記予測角度ステップにおいて求められた角度φ（ｈP ）が所定の許容角度値に等しくなる場合の前記ステップ幅ｈP を算出する予測ステップ幅算出工程と、を含み、前記予測ステップ幅算出工程で求められたステップ幅を用いて、解曲線上のｊ＋１番目の解が求められることを特徴とする連立方程式解法である。
【００３９】第四の本発明は、上記課題を解決するために、連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法であって、追跡中の各解の算出がニュートン・ラフソン法により算出される連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、ｊ＋１番目の解の算出の際に、ニュートン・ラフソン反復によって解が非収束であるときに、非収束に対処する対処ステップを実行する非収束判断ステップと、前記ステップ幅があらかじめ定められているステップ幅の下限値より大きい場合には、ステップ幅を小さくしてから、ｊ＋１番目の解を新たに算出し直し、あらかじめ定められているステップ幅の下限値と等しい場合には、既に実行したｊ番目の解を求めるステップをより小さなステップ幅で再実行し、ｊ番目の解を求め直す前記対処ステップと、を含むことを特徴とする連立方程式解法である。
【００４０】第五の本発明は、上記課題を解決するために、連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法であって、追跡中の各解の算出がニュートン・ラフソン法により算出される連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、ｊ＋１番目の解の算出において、ニュートン・ラフソン反復によって解が非収束である時に、ニュートン・ラフソン反復が発散する場合は、第一対処ステップを実行し、ニュートン・ラフソン反復が非発散非収束である場合には、第二対処ステップを実行する非収束判断ステップと、ｊ＋１番目の解の算出における前記ステップ幅があらかじめ定められているステップ幅の下限値より大きい場合は、ステップ幅を小さくしてから、ｊ＋１番目の解を求め直し、前記ステップ幅の下限値と等しい場合には、既に実行したｊ番目の解を求めるステップをより小さなステップ幅で再実行し、ｊ番目の解を求め直す前記第一対処ステップと、枢軸再選択を行ってから、再びニュートン・ラフソン反復を行ってｊ＋１番目の解を求める前記第二対処ステップと、を含むことを特徴とする連立方程式解法である。
【００４１】第六の本発明は連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法であって、追跡中の各解の算出がニュートン・ラフソン法により算出される連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、ｊ＋１番目の解の算出において、ニュートン・ラフソン反復によって解が非収束である時に、ニュートン・ラフソン反復が非発散非収束である場合に対処ステップを実行する非収束判断ステップと、枢軸再選択を行ってから、再びニュートン・ラフソン反復を行ってｊ＋１番目の解を求める前記対処ステップと、を含むことを特徴とする連立方程式解法である。
【００４２】第七の本発明は、上記課題を解決するために、上記第一又は第二の連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、さらに、ｊ番目の解を求める際における、ｊ番目の解からｊ＋１番目の解までのステップ幅が所定のステップ幅下限値以上に設定されるステップ幅設定ステップと、前記ステップ幅設定ステップで求められたステップ幅を用いてｊ＋１番目の解を求める求解ステップと、を含み、前記ステップ幅設定ステップは、固定の下限値と、ｊ番目の解を求める際に用いられたステップ幅の値と、に基づいて前記ｊ＋１番目の解を求める際の所定のステップ幅下限値を算出する下限値設定ステップ、を有することを特徴とする連立方程式解法である。
【００４３】第八の本発明は、上記課題を解決するために、上記第一又は第二の連立方程式解法において、前記ＢＤＦ法における解曲線の追跡における解曲線上のｊ＋１（ｊは整数）番目の解を求めるステップは、さらに、ｊ＋１番目の解を求める際における解の許容誤差を、ｊ＋２番目の解を求める際のステップ幅下限値より小さく設定する許容誤差設定ステップ、を含むことを特徴とする連立方程式解法である。
【００４４】第九の本発明は、上記課題を解決するために、連立方程式をＢＤＦ法を用いて解曲線を追跡することによって解く方法において、前記ＢＤＦ法の修正子方程式の各方程式のうちステップ幅による拘束条件を形成しているステップ幅方程式に対し、ステップ幅の値で両辺を除算するステップと、前記ステップ幅方程式に対し、所定の定数を両辺に乗ずることにより、前記修正子方程式の各方程式のうち前記ステップ幅方程式以外の方程式群に対応するヤコビ行列の絶対値の大きさと、前記ステップ幅方程式に対応するヤコビ行列の絶対値の大きさとを近づける絶対値近似ステップと、を含むことを特徴とする連立方程式解法である。
【００４５】
【作用】以下、本発明の作用について説明する。
【００４６】２. １ 疑似解収束の対策前述の疑似解収束現象は解曲線の方向が急激に変化する箇所で生ずることが多い。つまり、解曲線の方向の変化分が予測に反して急激に大きくなり、これに対してステップ幅が相対的に過大である場合に疑似解収束現象が生ずる。一般に、このような解曲線の方向の急激な変化を事前に予測することは難しい。本文ではこれを実用的な方法で事後評価して対策を講じるステップ輻制御アルゴリズムを提案する。疑似解収束現象を検出するために、次のような予測子と修正子の角度φを評価する。すなわち、ＢＤＦ法の場合のｊ＋１ステップにおける角度φj+1を【数２９】

で表す。ここで、‖・‖はユークリッドノルムを表し、【数３０】

は次のようなｋ階の差分商を表す（上記文献２３参照）。
【００４７】
【数３１】

ただし、【数３２】

とする。このような角度を評価して、その値が所定値【数３３】

より大きいときには疑似解収束現象が発生している可能性があると判定する。この判定は疑似解に収束した直後において、角度が急激に大きくなることが多いという経験に基づいたものである。
【００４８】このような予測子と修正子の角度が大きくなる典型的な例が図１２の（ａ）、（ｂ）〜（ｆ）に示されている。これらは次のような４つのケースに分類することができる。
【００４９】（ａ）真正解に収束しているが、現在のステップ幅が大きい。
【００５０】（ｂ）真正解に収束しているが、直前のステップ幅が大きい（従って、現在のステップ幅を小さくしても角度は小さくならない）。
【００５１】（ｃ）、（ｄ）現在のステップで疑似解に収束している。
【００５２】（ｅ）、（ｆ）直前のステップで疑似解に収束している。
【００５３】上述の角度判定によって、もちろん（ａ）や（ｂ）のような正常なケースも含まれるが、（ｃ）〜（ｆ）のような疑似解収束を検出することができる。従って、このような判定をもとに、解曲線追跡のすべてのステップでこの角度が常に所定値以下になるようにステップ幅を制御すれば、多くの疑似解収東現象を対策できると考えられる。本発明において特徴的な事項であるこの判定法ですべての疑似解収束が完全に検出されるわけではないが、後述するごとく有効性は高い。
【００５４】本発明によるアルゴリズムでは、現在のｊ＋１ステップを実行して、もし角度φj+1 が所定値【数３４】

以内に収まっていれば、直前のｊステップの受け入れを確定して、更に、ｊ＋１ステップを暫定的に受け入れて次のｊ＋２ステップに進む。ｊ＋１ステップで、もし角度が所定値より大きくなった場合、まず、（ａ）、（ｃ）、（ｄ）の可能性があると判定して、そのステップを棄却する。次に、ステツプ幅を適当に小さく取り直してｊ＋１ステップを再実行して、（ａ）、（ｃ）、（ｄ）を対策する。このとき、もしステップ幅を所定値ｈＬまで小さくしても角度が所定値以内に収まらなけれぱ、（ｂ）、（ｅ）、（ｆ）の可能性があると判定して、更に、直前のｊステップを棄却する。そして、ステップ幅を適当に小さく取り直してｊステップを再実行して、（ｂ）、（ｅ）、（ｆ）を対策する。
【００５５】本発明に係るアルゴリズムを用いた解曲線追跡の様子が図１に示されている。図１（ａ）は、現在のステップの予測子と修正子の角度を評価して、場合によっては現在のステップ（破線の小さな球面）に加えていったん暫定的に許容した直前のステップ（破線の大きな球面）も棄却する。そして、ステップ幅を小さく取り直して直前のステップ（実線の小さな球面）を再実行して、常に角度を所定値以内に収めるようにステップ幅を制御しながら解曲線を追跡する。このようにして、図１（ｂ）に示すように、従来のアルゴリズムでは失敗する図９のような解曲線を追跡することができる。
【００５６】疑似解収束の他の対策手段として、平面法に適したステップ幅制御アルゴリズムを球面法に適用することも考えられる。これは隣接するステップの接線ペクトルの角度を所定値以内に収めるようにステップ幅を制御するものである。しかし、このアルゴリズムは接線ペクトルの計算を要することに加え、図２に示すように隣接する接線ペクトル（破線）が平行に近い場合には対策できないという欠点がある。
【００５７】２. ２ ステップ輻の予測解曲線が緩やかに変化しているときには、適当なステップ幅を予測して、前述の角度判定によってステップが棄却されるのを可能な限り少なくするのが効率的である。本節では、前述の角度判定に適したステップ幅予測アルゴリズムに関する発明を提案する。提案する発明のアルゴリズムでは、既に得られた後方値を用いて解曲線を多項式で近似して、その近似解曲線に対して角度判定を満足するステップ幅を見積って、これを次のステッブの予測ステップ幅とする。
【００５８】まず、ｊステップのｓ＝ｓj の近傍ｓP ＝ｓj ＋ｈP で、解曲線がｊステップに至る後方値を用いて局所的に３次の多項式で近似できるものと仮定する。このとき、近似解曲線はニュートンの補間多項式（文献２３参照）を用いて【数３５】

と表すことができる。次に、この近似解曲線にＢＤＦ法を適用すると、ｊステップからｈP だけ進んた点における角度φ（ｈP ）は、式（１１）、（１２）と式（１４）より【数３６】

と表すことができる。従って、角度【数３７】

が与えられたとき、ステップ幅の予測値ｈP は次のようなｈP に関する４次方程式【数３８】

を解いて、その解集合の中から正の最小値を採用すればよい。
【００５９】なお、一度棄却されたｊ＋１ステップを再実行する際には、可能性の大きい図１２の（ａ）の場合を想定して、新たなｊ＋１ステップのステップ幅を予測するのが効率の点で望ましい。つまり、棄却したｊ＋１ステップの解を含むｊ＋１ステップに至る後方値を用いて解曲線を近似し、この近似解曲線をもとに、角度が所定値に収まるような新たなステップ幅ｈR （ｓR ＝ｓj ＋ｈR ）を予測する。このときの角度φ（ｈR ）は【数３９】

で表される。従って、ステップ幅ｈR は【数４０】

を解いて同様に求めることができる。
【００６０】実際のｊ＋１ステップの予測ステップ幅ｈj+1pは、上述の結果をもとに、更に、近似によるアルゴリズムのチャタリングを防ぐために、経験的な定数αを用いて、【数４１】

または、ｊ＋１ステップの再実行の場合には、【数４２】

とする。ここで、α1 ＞１であり、また、0 ＜α2 、α3 …＜１である。更に、ｈmax はステップ幅の上限値を表す。以下、このようなαi はすべて経験的な定数（0 ＜α1 ＜１、ｉ＝２、３、・・・）である。
【００６１】もし暫定的に受け入れたｊステップが棄却された場合は、解曲線が急激に変化している可能性が高く、上述のような近似によるステップ幅の予測は適当といえない。その際の新たなｊステップの予測ステップ幅ｈjPは、暫定的に受け入れた元のｊステップのステップ幅ｈjoldに対して、経験的にｈjp＝α5 ｈjold （２３）
とする。
【００６２】２. ３ ステップ幅の下限値の制御ステップ幅の下限値ｈL を従来のようにすべてある所定値ｈmin に一律に固定するという必要性はない。ＮＲ法の許容誤差を必要以上に厳しくせず、また、無解現象を発生させないために、本節では、ステップ幅の下限値とＮＲ反復の許容誤差を次のように各ステップごとに動的に制御する方法に関する発明を提案する。まず、ｊ＋１ステップのステップ幅の下限値ｈL j+1 は直前のステップ幅に応じて、 ｈL j+1 ＝max （α6 ｈj ，ｈmin ） （２４）
とする。次に、ｊ＋１ステップのＮＲ法の許容誤差ｈL j+1 を次のｊ＋２ステップのステップ幅の下限値に応じて【数４３】

とする。尚、【数４４】

についても同様である。このとき、【数４５】

となるようにα7 を設定すれば、次のｊ＋２ステップで無解現象は生しない。また、解曲線の追跡はその大部分がｈmin より十分大きなステップ幅で進められるので、ＮＲ反復の許容誤差を大幅に緩和することができる。
【００６３】２. ４ ＮＲ反復が非収束の際のステップ幅の制御ＮＲ反復が非収束の場合に、常にステップ幅を小さく取り直してＮＲ反復を再実行するという従来の方法では（文献１２及び１７参照）、ステップ幅を必要以上に小さくして、収束許容誤差も厳しくなり、収束が一層困難になることがある。このような問題を対策するために、本節では、ＮＲ反復が非収束となる場合を次のように２つの場合に分割してそれぞれの場合に応じてステップ幅を制御する方法の発明を提案する。
【００６４】（１）まず、ＮＲ反復の途中て発散となった場合およびヤコビ行列の行列式の符号が反転した場合は、従来と同様に、直ちに反復を停止して、ステップ幅をｈj+1 ＝max （α3 ｈj+1 ，ｈL j+1 ） （２６）
に取り直して、ＮＲ反復を再実行する。ステップ幅が下限値ｈL j+1 である場合においても発散するときは、前述と同様に直前のｊステップ（の結果）を放棄する。
【００６５】（２）次に、ＮＲ反復が非発散非収束の場合、同じステップ幅で枢軸選択を再実行した後、ＮＲ反復を再実行する。これは、ステップ幅を小さく取り直しても非発散非収束は対策できず、むしろ、枢軸再選択が有効であるという経験に基づいたものである。枢軸再選択が有効なのは、これによって数値的な悪条件の状態を改善できるためと考えられる。
【００６６】尚、これまでに述べた２. １〜２. ４までの発明を応用した解曲線追跡方法のフローチャートが図１４に示されている。
【００６７】図１４に示されているように、Ｓ１４−１においては、ｊやＳ0 などの値の初期値の設定が行われる。Ｓ１４−２においては、予測子ｙj+1Pの予測が行われる。Ｓ１４−３においては、ニュートン・ラフソン法による反復が行われる。そしてＳ１４−４においては、ニュートン・ラフソン法による反復が収束するか否かについて判定がなされる。
【００６８】この判定の結果、非発散非収束である場合には、Ｓ１４−５に処理が移行し、枢軸再選択が行われ、その後Ｓ１４−３において再びニュートン・ラフソン法による反復が行われる。
【００６９】一方、上記判定の結果、発散であると判定された場合には、Ｓ１４−６において現在のステップ（ｊ＋１ステップ）のステップ幅が、ステップ幅の下限値より大きいか否かが判断される。下限値より大きい場合には、Ｓ１４−７においてステップ幅が小さくされた後、Ｓ１４−２から再びニュートン・ラフソン法の反復が行われる。逆に、現在のステップ幅が下限値と等しい場合には、Ｓ１４−８において直前のステップ（ｊステップ）の予測ステップ幅が小さく設定された後、Ｓ１４−９において、ｊを一つ前に戻して（デクリメントし）、一つ前の解を求めることからやり直す。
【００７０】また、上記判定の結果、収束であると判定された場合には、Ｓ１４−１０において予測子と修正子とのなす角度の角度判定が行われる。この判定の結果、角度が所定の基準角度より大きい場合（すなわちＮＧの場合）には、Ｓ１４−１１において、現在のステップ幅が、ステップ幅の下限値より大きいか否かが判断される。下限値より大きい場合は、Ｓ１４−１２においてステップ幅が小さく設定された後、この予測ステップ幅をｊ＋１番目の解を求める際の実際のステップ幅として（Ｓ１４−１３）から、再び、ニュートン・ラフソン法の反復がＳ１４−２から行われる。
【００７１】逆に、現在のステップ幅が下限値と等しい場合には、Ｓ１４−８において直前のステップ幅が小さく設定された後、Ｓ１４−９において、ｊを一つ前に戻して（デクリメントし）、一つ前の解を求めることからやり直す。
【００７２】尚、上記角度の判定の結果、所定の基準角度より小さいと判断された場合には、Ｓ１４−１４において終了判断が行われる。Ｓ１４−１４において終了であると判断された場合には、最終的な解が得られたことになり、計算が終了するが、まだ終了ではないと判断された場合には、Ｓ１４−１５においてｊがインクリメントされた後、Ｓ１４−１６において、予測ステップ幅が算出され、Ｓ１４−１７においてステップ幅の下限値が更新され、Ｓ１４−１８において、ニュートン・ラフソン法における許容誤差が設定され、Ｓ１４−１９においてニュートン・ラフソン法における別の許容誤差も設定される。しかる後に、Ｓ１４−１３において、予測ステップ幅が実際に用いられるステップ幅として採用され、次の解を求める処理へ移行する。
【００７３】２. ５ 修正子方程式の尺度変換上記式（８）のヤコビ行列のｎ＋１行目の要素の絶対値は、ヤコビ行列の他の行の大部分の要素の絶対値（バイポーラアナログICでは、経験的にl0-2〜10-6のオーダ、または、それ以下）に比べて非常に大きく、更に、ステップ幅に依存して大きく変化する。このことは、ＢＤＦ法では【数４６】

となることから明らかである。この結果、枢軸選択がステップ幅の影響を受けて安定に行われず、場合によっては充填要素が多発する。これを対策するたあに、本節では、次のような修正子方程式の尺度変換を提案する。
【００７４】まず、修正子方程式として、式（５ｂ）の代わりにα8 ／ｈj+1 を掛けた次の式【数４７】

を用いる。このとき、ヤコビ行列のｎ＋１行目は【数４８】

となって、ステップ幅にほとんど依存しないほほ一定値に保つことができる。また、α8 を適当な値に設定すれば、ヤコビ行列の他の行の要素の絶対値のオーダに調整できる。このように尺度変換することにより、汎用回路シミュレータで広く用いられているＭａｒｋｏｗｉｔｚの方法を用いて枢軸選択を安定に行うことができる。
【００７５】要約すれば、以下のようになる。
【００７６】すなわち、第一の本発明の判定ステップは、予測子と修正子とのなす角度が所定値より大きければ追跡に異常が生じたと判断するのである。そして、異常が生じた場合には、最近に求められた解を放棄し、一つ前に求められた解から、再び解の追跡をやり直すのである。これによって、追跡が失敗したことが高い確率で判定可能である。
【００７７】第二の本発明によれば、予測子と修正子とのなす角度が上記式により容易に求めることが可能となる。
【００７８】第三の本発明によれば、予測子と修正子とのなす角度を予測して、この角度が所定角度以下になるように、予測ステップ幅を算出されるので、追跡の失敗を未然に防止できる確立が高くなる。
【００７９】第四の本発明によれば、ニュートン・ラフソン反復が非収束であって、かつステップ幅が下限値の時は、既に求めた一つ前の解から新たに求め直す。そのため、非収束の場合でも、より高い確率で収束状態とすることが可能である。
【００８０】第五の本発明によれば、ニュートン・ラフソン反復が非収束の場合を、発散の場合と、非発散非収束の場合に分けてそれぞれの場合に応じて対処した。従って、より的確な対処が可能となる。
【００８１】第六の本発明によれば、一定の場合に枢軸再選択を行うので、数値的悪条件を改善することができる。
【００８２】第七の本発明によれば、ステップ幅の下限値を、直前に解が求められた際のステップ幅の値に応じて変化させたので、より効率のよい解曲線の追跡が行える。第八の本発明によれば、解の許容誤差を、直前に解が求められた際のステップ幅の値に応じて変化させたので、無解となってしまうことを未然に防止可能である。
【００８３】第九の本発明によれば、ステップ幅の値で、ステップ幅方程式の方程式の両辺を除算したので、ヤコビ行列の要素の値をステップ幅の値とは独立な一定値に設定することができる。
【００８４】
【実施例】本文で提案したアルゴリズムを汎用回路シミュレー夕に実装し、アナログＬＳＩを含むいくつかの回路の解折に適用した例を示してその有効性を確認する。使用する汎用回路シミュレータ三洋ＳＰＩＣＥには既にＢＤＦ法（球面法の一種である斜円柱の方法）が実装されているので（文献１７を参照）、これをもとにしてステップ幅制御アルゴリズムを実装した。また、コンピュータはＩＢＭ ＲＳ６０００−５９０を使用した。
【００８５】以下の数値例において、アルゴリズムの各パラメータは経験的に次のような値を用いた。
【００８６】
ｈmin ＝１０-4 ｈmax ＝１０ ｈ1 ＝１０-2【数４９】

α1 ＝２ α2 ＝０. ５ α3 ＝０. ３３ α4 ＝０. ８α5 ＝０. ２５ α6 ＝０. １２５ α7 ＝０. １ α8 ＝１０-4（１）トンネルダイオード直列回路文献（１２）、（１３）に記載されているトンネルダイオード直列回路に本発明を適用したときの解曲線追跡の様子が図３に示されている。図３は、各ステップの近似解をプロットしたもので、この中には暫定的に受け入れられて角度判定で棄却されたものも含む（以下同様）。図３に示されているグラフは、縦軸が電圧Ｖ２であり、横軸は印加電圧Ｖｉｎが示されている。解曲線の追跡に要した総ステップ数は６８（文献（１２）では６８）であり、受け入れステップ数は６３（同４７）、棄却ステップ数は５（同２１）、また、ＮＲ反復回数の総数は３１７（同４２８）である。解曲線の方向の変化の緩急に応してステップ幅が効果的に変化している様子が図３から理解されよう。図３に示されているように、解曲線が比較的緩やかに変化する場合には、棄却ステップ数が少なく、ステップ幅予測アルゴリズムが有効に作動している。
【００８７】（２）定電圧電源用ＩＣ低電圧電源用ＩＣとして、総素子数３７８（ＢＪＴ １８１）のバイポーラアナログＩＣを用いた。図４にはこのＩＣの回路の直流動作点解折の際の解曲線追跡の様子をある節点電圧の変化を、２次元平面上にプロットしたグラフが示されている。図４に示されているグラフにおいて、縦軸は電圧Ｖ２であり、横軸は電圧Ｖ１である。また、図４の原点付近の長方形で囲んだ部分を拡大した拡大図が図５である。図５においても縦軸は電圧Ｖ２であり、横軸は電圧Ｖ１である。解曲線は少なくとも二つの要素からなっており、迫跡の途中でいったん他の解曲線（破線で示したループ）に乗り移った（暫定的に受け入れられた）直後に、角度判定で棄却されて元の解曲線に戻っている様子が理解されよう。また、実用的なＩＣでもこのようなループを含む解曲線をもつ回路が実在する典型的な例と言えよう。本願発明者らがこのような回路に遭遇したことは、本発明をなすための契機となったものである。
【００８８】（３）大規模回路の直流動作点解析最後の例は１０，０００素子級の映像信号処理用バイポーラアナログＬＳＩである。これらの回路の直流動作点解析の結果が図６に示されている。図６には、ＡからＥまで５つのＬＳＩの例が示されている。従来のアルゴリズムでは、疑似解収束や充填要素の多発、更に、非発散非収束などによって、いずれの回路も解析不可能であった。また、従来のアルゴリズムに方程式の尺度変換だけを適用した場合、ステップ幅制御のパラメータを各回路ごとに微調整することにより、回路Ａ、Ｂ、Ｅに関しては解析可能とすることができた。しかし、その他の回路Ｃ、Ｄはそのような工夫した従来の手段を採用しても解析不可能であった。本文で提案した発明のアルゴリズムでは、回路ごとにステップ幅制御のパラメータを変更することなしに、いずれも実用的な時間で解析可能である。以上のように、本文で提案するステップ幅制御アルゴリズムは、球面法における疑似解収束やＮＲ反復の非収束などの対策に有効で、実用的であることが理解されよう。
【００８９】以上のように、本実施例においては、解曲線追跡のための確実性の高い実用的なステップ幅制御アルゴリズムを提案し、大規模回路に適用してその有効生が確認された。すなわち、本発明によれば、効果的な解曲線の追跡を行うことができ、結果として連立方程式を効率的に解くことが可能となる。
【００９０】さらに、本発明のアルゴリズムに加えて次のような既知の技法を併用すると、より効果的に解曲線を追跡できることが経験的に判明している。すなわち、変数ｔの尺度変換（文献１３参照）やホモトピー方程式の摂動（perturbation）（文献３及び４参照）は解曲線を滑らかにする上で効果がある。また、高次の予測子はＮＲ反復数を滅少させる上で効果がある。提案するアルゴリズムはl 次のＢＤＦ解曲線追跡アルゴリズムを想定しているが、ステップ幅の予測で既に解曲線を高次（３次）の多項式で近似していることから、予測子も高次にするのは容易である。
【００９１】
【発明の効果】以上述べたように、第一の本発明によれば、角度によって、追跡が失敗したか否かを確認するので、追跡の正否をより高い確率で判断可能な解曲線の追跡が可能となる。
【００９２】第二の本発明によれば、角度を簡易な方法で算出可能であるので、より高速に追跡の正否の判断が可能となる。
【００９３】第三の本発明によれば、あらかじめ予測により、追跡が失敗しないようにステップ幅を設定したので、失敗を未然に防止し、より正確でかつ迅速な解曲線の追跡が可能となる。
【００９４】第四の本発明によれば、ニュートン・ラフソン反復において非収束となった場合でもより高い確率で収束状態に移行できるので、より多くの場合に対して解曲線の追跡が可能となる。
【００９５】第五の本発明によれば、非発散非収束の場合には枢軸再選択を行うので、より的確に非収束状態を解消できる。そのため、より多くの場合に対して解曲線の追跡が可能となる。
【００９６】第六の本発明によれば、一定の場合にのみ、枢軸再選択を行ったので、数値的悪条件を脱し、非収束状態を解消し、解曲線の追跡が可能となる。
【００９７】第七の本発明によれば、ステップ幅の下限値を解が求められたときのステップ幅の値に応じて変化させたので、より効率の良い、解曲線の追跡が可能となる。第八の本発明によれば、解を求める際の許容誤差を適応的に変化させたので、無解となることを未然に防止し、より多くの場合に対して解曲線の追跡が可能となる。
【００９８】第九の本発明によれば、修正子方程式の各方程式に対応するヤコビ行列の絶対値の大きさを近づけることが可能となるので、枢軸選択を安定に行うことができて、いわゆる数値的な悪条件を回避することが可能となり、より正確でかつ迅速な解曲線の追跡が可能となる。
【出願人】 【識別番号】０００００１８８９
【氏名又は名称】三洋電機株式会社
【出願日】 平成６年（１９９４）１１月２５日
【代理人】 【弁理士】
【氏名又は名称】吉田 研二 （外２名）
【公開番号】 特開平８−１４７２６７
【公開日】 平成８年（１９９６）６月７日
【出願番号】 特願平６−２９０９９１