JP2005025651A - 微分方程式の数値解を求める方法およびその方法を用いたプログラム - Google Patents

Abstract

【課題】台形公式などの差分法で微分方程式の数値的な解を求める場合、解が周波数の比較的高い成分と低い成分を含むときは、区間分割を高い周波数成分の短い周期（波長）に比べて十分小さく設定しなくてはならなかった。そのため低い周波数成分が関与する大きい（広い）領域に関して解を求めるためには非常に大きな計算負荷が必要となるという課題が存在した。
【課題を解決するための手段】本発明では上記の従来の方法における課題を解決するための手段として、微分方程式を解く際に、領域を複数の区間に分割し、未知の変数をその区間における（離散）フーリエ級数により近似する。さらに微分方程式を、未知数の（離散）フーリエ級数の係数間の差分方程式へと変換し、この差分方程式を解くことにより、微分方程式を数値的に解く。
【選択図】 図１

Description

【発明の属する技術分野】
本発明は微分方程式を数値的に解く技術の一つである。特に、解となる変数がある周波数近傍の成分を特に多く含むことが判明していて、解のその周波数周辺の成分を効率的に求めつつも、他の周波数成分も考慮に入れたいときに有効である。
【従来の技術】
微分方程式の数値的に解く代表的な方法として、台形公式、前進積分法、バックワードオイラー法などの差分法がある。電気回路が満たすべき微分方程式を、台形公式を用いて数値的に解き、その過渡現象を解析するプログラムとしてＥＭＴＰ（参考文献Ｈ．Ｗ．Ｄｏｍｍｅｌ，”Ｎｏｎｌｉｎｅａｒ ａｎｄ ｔｉｍｅ−ｖａｒｙｉｎｇ ｅｌｅｍｅｎｔｓ ｉｎ ｄｉｇｉｔａｌ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ ｔｒａｎｓｉｅｎｔｓ”，ＩＥＥＥ Ｐｏｗｅｒ Ａｐｐ．Ｓｙｓｔ．，Ｖｏｌ．ＰＡＳ−９０，ｐｐ．２５６１−２５６７，Ｊｕｎｅ １９７１）が有名である。
また、微分方程式を一定の条件の下で数値的に解く別の方法として次のようなものがある。求めようとする変数の解が有る周波数成分（基本周波数）を非常に多く含み、解析の対象が主にその周波数近傍成分に限定されるときは、解がその周波数の成分の波で表されると仮定し，その振幅と位相の変化すると仮定して、解を求める方法が存在する。このような方法の例として、電力系統解析における動的安定度解析をあげることができる。電力系統は電気回路として表現され、電気回路の満たすべき方程式はもともと微分方程式で表されるわけであるから、これらの解法は、もともとの微分方程式を特殊な条件の下に解いていると考えられる。
【発明が解決しようとする課題】
台形公式などの差分法で微分方程式の数値的な解を求める場合，解が周波数の比較的高い成分と、低い成分を含むとき（Ｓｔｉｆｆな系と呼ばれる）は区間分割を、高い周波数成分の短い周期（波長）に比べて十分に小さく設定しなくてはならなかった。そのため低い周波数成分が関与する大きい（広い）領域に関して解を求めるためには非常に大きな計算負荷が必要となるという課題が存在した。
一方、電力系統解析における動的安定度解析では、解がある周波数の波（基本周波数）で近似されることを前提に方程式が構成されている。そのため，方程式の解がその周波数以外の成分を含むような場合は、求まる解の精度が悪化するという問題点があった。また、解析精度を上げることを目的として、区分（時間刻み）を細かくすることも試みられているが、前提として解が基本周波数成分により表されるため、区分を基本周波数の波長（周期）以下に縮めても、その精度が上がるとは限らないという課題が存在した。
【課題を解決するための手段】本発明では上記の従来の方法における課題を解決するための手段として、微分方程式を解く際に、領域を複数の区間に分割し、未知の変数をその区間における、（離散）フーリエ級数により近似する。さらに微分方程式を、未知数の（離散）フーリエ級数の係数間の差分式へと変換し、この差分式を解くことにより、微分方程式を数値的に解く。
【発明実施の形態】
【実施例１】
ｙに関する微分方程式
【数３８】

