(57)【要約】【課題】ＣＰＥを含む電気回路を、任意の回路構成及び任意の入力電圧に対して、時間領域で解くことのできる技術を提供する。【解決手段】コンピュータが、境界条件として少なくとも前記電気回路に印加される入力電圧を取得する取得ステップと、前記境界条件に基づいてＣＰＥに流れる枝電流を計算する第１の計算ステップと、前記境界条件に基づいてＣＰＥ以外の回路素子に流れる枝電流を計算する第２の計算ステップと、前記境界条件と前記第１及び第２の計算ステップでの計算結果とから、キルヒホッフの電流則及びキルヒホッフの電圧則を用いて時間領域で未知変数を求めるステップと、を実行する。【選択図】図１
2011145944.tif 000002
【特許請求の範囲】
【請求項１】
回路素子としてＣＰＥ（Constant Phase Element）を含む電気回路を時間領域で解くための数値計算方法であって、
コンピュータが、境界条件として少なくとも前記電気回路に印加される入力電圧を取得する取得ステップと、
コンピュータが、前記境界条件に基づいてＣＰＥに流れる枝電流を計算する第１の計算ステップと、
コンピュータが、前記境界条件に基づいてＣＰＥ以外の回路素子に流れる枝電流を計算する第２の計算ステップと、
コンピュータが、前記境界条件と前記第１及び第２の計算ステップでの計算結果とから、キルヒホッフの電流則及びキルヒホッフの電圧則を用いて時間領域で未知変数を求めるステップと、を有することを特徴とする数値計算方法。
【請求項２】
請求項１に記載の数値計算方法において、
ＣＰＥが、電位がＶ０である端子０と電位がＶ１である端子１との間に存在し、物性値α及びＴを持つ場合に、
前記第１の計算ステップにおいて、ａを計算開始時刻とし、端子０から端子１に流れる向きを正として、時刻ｔにおけるＣＰＥに流れる枝電流Ｉを分数階微分ａＤｔαを用いて、
Ｉ＝ａＤｔα{Ｔ（Ｖ０−Ｖ１）}
と計算することを特徴とする数値計算方法。
【請求項３】
請求項１に記載の数値計算方法において、
ＣＰＥが、電位がＶ０である端子０と電位がＶ１である端子１との間に存在し、物性値α及びＴを持つ場合に、
前記第１の計算ステップにおいて、ａを計算開始時刻とし、端子０から端子１に流れる向きを正として、時刻ｔにおけるＣＰＥに流れる枝電流Ｉを分数階微分ａＤｔαを用いて、
Ｉ＝ＴａＤｔα（Ｖ０−Ｖ１）
と計算することを特徴とする数値計算方法。
【請求項４】
請求項１〜３のうちいずれか１項に記載の数値計算方法の各ステップをコンピュータに実行させることを特徴とするプログラム。
【請求項５】
請求項４に記載のプログラムを記録したことを特徴とするコンピュータ読み取り可能な記録媒体。
【請求項６】
回路素子としてＣＰＥ（Constant Phase Element）を含む電気回路を時間領域で解くための数値計算装置であって、
境界条件として少なくとも前記電気回路に印加される入力電圧を取得する取得手段と、
前記境界条件に基づいてＣＰＥに流れる枝電流を計算する第１の計算手段と、
前記境界条件に基づいてＣＰＥ以外の回路素子に流れる枝電流を計算する第２の計算手段と、
前記境界条件と前記第１及び第２の計算手段での計算結果とから、キルヒホッフの電流則及びキルヒホッフの電圧則を用いて時間領域で未知変数を求める手段と、
を有することを特徴とする数値計算装置。
【請求項７】
請求項６に記載の数値計算装置において、
ＣＰＥが、電位がＶ０である端子０と電位がＶ１である端子１との間に存在し、物性値
α及びＴを持つ場合に、
前記第１の計算手段が、ａを計算開始時刻とし、端子０から端子１に流れる向きを正として、時刻ｔにおけるＣＰＥに流れる枝電流Ｉを分数階微分ａＤｔαを用いて、
Ｉ＝ａＤｔα{Ｔ（Ｖ０−Ｖ１）}
と計算することを特徴とする数値計算装置。
【請求項８】
請求項６に記載の数値計算装置において、
ＣＰＥが、電位がＶ０である端子０と電位がＶ１である端子１との間に存在し、物性値α及びＴを持つ場合に、
前記第１の計算手段が、ａを計算開始時刻とし、端子０から端子１に流れる向きを正として、時刻ｔにおけるＣＰＥに流れる枝電流Ｉを分数階微分ａＤｔαを用いて、
Ｉ＝ＴａＤｔα（Ｖ０−Ｖ１）
と計算することを特徴とする数値計算装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明はConstant Phase Element（ＣＰＥ）を含む電気回路を時間領域で解くための数値計算方法、数値計算装置、プログラム、及び、プログラムを記憶した記録媒体に関するものである。
【背景技術】
【０００２】
近年では、電気製品のような電気で駆動又は制御される装置を設計する際に、装置に組み込む回路・モジュール・材料等の電気的な挙動を数値計算（シミュレーション）によって予測することにより、開発の効率を向上させることが行われている。
数値計算を行うためには、まず、解析対象となる物の電気特性を正確に把握する必要がある。電気特性を調査する方法として、例えば、交流インピーダンス測定が知られている。交流インピーダンス測定とは、周波数を変えながら交流電圧を印加し解析対象のインピーダンスを測定する方法であり、その測定結果から等価回路の同定と回路素子の物性値の決定を行うことができる。
【０００３】
交流インピーダンス測定では等価回路同定のためにいくつもの回路素子を用いてフィッティングを行うが、その中のよく知られた回路素子としてConstant Phase Element（ＣＰＥ）がある。ＣＰＥは抵抗と容量の中間的性質を持つ回路素子であり、そのインピーダンスは次の式（１）で表わされる。
ＺＣＰＥ＝１／｛（ｉω）αＴ｝    (1)
ここでＺＣＰＥはＣＰＥのインピーダンスであり、ｉは虚数単位であり、Ｔ及びαはフィッティングパラメータでかつＣＰＥで表現される物質の物性値である。
【０００４】
数値計算の対象となる等価回路がＣＰＥを含むものである場合には、以下に述べるような課題が生じる。
【０００５】
入力波形を周波数展開して数値計算する時にはＣＰＥを式（１）のままの形で解くことができる。しかし周波数展開が可能なのはあくまでも線形の範囲内であり、実際の物質では物性値はしばしば電界依存性など非線形な振る舞いを示すため適用できない。非線形な振る舞いを示す物質を数値計算する際、一般的には時間発展の時間刻み幅を小さくすることで、区分的には線形に振る舞うとして数値計算を行う。そのため時間領域の方程式に書き換える必要がある。
また電気回路などでは入力電圧として矩形波を用いることが多く、周波数展開すると無限大の周波数成分を必要とし、計算コストが非常に多い。よって任意の入力波形に対応するためには周波数ではなく時間領域での方程式に書き換える必要がある。
【０００６】
非特許文献１によれば、ＣＰＥに流れる電流Ｉは分数階微分を用いて式（２）のように書くことができる。
【数１】
2011145944.tif 000003
ここでｄα／ｄｔαは分数階微分を示す記号であり、ΔＶはＣＰＥの両端にかかる電位差である。また、非特許文献１では、ＣＰＥもしくは抵抗が２つ直列に並んだ回路に対して
特定の種類の入力電圧を印加した場合に、回路素子間の電位がMittag-Leffler’s functionを用いて書くことができると述べている。しかしながら、非特許文献１で述べられている数値計算方法は、特定の回路構成と特定の入力電圧に限定されており、この方法を任意の回路構成、任意の入力電圧の場合に汎用的に適用することは難しい。
【０００７】
また実際の物質では物性値はしばしば電界依存性を持つ。即ち物性値Ｔは電界の関数であり、電界が時間に伴い変化する場合、物性値Ｔは時間に依存する。
さて普通のコンデンサでは静電容量Ｃは式（３）のように定義される。
【数２】
2011145944.tif 000004
回路の方程式を導く場合、この関係はいかなる時刻によっても成り立つと考える。そのため電流は式（４）のようになる。
【数３】
2011145944.tif 000005
Ｃが時間に依存しないならば式（５）のようになる。
【数４】
2011145944.tif 000006
以上の事柄とαが１の時にＣＰＥはコンデンサと一致する必要があることとの両者から推察するとＣＰＥの電流の表式（２）は物性値Ｔが時間に依存するときには適切ではなく、このままでは物質の電流電圧特性を正確に再現することができない。
【０００８】
なお、分数階微分に関する参考文献としては、非特許文献２〜４があげられる。非特許文献２では均質な粘弾性体を対象として時間領域で分数階微分の方程式を解いている。また非特許文献３では分数階微分の定義等が示されており、非特許文献４では分数階微分の数値計算をする時の離散化手法について幾つか紹介されている。また、一般的な電気回路の数値計算（シミュレーション）については、非特許文献５に詳しい。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】S. Westerlund, “Capacitor Theory”, IEEE Transactions on Dielectrics and Electrical Insulation 1 5 (1994) 826-839.
【非特許文献２】L. B. Eldred, W. P. Baker, A. N. Palazotto, “Numerical Application of Fractional Derivative Model Constitutive Relations for Viscoelastic Matericals”, Computer & Structures 60 6 (1996) 875-882.
【非特許文献３】I. Podlubny, “Fractional Differential Equations”, Academic Press, (1999) ISBN0-12-558840-2.
【非特許文献４】K. B. Oldham, J. Spanier, “The Fractional Calculus”, Dover Publications Inc. (2006) ISBN0-486-45001-5.
【非特許文献５】浅井秀樹、渡部貴之, 「電子回路シミュレーション技法」, 科学技術出版 (2002) ISBN4-87653-200-1.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
本発明は上記実情に鑑みてなされたものであって、その目的とするところは、ＣＰＥを
含む電気回路を、任意の回路構成及び任意の入力電圧に対して、時間領域で解くことのできる技術を提供することにある。
【課題を解決するための手段】
【００１１】
本発明の第１態様は、回路素子としてＣＰＥ（Constant Phase Element）を含む電気回路を時間領域で解くための数値計算方法であって、コンピュータが、境界条件として少なくとも前記電気回路に印加される入力電圧を取得する取得ステップと、コンピュータが、前記境界条件に基づいてＣＰＥに流れる枝電流を計算する第１の計算ステップと、コンピュータが、前記境界条件に基づいてＣＰＥ以外の回路素子に流れる枝電流を計算する第２の計算ステップと、コンピュータが、前記境界条件と前記第１及び第２の計算ステップでの計算結果とから、キルヒホッフの電流則及びキルヒホッフの電圧則を用いて時間領域で未知変数を求めるステップと、を有することを特徴とする数値計算方法である。
【００１２】
本発明の第２態様は、回路素子としてＣＰＥ（Constant Phase Element）を含む電気回路を時間領域で解くための数値計算装置であって、境界条件として少なくとも前記電気回路に印加される入力電圧を取得する取得手段と、前記境界条件に基づいてＣＰＥに流れる枝電流を計算する第１の計算手段と、前記境界条件に基づいてＣＰＥ以外の回路素子に流れる枝電流を計算する第２の計算手段と、前記境界条件と前記第１及び第２の計算手段での計算結果とから、キルヒホッフの電流則及びキルヒホッフの電圧則を用いて時間領域で未知変数を求める手段と、を有することを特徴とする数値計算装置である。
【００１３】
本発明の第３態様は、上記数値計算方法の各ステップをコンピュータに実行させることを特徴とするプログラムである。
【００１４】
本発明の第４態様は、上記プログラムを記録したことを特徴とするコンピュータ読み取り可能な記録媒体である。
【発明の効果】
【００１５】
本発明によれば、ＣＰＥを含む電気回路を、任意の回路構成及び任意の入力電圧に対して、時間領域で解くことが可能である。
【図面の簡単な説明】
【００１６】
【図１】本発明の実施形態に係る数値計算方法のフローチャート。
【図２】本発明の実施形態に係る数値計算装置の構成例。
【図３】実施例１、５及び６で計算対象とする回路。
【図４】実施例１での計算結果。
【図５】実施例２で計算対象とする回路。
【図６】実施例２での計算結果。
【図７】実施例３で計算対象とする回路。
【図８】実施例３での計算結果（α＝１の時のものとα＝０．５の時のもの）。
【図９】実施例４で計算対象とする回路。
【図１０】実施例４での計算結果（α＝１の時のものとα＝０．５の時のもの）。
【図１１】実施例５で計算対象とする回路。
【図１２】実施例５及び６での印加電圧を示すグラフ。
【図１３】実施例５での計算結果。
【図１４】実施例６での計算結果（Ｔが時間依存しない場合と時間依存する場合）。
【発明を実施するための形態】
【００１７】
本発明の実施形態は、回路素子としてＣＰＥ（Constant Phase Element）を含む電気回路を時間領域で解くための数値計算方法を提供するものである。回路素子としては、ＣＰ
Ｅ以外に、抵抗、コンデンサ、コイルなどがあるが、本方法は任意の回路素子を含む任意の回路構成に適用することができる。このとき、各種物性値（ＣＰＥの物性値α及びＴ、抵抗値、キャパシタンス、インダクタンスなど）が回路素子ごとに異なっていてもかまわない。また、本方法は、解析対象自体が電気回路であることに限定されない。物の電気的な振る舞いが等価な電気回路（等価回路）で近似できる場合には、本方法は、そのような等価回路の数値計算にも好ましく適用できる。
【００１８】
図１に、本実施形態に係る数値計算方法のフローチャートを示す。本方法は、概略、計算に必要なデータを取得するステップ（Ｓ１００）、ＣＰＥに流れる枝電流を計算するステップ（Ｓ１０１）、ＣＰＥ以外の回路素子に流れる枝電流を計算するステップ（Ｓ１０２）、未知変数を求めるステップ（Ｓ１０３）からなる。後述するが、これらのステップはコンピュータにより実行されるものである。以下、各ステップについて説明する。
【００１９】
（１）ステップＳ１００（取得ステップ）
ステップＳ１００では、数値計算に必要な各種のデータが取得される。例えば、回路素子の物性値や接続関係などの回路構成、電気回路に印加する入力電圧などの境界条件、計算終了時刻や時間刻みなどの計算条件、などのデータが取得される。なお、これらのデータをユーザに都度入力させてもよいし、記憶媒体に記憶されているデータを読み込んでもよい。
【００２０】
（２）ステップＳ１０１（第１の計算ステップ）
ステップＳ１０１では、ステップＳ１００で得られた境界条件に基づいてＣＰＥに流れる枝電流を計算する。
【００２１】
前述のように、ＣＰＥのインピーダンスＺＣＰＥは下記式（１）で表わされる。
ＺＣＰＥ＝１／｛（ｉω）αＴ｝    (1)
ここで、α及びＴはＣＰＥの物性値であり、ｉは虚数単位であり、ωは角振動数である。但し当然のことながらこの式は線形の範囲内の表記である。例えばＴが電界依存性を持つ場合には、線形性が保てる程度の微小交流を用いなければ式（１）を用いることができない。Ｔの電界依存性の値を求めるためには直流バイアスをかける必要がある。
上記特性を持つＣＰＥが、電位がＶ０である端子０と電位がＶ１である端子１との間に存在するとした場合、時刻ｔにおけるＣＰＥの枝電流Ｉは、分数階微分ａＤｔαを用いて下記式（６）もしくは下記式（７）のように求める。どちらの式を用いたかは各実施例で異なるため、各実施例で述べる。
Ｉ＝ＴａＤｔα（Ｖ０−Ｖ１）    (6)
Ｉ＝ａＤｔα{Ｔ（Ｖ０−Ｖ１）}    (7)
ここで、ａは計算開始時刻であり、端子０から端子１に流れる向きを正とする。
【００２２】
分数階微分の定義には幾つか種類があり、ここではGrunwald-Letnikovの定義を元にし
て非特許文献４に示されるＧ１アルゴリズムを用いて計算することとした。詳細を以下に示す。
Grunwald-Letnikovの分数階微分の定義では関数ｆに作用するα階の微分を式（８）と
定義する。
【数５】
2011145944.tif 000007
但しαは実数、ａは計算開始時刻、ｔは現在の時間、ｍΔｔ＝ｔ−ａであり、ｍは整数で、Δｔは時間刻み幅である。また式(９）の定義を用いている。
【数６】
2011145944.tif 000008
上記の定義を元に離散化を行い、数値計算を行った。具体的にはGrunwald-Letnikovの
分数階微分の定義で極限をとらずにそのまま計算した。つまり以下の実施例では式（１０）のように数値計算を行った。
【数７】
2011145944.tif 000009
【００２３】
このように本実施形態では、ＣＰＥの枝電流を分数階微分を用いて時間領域で計算することとしたので、任意の波形の入力電圧に対し、少ない計算コストで精度良い計算結果を得ることができる。なお、分数階微分の離散化手法には上記したものの他にも様々な計算手法が考えられる（例えば非特許文献４参照）が、ステップＳ１０１の計算にはいずれの手法を用いてもよい。
【００２４】
（３）ステップＳ１０２（第２の計算ステップ）
ステップＳ１０２では、ステップＳ１００で得られた境界条件に基づいてＣＰＥ以外の回路素子に流れる枝電流を計算する。ＣＰＥ以外の回路素子（抵抗、コンデンサ、コイルなど）に流れる枝電流の計算方法は公知であり（例えば非特許文献５参照）、ステップＳ１０２ではそれらの計算方法を利用すればよい。
以下、抵抗とコンデンサに流れる枝電流の計算方法を例示する。
抵抗値Ｒを持つ抵抗が、電位がＶ０である端子０と電位がＶ１である端子１との間に存在するとした場合、時刻ｔにおける抵抗に流れる枝電流Ｉは、下記式（１１）のように求まる。
【数８】
2011145944.tif 000010
但し端子０から端子１に流れる向きを正とする。
静電容量Ｃを持つコンデンサが、電位がＶ０である端子０と電位がＶ１である端子１との間に存在するとした場合、時刻ｔにおけるコンデンサに流れる枝電流Ｉは、下記式（１２）のように求まる。静電容量が時間依存性を持つことを想定しているので、静電容量は時間の関数である。
【数９】
2011145944.tif 000011
但し端子０から端子１に流れる向きを正とする。
【００２５】
（４）ステップＳ１０３
ステップＳ１０３では、ステップＳ１００で得られた境界条件とステップＳ１０１及び
Ｓ１０２での計算結果とから、キルヒホッフの電流則及びキルヒホッフの電圧則を用いて時間領域で未知変数を求める。回路構成や数値計算の目的などにより未知変数は変わるが、例えば、回路の各節点における電位、回路全体に流れる枝電流などが未知変数となる。
【００２６】
（装置構成）
図２に、上述した数値計算を実行する数値計算装置の構成例を示す。数値計算装置は、ＣＰＵ２０２、入力装置２０３、一次記憶装置２０４、出力装置２０５、通信装置２０６、二次記憶装置２０７などを備えるコンピュータにより構成される。ＣＰＵ２０２は、演算及び制御を担う処理装置である。入力装置２０３は、マウス、キーボード、スキャナなどの情報入力装置である。ステップＳ１００の境界条件や計算条件の入力などに利用される。一次記憶装置２０４は、メモリ（ＲＡＭ）などの一時的な記憶装置である。出力装置２０５は、ディスプレイ、プリンタなどの情報出力装置である。数値計算の結果の表示やプリントアウトに利用される。通信装置２０６は、ＬＡＮボードなどであり、他のコンピュータやネットワーク機器との間で有線又は無線による通信を可能にする装置である。二次記憶装置２０７は、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭなどのデータ保存のための装置を指す。符号２０８はバスラインであり、コンピュータ内の装置群間の情報を伝達する役割を持つ。符号２０９はＬＡＮケーブルであり、異なるコンピュータ間の通信装置に接続し、情報を伝達する。
【００２７】
上述した数値計算方法を実現するためのプログラムは、二次記憶装置群２０７に保存されている。当該プログラムが一次記憶装置群２０４に読み込まれ、ＣＰＵ２０２で実行されることにより、上述したステップＳ１００〜Ｓ１０３の処理が実現される。なお、このプログラムを利用者に提供する場合には、フラッシュメモリ、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭなどのコンピュータ読み取り可能な記録媒体にプログラムを記録して提供してもよいし、ネットワークを通じてプログラムを提供してもよい。
【００２８】
［実施例１］
実施例１として、図３に示すような抵抗とＣＰＥが並列に並んだ回路の数値計算例を示す。本実施例では式（６）を用いた。端子３０１にＶＡsin(２πｆｔ)［Ｖ］なる電位を
かけるものとする。ＶＡ＝１００とする。ｆは振動数であり、本実施例ではCole-Cole Plot（コールコールプロット）を作成するために約３０［Ｈｚ］から１×１０８［Ｈｚ］の範囲で振って計算するものとする。Cole-Cole Plotとは、周波数毎のインピーダンスの実数成分と虚数成分を複素平面にプロットしたグラフである。端子３０２の電位は０［Ｖ］で固定とする。符号３０３は抵抗を示し、抵抗値Ｒ＝１［Ω］とする。符号３０４はＣＰＥを示し、Ｔ＝１．５９×１０−３[Ｆ／ｓｅｃ０．５]、α＝０．５とする。電流は端子３０１から端子３０２へ流れる向きを正とする。
【００２９】
本実施例の場合には全ての電位は境界条件として与えられているため、ステップＳ１０３では回路の枝電流（つまり、抵抗３０３に流れる枝電流とＣＰＥ３０４に流れる枝電流の和）を求めるだけでよい。
【００３０】
上述した数値計算方法を用いて図３に示す回路の枝電流を計算し、インピーダンスを求め、Cole-Cole Plotを求めたものを図４に示す。Cole-Cole Plotの計算に際して、電流は３０３及び３０４を流れる枝電流の和を取った。本方法の計算結果は図４のプロット（黒丸）で示される。
ここで、計算結果の検証のために、インピーダンスの理論式による予測を行った。抵抗とＣＰＥが並列に並んだ回路のインピーダンスは、下記式（１３）のように書くことができるので、これを用いて理論値を計算した。
【数１０】
2011145944.tif 000012
図４に示すように理論予測と計算結果はよい一致をしており、本方法を用いて図３に示す回路が適切に計算できていることがわかる。
【００３１】
［実施例２］
実施例２として、図５に示すような図３の回路が２つ直列に並んだ回路の数値計算例を示す。本実施例では式（６）を用いた。端子５０１にＶＡsin(２πｆｔ)［Ｖ］なる電位
をかけるものとする。ＶＡ＝１００とする。ｆは振動数であり、本実施例ではCole-Cole Plotを作成するために約１［Ｈｚ］から１×１０８［Ｈｚ］の範囲で振って計算するものとする。端子５０３の電位は０［Ｖ］で固定とする。符号５０４、５０６は双方とも抵抗を示し、双方とも抵抗値Ｒ＝１［Ω］とする。符号５０５、５０７は双方ともＣＰＥを示し、ＣＰＥ５０５はＴ＝１．５９×１０−３[Ｆ]、α＝１．０とし、ＣＰＥ５０７はＴ＝１．５９×１０−３[Ｆ／ｓｅｃ０．５]、α＝０．５とする。電流は端子５０１から端子５０３へ流れる向きを正とする。
ステップＳ１０３ではまず節点５０２の電位を変数として求め、その後に回路全体の枝電流を求める。
【００３２】
上述した数値計算方法を用いて図５に示す回路の枝電流を計算し、インピーダンスを求め、Cole-Cole Plotを求めたものを図６に示す。電圧は端子５０１と端子５０３の差をとり、電流は抵抗５０４及びＣＰＥ５０５を流れる枝電流の和を取ってインピーダンスを求めた。本方法の計算結果は図６のプロット（黒丸）で示される。また計算結果の検証のために理論式による予測結果を図６の実線で示す。図５の回路のように２つの回路が直列に並んでいる場合の合成インピーダンスは、２つの回路のインピーダンスの和となる。
【００３３】
図６に示すように理論予測と計算結果はよい一致をしており、本方法を用いて図５に示す回路が適切に計算できていることがわかる。つまり本方法を用いることで、αやＴの値が異なる複数のＣＰＥが存在し、さらに抵抗などＣＰＥ以外の素子が含まれる回路を解くことができることが示された。
【００３４】
［実施例３］
次に実施例３として時間領域で解いていることが分かりやすい例を示す。本実施例では式（６）を用いた。解くべき回路は図７に示すようなＣＰＥ一つからなるものとする。端子７０１には、時刻０［ｓｅｃ］≦ｔ＜３［ｓｅｃ］では電位０［Ｖ］であり、時刻３［ｓｅｃ］≦ｔ≦１０［ｓｅｃ］では電位１００［Ｖ］となる電位をかけるものとする。端子７０２の電位は０［Ｖ］で固定とする。符号７０３はＣＰＥを示す。電流は端子７０１から端子７０２へ流れる向きを正とする。
本実施例の場合には全ての電位は境界条件として与えられているため、ステップＳ１０３では回路の枝電流を求めるだけでよい。
開始時刻はｔ＝０［ｓｅｃ］とし、終了時刻はｔ＝１０［ｓｅｃ］とする。時間刻みはΔｔ＝１．０×１０−２［ｓｅｃ］とする。
【００３５】
ＣＰＥ７０３をＴ＝０．３１８[Ｆ]、α＝１．０とした時の数値計算例を図８Ａに示す。図８Ａの横軸は時刻［ｓｅｃ］であり、縦軸は電流［Ａ］（変位電流と導電電流の和）である。
ＣＰＥ７０３をＴ＝０．３１８[Ｆ／ｓｅｃ０．５]、α＝０．５とした時の数値計算例を図８Ｂに示す。図８Ｂの横軸と縦軸は図８Ａと同様である。
ここで計算結果の検証のために、理論式による予測を行った。ｆ（ｔ）がヘビサイト関
数Ｈ（ｔ）の時、式（１４）と書くことができるのでこれを用いた。但しΓはガンマ関数を示す。
【数１１】
2011145944.tif 000013
本方法の計算結果をプロット（黒丸）で示し、式（１４）を用いた理論式による予測結果を実線で示す。
【００３６】
αが１の時にはＣＰＥは普通のコンデンサと一致する。よって本実施例のような入力電位を与えると電流はデルタ関数となる。数値計算では時間刻みが離散化されているために、厳密にデルタ関数と一致するわけではないが、図８Ａに示す電流はデルタ関数と似た形状をしており、正しく計算されていることがわかる。
それに対してαが０．５の時には分数階微分によって電流はデルタ関数とならず、テールを引くことが予想される。図８Ｂでは予想どおり、電流がテールを引いている計算結果が得られている。また、理論予測と計算結果はよい一致をしており、正しく計算されていることが分かる。
このように本方法を用いることによって、入力波形がヘビサイト的な振る舞いをする場合、即ち周波数領域で解く場合には非常に多くの周波数成分を必要とするような場合にも簡単に解が得られることが示された。
【００３７】
［実施例４］
実施例３では電位を求める必要が無かったが、実施例４では電位を求める必要がある場合を示す。本実施例では式（６）を用いた。解くべき回路は図９に示すようなＣＰＥと抵抗の直列回路とする。端子１００１には、時刻０［ｓｅｃ］≦ｔ＜３［ｓｅｃ］では電位０［Ｖ］であり、時刻３［ｓｅｃ］≦ｔ≦１０［ｓｅｃ］では電位１００［Ｖ］となる電位をかけるものとする。端子１００３の電位は０［Ｖ］で固定とする。１００４は抵抗を示し、抵抗値Ｒ＝１［Ω］とする。１００５はＣＰＥを示す。電流は端子１００１から端子１００３へ流れる向きを正とする。
ステップＳ１０３ではまず節点１００２の電位を変数として求め、その後に回路の枝電流を求める。
開始時刻はｔ＝０［ｓｅｃ］とし、終了時刻はｔ＝１０［ｓｅｃ］とする。時間刻みはΔｔ＝１．０×１０−２［ｓｅｃ］とする。
【００３８】
ＣＰＥ１００５をＴ＝０．３１８[Ｆ]、α＝１．０とした時の数値計算例を図１０Ａに示す。図１０Ａの横軸は時刻［ｓｅｃ］であり、縦軸は電流［Ａ］（変位電流と導電電流の和）である。本方法の計算結果をプロット（黒丸）で示し、理論式による予測結果を実線で示す。αが１の時にはＣＰＥは普通のコンデンサと一致する。単なる抵抗とコンデンサの直列回路であり、時定数τ＝ＲＴの振る舞いを示す。これを用いて理論式による予測を行った。理論予測と計算結果はよい一致をしており、正しく計算されていることが分かる。
【００３９】
ＣＰＥ１００５をＴ＝０．３１８[Ｆ０．５／Ω０．５]、α＝０．５とした時の数値計算例を図１０Ｂに示す。図１０Ｂの横軸と縦軸は図１０Ａと同様である。αが０．５の時には分数階微分によって、電流はα＝１の場合と比べてテールを引くことが予想される。図１０Ｂでは予想どおり、電流がテールを引いている計算結果が得られ、正しく計算されていることが確認できた。
このように本方法を用いることによって、入力波形がヘビサイト的な振る舞いをする場合、即ち周波数領域で解く場合には非常に多くの周波数成分を必要とするような場合にも、簡単に解が得られ、かつ抵抗などＣＰＥ以外の素子が含まれる回路を解くことができる
ことが示された。
【００４０】
［実施例５］
本実施例では式（６）と式（７）の違いについて示す。本実施例では図３に示すような抵抗とＣＰＥが並列に並んだ回路と、図１１に示すような抵抗とコンデンサが並列に並んだ回路を用いる。端子３０１及び端子１３０１には図１２で示されるようにＶ３０１＝１０ｔ［Ｖ］なる電位をかけるものとする。端子３０２及び端子１３０２の電位は０［Ｖ］で固定とする。３０３及び１３０３は抵抗を示し、抵抗値Ｒ＝１［Ω］とする。３０４はＣＰＥを示し、Ｔは両端にかかる電位差の関数となっているものとし、Ｔ＝ｅｘｐ（０．０１ΔＶ３０４）[Ｆ]、α＝１．０とする。但しΔＶ３０４はＣＰＥ３０４にかかる電位差の絶対値である。
１３０４はコンデンサを示し、静電容量Ｃは両端にかかる電位差の関数となっているものとし、Ｃ＝ｅｘｐ（０．０１ΔＶ１３０４）[Ｆ]とする。但しΔＶ１３０４はコンデンサ１３０４にかかる電位差の絶対値である。
電流は端子３０１から端子３０２へ流れる向き及び端子１３０１から端子１３０２へ流れる向きを正とする。
本実施例の場合には全ての電位は境界条件として与えられているため、工程Ｃでは枝電流を求めるだけでよい。
開始時刻はｔ＝０［ｓｅｃ］とし、終了時刻はｔ＝１０［ｓｅｃ］とする。時間刻みはΔｔ＝１．０×１０−２［ｓｅｃ］とする。
【００４１】
式（６）を用いた場合と式（７）を用いた場合のそれぞれについて、上記条件の元で図３に示す回路の数値計算例を図１３に示す。また上記条件の元で図１１に示す回路の数値計算例も同時に図１３に示す。図１３の横軸は時刻［ｓｅｃ］であり、縦軸は電流［Ａ］（変位電流と導電電流の和）である。
既に考察したようにＣＰＥはα＝１の時にはコンデンサと一致する必要がある。図１３を見ると、式（７）を用いた場合にはコンデンサと一致するが、式（６）を用いた場合にはコンデンサと一致しない。従ってＴが両端にかかる電位差の関数となっている場合、即ち時間に依存する関数となっている場合には式（７）を用いた方がより正確に電流を計算することができることがわかる。
【００４２】
［実施例６］
本実施例ではαが１ではない場合についての、式（６）と式（７）の違いについて示す。基本的には実施例５で述べたような回路構成と計算条件を用いるものとするが、ＣＰＥの物性値を変えて計算することとする。
【００４３】
まず、ＣＰＥ３０４の物性値をＴ＝１．０[Ｆ／ｓｅｃ０．２]、α＝０．８とする。式（６）を用いた場合と式（７）を用いた場合のそれぞれについて、上記条件の元で図３に示す回路の数値計算例を図１４Ａに示す。図１４Ａの横軸は時刻［ｓｅｃ］であり、縦軸は電流［Ａ］（変位電流と導電電流の和）である。
図１４Ａに示されるようにＴが時間に依存しない場合には式（６）と式（７）の計算結果は一致し、どちらを用いて計算してもよいことがわかる。
【００４４】
次にＣＰＥ３０４の物性値をＴ＝ｅｘｐ（０．０１ΔＶ３０４）[Ｆ／ｓｅｃ０．２]、α＝０．８とする。式（６）を用いた場合と式（７）を用いた場合のそれぞれについて、上記条件の元で図３に示す回路の数値計算例を図１４Ｂに示す。図１４Ｂの横軸は時刻［ｓｅｃ］であり、縦軸は電流［Ａ］（変位電流と導電電流の和）である。
図１４Ｂに示されるようにＴが時間に依存する場合には式（６）と式（７）の計算結果は異なり、α＝１の時の場合を参考にすれば、式（７）を用いた方がより正確に電流を計算することができると思われる。
【符号の説明】
【００４５】
３０３、５０４、５０６、１００４、１３０３：抵抗
３０４、５０５、５０７、７０３、１００５：ＣＰＥ
１３０４：コンデンサ
【図１】
2011145944.tif 000014
【図２】
2011145944.tif 000015
【図３】
2011145944.tif 000016
【図４】
2011145944.tif 000017
【図５】
2011145944.tif 000018
【図６】
2011145944.tif 000019
【図７】
2011145944.tif 000020
【図８】
2011145944.tif 000021
【図９】
2011145944.tif 000022
【図１０】
2011145944.tif 000023
【図１１】
2011145944.tif 000024
【図１２】
2011145944.tif 000025
【図１３】
2011145944.tif 000026
【図１４】
2011145944.tif 000027