を初期条件
【数３９】

のもとで領域
【数４０】

について本発明を用いて数値的に解く。
微分方程式を解くためにｘを、区分幅１を用いて区分
【数４１】

を考え、それぞれの区分において、ｙ，ｕを
【数４２】

【数４３】

のようにフーリエ級数で近似する。またこれらの係数を並べてベクトル
【数４４】

【数４５】

を考える。
微分方程式
【数３８】は本発明により次の差分方程式へと置き換えられる。
【数４６】

【数４７】

ここで初期条件は
【数４８】

であり、Ｕ_ｎは
【数４９】

となる。また係数Ａ，Ｋは
【数５０】

【数５１】

により与えられる。
【数４６】は
【数５２】

と変換できるのでこれを
【数４７】と併せて用いると、順次Ｙ_ｎを求めていくことができる。
得られた数値解を
【表１】に示す。
この数値解をグラフで表すと
【図１】の様になる。このなかで、実線が求まった本発明による方法で求まった数値解を
【数４２】を適用してグラフに描いたものである。点線が微分方程式
【数４２】を解析的に解いて求まる解をグラフを描いたものである。グラフから得られた数値解が妥当なものであることが判別できる。
【実施例２】
【図３】のような回路の過渡現象を考える。ここで回路のパラメータとして
【数５３】

【数５４】

【数５５】

とし、初期条件として
【数５６】

とする。
過渡現象を解くためにｘについて区分幅１を用いた区分
【数５７】

を考える。また回路を流れる電流ｉ（ｘ）とノードＱの電圧ｖ_Ｑ（ｘ）は
【数５８】

【数５９】

により近似されるとする。
インダクタンスの電流と電圧の間にはｖ（ｔ）の間には
【数６０】

なる関係が成立する。この関係式は本発明の方法により差分方程式
【数６１】

【数６２】

に変換される。ただし
【数６３】

【数６４】

【数６５】

【数６６】

である。
【数６１】によりインダクタンスは仮想的な区分Ｔ_ｎにおいて
【数６７】、
【数６８】により与えられる仮想的なアドミタンスＧ_Ｌとそれに並列な仮想的電流源Ｈ_ｎにより表される。
【数６７】

【数６８】

一方回路網方程式はノードＱから流れ出る電流の合計が０であることを利用すると
【数６９】

を得る。したがって各区分Ｔ_ｎにおいてＶ_Ｅｎ、Ｈ_ｎが与えられれば
【数７０】

によりＶ_Ｑｎを求めることができる。Ｇ_Ｒは
【数７１】

である。
過渡現象の数値解は次の手続き（ｉｖ）〜（ｖｉｉ）を実施することにより順次求められる。
（ｉｖ）ｎ＝０とし、初期条件よりｈ_ｎ＝０とする。
（ｖ）
【数６８】からＨ_ｎを求める。
（ｖｉ）問題の条件
【数５３】から
【数７２】

であるので、
【数７０】を使用してＶ_Ｑｎを計算する。
（ｖｉｉ）
【数６２】を使用してｈ_ｎ＋１を計算する。ｎを一つ増加させ（ｖ）へ戻る。
実際に数値計算を実施して求まった数値解を
【表２】に示す。
【実施例３】
ｙに関する微分方程式
【数７３】

を初期条件
【数７４】

のもとで領域
【数７５】

について本発明を用いて数値的に解く。解を求める領域を刻み幅Ｔ＝１を用いて
【数７６】

に分割する。請求項３の中で使用されている変数Ｐ，Ｍについてそれぞれ１，６０として各区分の中でｙ（ｘ）を離散フーリエ級数によって近似する。
【数７７】

ここでθ、Ｔ_ｎ，ｍは
【数７８】

【数７９】

とする。
離散フーリエ級数を並べてベクトル
【数８０】

を考えると、本発明により微分方程式
【数７３】は差分方程式
【数８１】

【数８２】

へと変換できる。
ここで初期条件は
【数８３】

であり、Ｕ_ｎは
【数８４】

となる（ｕ（ｘ）を
【数７９】の領域で平均して求まる値に対して、離散フーリエ級数を求める）。また係数Ａ_Ｄ，Ｋ_Ｄは
【数８５】

【数８６】

により与えられる。
【数８１】は
【数８７】

と変換できるので、これを
【数８２】と併せて用いると、順次Ｙ_Ｄｎを求めていくことができる。数値計算を実施して求まった数値解を
【表３】に示す。
【実施例４】
【図３】のような回路の過渡現象を考える。ここで回路のパラメータとして
【数５３】（再掲）

【数５４】（再掲）

【数５５】（再掲）

とし、初期条件として
【数５６】（再掲）

とする。
過渡現象を解くためにｘについて区分幅１を用いた区分
【数５７】（再掲）

を考える。さらに本発明で使用するパラメータＭ，Ｐをここでは６０，１と選んで
【数７９】（再掲）

とする。
回路を流れる電流ｉ（ｘ）とノードＱの電圧ｖ_Ｑ（ｘ）は
【数８８】

【数８９】

により近似する。
インダクタンスの電流と電圧の間にはｖ（ｔ）の間には
【数６０】（再掲）

なる関係が成立する。この関係式は本発明の方法により差分方程式
【数９０】

【数９１】

に変換される。ただし
【数９２】

【数９３】

【数９４】

【数９５】

である。
【数９０】によりインダクタンスは仮想的な区分Ｔ_ｎにおいて
【数９６】、
【数９７】により与えられる仮想的なアドミタンスＧ_Ｌとそれに並列な仮想的電流源Ｈ_ｎにより表される。
【数９６】

【数９７】

一方回路網方程式はノードＱから流れ出る電流の合計が０であることを利用すると
【数９８】

を得る。したがって各区分Ｔ_ｎにおいてＶ_ＥＤｎ、Ｈ_Ｄｎが与えられれば
【数９９】

によりＶ_Ｑｎを求めることができる。Ｇ_ＲＤは
【数１００】

である。
過渡現象の数値解は次の手続き（ｉｖ）〜（ｖｉｉ）を実施することにより順次求められる。
（ｉｖ）ｎ＝０とし、初期条件よりｈ_Ｄｎ＝０とする。
（ｖ）
【数９７】からＨ_Ｄｎを求める。
（ｖｉ）問題の条件
【数５３】から
【数１０１】

であるので、
【数９９】を使用してＶ_ＱＤｎを計算する。
（ｖｉｉ）
【数９１】を使用してｈ_Ｄｎ＋１を計算する。ｎを一つ増加させ（ｖ）へ戻る。
実際に数値計算を実施して求まった数値解（Ｉ_Ｄｎ）を
【表４】に示す。Ｉ_{Ｄｎ，−１}はＩ_Ｄｎ，１に複素共役なので表からは省略してある。
【実施例５】
ｙに関する微分方程式
【数１０２】

を初期条件
【数１０３】

のもとで領域
【数１０４】

について本発明を用いて数値的に解く。解を求める領域を刻み幅Ｔ＝１０を用いて
【数１０５】

に分割する。それぞれの区分において、ｙ，ｕを
【数１０６】

【数１０７】

のようにフーリエ級数で近似する。またこれらの係数を並べてベクトル
【数１０８】

【数１０９】

を考える。
微分方程式
【数１０２】は本発明により次の差分方程式へと置き換えられる。
【数１１０】

【数１１１】

ここで初期条件は
【数１１２】

であり、Ｕ_ｎは
【数１１３】

となる。また係数Ａ，Ｋは
【数１１４】

【数１１５】

により与えられる。
【数１１０】は
【数１１６】

と変換できるのでこれを
【数１１１】と併せて用いると、順次Ｙ_ｎを求めていくことができる。
得られた数値解を
【表５】に示す。
【実施例６】
【図３】のような回路の過渡現象を考える。ここで回路のパラメータとして
【数１１７】

【数１１８】

【数１１９】

とし、初期条件として
【数１２０】

とする。
過渡現象を解くためにｘについて区分幅１０を用いた区分
【数１２１】

を考える。また回路を流れる電流ｉ（ｘ）とノードＱの電圧ｖ_Ｑ（ｘ）は
【数１２２】

【数１２３】

により近似されるとする。
インダクタンスの電流と電圧の間にはｖ（ｔ）の間には
【数１２４】

なる関係が成立する。この関係式は本発明の方法により差分方程式
【数１２５】

【数１２６】

に変換される。ただし
【数１２７】

【数１２８】

【数１２９】

【数１３０】

である。
【数１２５】によりインダクタンスは仮想的な区分Ｔ_ｎにおいて
【数１３１】、
【数１３２】により与えられる仮想的なアドミタンスＧ_Ｌとそれに並列な仮想的電流源Ｈ_ｎにより表される。
【数１３１】

【数１３２】

一方回路網方程式はノードＱから流れ出る電流の合計が０であることを利用すると
【数１３３】

を得る。したがって各区分Ｔ_ｎにおいてＶ_Ｅｎ、Ｈ_ｎが与えられれば
【数１３４】

によりＶ_Ｑｎを求めることができる。Ｇ_Ｒは
【数１３５】

である。
過渡現象の数値解は次の手続き（ｉｖ）〜（ｖｉｉ）を実施することにより順次求められる。
（ｉｖ）ｎ＝０とし、初期条件よりｈ_ｎ＝０とする。
（ｖ）
【数１３２】からＨ_ｎを求める。
（ｖｉ）問題の条件
【数１１７】から
【数１３６】

であるので、
【数１３４】を使用してＶ_Ｑｎを計算する。
（ｖｉｉ）
【数１２６】を使用してｈ_ｎ＋１を計算する。ｎを一つ増加させ（ｖ）へ戻る。
実際に数値計算を実施して求まった数値解を
【表６】に示す。
【実施例７】
請求項１の方法において
【数１３７】

【数１３８】

なる変数変換を施すことを考える。ここでＷは例として

により与える。すると、
【数６】、
【数７】はそれぞれ
【数１３９】

【数１４０】

と変換される。ただし、
【数１４１】

【数１４２】

【数１４３】

である。得られたこれらの差分式を用いて微分方程式を解く方法は実施例１と同様である。
【実施例８】実施例１の方法で
【数１４４】

【数１４５】

とする。
【数１４６】

【数１４７】

【数１４８】

【数１４９】

とすれば
【数１５０】

【数１５１】

【数１５２】

となり差分方程式
【数６】、
【数７】はそれぞれ
【数１５３】

【数１５４】

と変換される。ただし、
【数１５５】

【数１５６】

【数１５７】

である。得られたこれらの差分式を用いて微分方程式を解く方法は実施例１と同様である。
【発明の効果】微分方程式の解となる変数がある周波数近傍の成分を特に多く含むことが判明していている場合に、本発明により、解のその周波数周辺の成分を効率的に求めつつも、他の周波数成分も考慮に入れて効率的に数値的に解くことができる。
【図面の簡単な説明】
【図１】実施例１の結果得られた数値解をグラフ化したもの
【図２】実施例１の結果得られた数値解
【図３】過渡現象を数値解析する回路例
【図４】実施例２の結果得られた数値解
【図５】実施例３の結果得られた数値解
【図６】実施例４の結果得られた数値解
【図７】実施例５の結果得られた数値解
【図８】実施例６の結果得られた数値解
【符号の説明】
１−１は本発明による数値解をグラフ化したもの
１−２は解析解
２−１は実施例１の解析結果
３−１は電圧源
３−２は抵抗
３−３はインダクタンス
４−１は実施例２の解析結果
５−１は実施例３の解析結果
６−１は実施例４の解析結果
７−１は実施例５の解析結果
８−１は実施例６の解析結果

Claims (8)

1. 微分方程式を数値的に解く際に次の処理（ｉ）〜（ｉｖ）を持つことを特徴とする方法、その方法を用いた微分方程式の数値解をもとめるプログラム。
   （ｉ）解を求めたい領域を刻み幅Ｔを用いて複数の区間
   【数１】に分割する。
   

   

   （ｉｉ）それぞれの区間で、解析の対象となる変数を、フーリエ級数によって近
   似する（
   【数２】、
   【数３】）。
   

   

   

   ただし、ここでＮは１以上の整数、ωは
   

   

   で与えられる。
   （ｉｉｉ）二つの変数に関係
   

   

   が存在するとき、ステップ（ｉｉ）で表したフーリエ級数の間に
   

   

   

   なる関係式が存在するとして微分方程式を差分化する。ただし、ここで
   

   

   

   

   

   である。ｈ_ｎは前の区分までの解に由来する履歴項である。また右肩の^＊は転置をあらわす。
   （ｉｖ）（ｉｉｉ）で求めた関係式を順次解くことにより微分方程式の数値解をもとめる。
2. 電気回路の過渡現象を数値的に解く際に次の処理（ｉ）〜（ｖｉｉ）を持つことを特徴とする方法、およびその方法を用いた電気回路の過渡現象を数値的に解くプログラム。
   （ｉ）解を求めたい時間を時間刻みＴにより区間
   【数１】に分割する（ｘは時間を表す）。
   （ｉｉ）それぞれの区間において、回路の電気的な諸量（電圧、電流など）をその区間におけるフーリエ級数の係数により近似する。（
   【数２】、
   【数３】）（ｉｉｉ）個々の電気回路の素子の２つの諸量の間に関係
   【数５】が存在するとき、ステップ（ｉｉ）で表したフーリエ級数の間に
   【数６】、
   【数７】で表される関係式が存在するとして微分方程式を差分化する。ただし、
   【数６】、
   【数７】で使用されるＦ_ｎ、Ｇ_ｎ、Ａ、Ｋは
   【数８】、
   【数９】、
   【数１０】、
   【数１１】により与えられる。ｈ_ｎは前の区分までの解に由来する履歴項である。
   （ｉｖ）適切な初期化方法により開始時間刻みにおける個々の電気素子に対する履歴項ｈ_ｎをもとめる。
   （ｖ）（ｉｉｉ）により求められた差分式と履歴項ｈ_ｎを用いて、個々の電気素子の両端間の電圧が
   

   

   素子を流れる電流が
   

   

   により与えられるとき、その関係を、仮想的なアドミタンスＧとそれに並列な電流源
   

   

   を用いて
   

   

   による等価回路により表現する。ここで
   

   

   

   

   である。
   （ｖｉ） 解析しようとする電気回路中のすべての素子を（ｖ）の方法で等価回路として表現した後、その時間刻みについて回路の各ノードの電圧を求める。ここでノードの電圧は
   

   

   として与えられる。
   （ｖｉｉ）（ｖｉ）にて求まった回路のノード電圧を使用し、各電気素子の履歴項ｈ_ｎの次の時間区分における値を求め、時間刻みを一つ進める。その後（ｖ）へ戻る。
3. 微分方程式を数値的に解く際に次の処理（ｉ）〜（ｉｖ）を持つことを特徴とする方法、その方法を用いた微分方程式の数値解をもとめるプログラム。
   （ｉ）解を求めたい領域を刻み幅Ｔを用いて複数の区間
   【数１】に分割する。
   【数１】（再掲）
   

   

   

   象となる変数を、離散フーリエ級数によって近似する（
   【数２０】、
   【数２１】）。
   

   

   

   

   Ｔ_ｎ，ｍは
   

   

   で与えられる。
   （ｉｉｉ）二つの変数に関係
   【数５】（再掲）
   

   

   が存在するとき、ステップ（ｉｉ）で表したフーリエ級数の間に
   

   

   

   なる関係式が存在するとして微分方程式を差分化する。ただし、ここで
   

   

   

   

   

   である。ｈ_Ｄｎは前の区分までの解に由来する履歴項である。
   （ｉｖ）（ｉｉｉ）で求めた関係式を順次解くことにより微分方程式の数値解をもとめる。
4. 電気回路の過渡現象を数値的に解く際に次の処理（ｉ）〜（ｖｉｉ）を持つことを特徴とする方法、およびその方法を用いた電気回路の過渡現象を数値的に解くプログラム。
   （ｉ）解を求めたい時間を時間刻みＴにより区間
   【数１】に分割する（ｘは時間を表す）。
   

   

   路の電気的な諸量（電圧、電流など）をその区間における離散フーリエ級数の係数により近似する。（
   【数２０】、
   【数２１】、
   【数２２】、
   【数２３】）（ｉｉｉ）個々の電気回路の素子の２つの諸量の間に関係
   【数５】が存在するとき、ステップ（ｉｉ）で表した離散フーリエ級数の間に
   【数２４】、
   【数２５】で表される関係式が存在するとして微分方程式を差分化する。ただし、
   【数２４】、
   【数２５】で使用されるＦ_Ｄｎ、Ｇ_Ｄｎ、Ａ_Ｄ、Ｋ_Ｄは
   【数２６】、
   【数２７】、
   【数２８】、
   【数２９】により与えられる。ｈ_Ｄｎは前の区分までの解に由来する履歴項である。
   （ｉｖ）適切な初期化方法により開始時間刻みにおける個々の電気素子に対する履歴項ｈ_Ｄｎをもとめる。
   （ｖ）（ｉｉｉ）により求められた差分式と履歴項ｈ_Ｄｎを用いて、個々の電気素子の両端間の電圧が
   

   

   素子を流れる電流が
   

   

   により与えられるとき、その関係を、仮想的なアドミタンスＧ_Ｄとそれに並列な電流源
   

   

   を用いて
   

   

   による等価回路により表現する。ここで
   

   

   

   

   である。
   （ｖｉ） 解析しようとする電気回路中のすべての素子を（ｖ）の方法で等価回路として表現した後、その時間刻みについて回路の各ノードの電圧を求める。ここでノードの電圧は
   

   

   として与えられる。
   （ｖｉｉ）（ｖｉ）にて求まった回路のノード電圧を使用し、各電気素子の履歴項ｈ_ｎの次の時間区分における値を求め、時間刻みを一つ進める。その後（ｖ）へ戻る。
5. 請求項１、請求項３の方法において最高次数の高調波以外の成分のうち一つもしくはいくつかを０と仮定し、縮小した行列、ベクトルを用いて、微分方程式を解く方法、及びその方法を用いた微分方程式の数値解をもとめるプログラム。
6. 請求項２、請求項４、の方法において最高次数の高調波以外の成分のうち一つもしくはいくつかを０と仮定し、縮小した行列、ベクトルを用いて、電気回路の過渡現象を数値的に解く方法、およびその方法を用いた電気回路の過渡現象を数値的に解くプログラム。
7. 請求項１、請求項２、請求項３、請求項４、請求項５、請求項６の方法、プログラムについて（離散）フーリエ級数を並べたベクトルに対して、逆行列をもつ行列を乗じて、変換を施したもの。
8. 請求項１、請求項２、請求項３、請求項４、請求項５、請求項６、請求項７の方
   法、プログラムについて
   【数１】で表される区分をシフトさせて適用したもの。

Images