JPH09160903A

JPH09160903A - 物理量分布推定方法

Info

Publication number: JPH09160903A
Application number: JP31682195A
Authority: JP
Inventors: Kiyoshi Yoda; 潔依田
Original assignee: Mitsubishi Electric Corp
Current assignee: Mitsubishi Electric Corp
Priority date: 1995-12-05
Filing date: 1995-12-05
Publication date: 1997-06-20

Abstract

(57)【要約】【課題】密行列を解く必要があるため、行列方程式の
計算に長い時間を必要とする。【解決手段】与えられた支配方程式を境界積分方程式
に変換し（ＳＴ２）、その境界積分方程式を離散化方程
式に変換し（ＳＴ３）、その離散化方程式を行列方程式
Ａ・ｘ＝ｂに変換し（ＳＴ４）、得られた行列方程式Ａ
・ｘ＝ｂをウェーブレット変換し（ＳＴ４）、変換され
た行列方程式Ａ_t ・ｘ_t ＝ｂ_t を、その行列Ａ_t 中の絶
対値が閾値以下の行列要素をゼロに置き換えて近似的に
解き（ＳＴ６）、その近似解ｘ_t をウェーブレット逆変
換して近似的な物理量分布ｘ^* を求め（ＳＴ７）、繰り
返し計算によってその近似的な物理量分布ｘ^* の補正を
行う（ＳＴ８）。

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】この発明は、電位分布、電磁
界分布、応力分布、変位分布、温度分布などをはじめと
する各種物理量分布を高速に推定する方法に関するもの
であり、さらに詳しくは、境界要素法と呼ばれる物理量
分布の推定方法を著しく高速化した物理量分布推定方法
に関するものである。

【０００２】

【従来の技術】図９は、例えば「境界要素解析」（榎薗
正人著、昭和６１年培風館発行）の１９頁に示され
た、従来の境界要素法による物理量分布推定方法をわか
りやすく改変して示す説明図である。図において、１は
対象システムの物理現象を表す支配方程式を与えるため
の支配方程式設定部、２は与えられた支配方程式を境界
積分方程式に変換する境界積分方程式変換部であり、３
はこの境界積分方程式への変換に用いられる積分定理で
ある。４はその境界積分方程式を離散化方程式に変換す
る離散化方程式変換部であり、５はこの離散化方程式へ
の変換の際の、境界を複数の要素に分割した境界要素で
ある。６はその離散化方程式より連立一次方程式（実際
には行列方程式Ａ・ｘ＝ｂ）を導出する連立一次方程式
変換部であり、７はこの行列方程式Ａ・ｘ＝ｂをガウス
の消去法で解いて未知ベクトルｘを得る行列方程式計算
部である。

【０００３】以下、この図９に従って従来の物理量分布
推定方法について説明する。まず最初に、支配方程式設
定部１において、対象システムの物理現象を表す支配方
程式を微分方程式および境界条件式によって与える。次
に境界積分方程式変換部２において、この支配方程式設
定部１にて与えられた支配方程式を、積分定理３を用い
て対象システムの境界領域における境界積分方程式に変
換する。次に境界をたくさんの境界要素５に分割し、境
界積分方程式変換部２で得られた境界積分方程式より離
散化方程式変換部４にて離散化方程式を作る。次に連立
一次方程式変換部６において、その離散化方程式変換部
４で得られた離散化方程式より連立一次方程式（実際に
は行列方程式Ａ・ｘ＝ｂ）を導出する。そして行列方程
式計算部７において、連立一次方程式変換部６で得られ
た行列方程式Ａ・ｘ＝ｂより、ガウスの消去法などを用
いて未知ベクトルｘを解く。ここで、上記離散化方程式
変換部４、連立一次方程式変換部６、および行列方程式
計算部７は、コンピュータ・ブログラミングによって実
現されている。

【０００４】なお、上記行列方程式Ａ・ｘ＝ｂ中の、Ａ
は対象システムに依存するシステム行列、ｘは前記境界
上の未知の物理量分布を表す未知ベクトル、ｂは定数ベ
クトルをそれぞれ示している。

【０００５】ここで、実際の計算においては、対象シス
テムに依存するシステム行列Ａを導出するための時間に
比べて、行列方程式Ａ・ｘ＝ｂを解くための時間が大変
長くなる。これは、ゼロ要素がほとんどない行列（通常
密行列と呼ばれている）に関する行列方程式を解く必要
があり、さらにこの行列は非対称行列であるため、計算
効率の悪いガウスの消去法、またはＬＵ分解という方法
を用いる必要があるからである。

【０００６】

【発明が解決しようとする課題】従来の物理量分布推定
方法は以上のように構成されているので、ゼロ要素がほ
とんどない密行列と呼ばれる行列に関する行列方程式を
解く必要があり、さらにこの行列は非対称行列であるた
め、計算効率の悪いガウスの消去法やＬＵ分解などの方
法により解を求めざるを得ず、その計算時間が非常に長
くなるという課題があった。

【０００７】この発明は上記のような課題を解決するた
めになされたもので、ウェーブレット変換処理により、
密行列をゼロ要素が多く含まれる疎行列に変換して行列
方程式を解くことによって、その計算時間を大幅に短縮
し、高速に物理量分布を推定することができる物理量分
布推定方法を得ることを目的とする。

【０００８】請求項１記載の発明に係る物理量分布推定
方法は、与えられた支配方程式を境界積分方程式に、そ
の境界積分方程式を離散化方程式に、その離散化方程式
を行列方程式Ａ・ｘ＝ｂに順次変換し、得られた行列方
程式Ａ・ｘ＝ｂをウェーブレット変換したＡ_t ・ｘ_t ＝
ｂ_t を、その行列Ａ_t 中の絶対値が閾値以下の行列要素
をゼロに置き換えることによって近似的に解き、得られ
た近似解ｘ_t をウェーブレット逆変換することで近似的
な物理量分布ｘ^* を求め、その近似的な物理量分布ｘ^*
を繰り返し計算によつて補正することによって解の精度
を向上させたものである。

【０００９】請求項２記載の発明に係る物理量分布推定
方法は、ドビッシーの４係数型ウェーブレットを用いて
ウェーブレット変換を行うようにしたものである。

【００１０】請求項３記載の発明に係る物理量分布推定
方法は、ウェーブレット変換された行列方程式Ａ_t ・ｘ
_t ＝ｂ_t を近似的に解く際の閾値を、行列Ａ_t の各行列
要素の絶対値の最大値の０．０１％ないし０．００００
１％の範囲に選定したものである。

【００１１】請求項４記載の発明に係る物理量分布推定
方法は、前処理付き双共役勾配法を用いて行列方程式Ａ
_t ・ｘ_t ＝ｂ_t を近似的に解くようにしたものである。

【００１２】請求項５記載の発明に係る物理量分布推定
方法は、境界積分方程式の離散化方程式への変換に際し
て、離散化したときの境界要素数が２のべき乗となるよ
うにしたものである。

【００１３】請求項６記載の発明に係る物理量分布推定
方法は、繰り返し計算による近似的な物理量分布ｘ^* の
補正を、近似的な物理量分布ｘ^* とシステム行列Ａおよ
び定数ベクトルｂから算出された誤差ｄ_b より、誤差行
列方程式Ａ・ｄ_x ＝ｄ_b に基づいて誤差ｄ_x を求め、近
似的な物理量分布ｘ^* からその誤差ｄ_x を差し引いて得
られた真の解ｘ₁ を近似的な物理量分布ｘ^* に代入し、
上記処理を再度実行することによって解ｘ₂ を求め、こ
の解ｘ₂ と真の解ｘ₁ との相対誤差が所定の設定値より
も小さくなるまで、その解ｘ₂ を真の解ｘ₁ に代入して
解ｘ₂ の計算を繰り返すことによって行うようにしたも
のである。

【００１４】請求項７記載の発明に係る物理量分布推定
方法は、行列方程式をウェーブレット変換して近似的に
解き、得られた近似解をさらにウェーブレット逆変換す
るという手順に従って、誤差行列方程式Ａ・ｄ_x ＝ｄ_b
から誤差ｄ_x を求めるようにしたものである。

【００１５】請求項８記載の発明に係る物理量分布推定
方法は、真の解ｘ₁ と解ｘ₂ との相対誤差の検定のため
の所定の設定値を、０．０１％ないし０．０００１％の
範囲で選定するようにしたものである。

【００１６】

【発明の実施の形態】以下、この発明の実施の一形態を
説明する。実施の形態１．図１はこの発明の実施の形態１による物
理量分布推定方法を示すフローチャートである。図にお
いて、ＳＴ１は対象システムの支配方程式を与える第１
のステップ、ＳＴ２はその支配方程式を境界積分方程式
に変換する第２のステップ、ＳＴ３はその境界積分方程
式を離散化方程式に変換する第３のステップ、ＳＴ４は
その離散化方程式を行列方程式Ａ・ｘ＝ｂに変換する第
４のステップである。なお、これら第１のステップＳＴ
１ないし第４のステップＳＴ４による、支配方程式から
行列方程式Ａ・ｘ＝ｂを得るまでの処理は、図９に示し
た従来の物理量分布推定方法における、支配方程式設定
部１、境界積分方程式変換部２、離散化方程式変換部
４、および連立一次方程式変換部６による処理と同等で
ある。

【００１７】また、ＳＴ５は行列方程式Ａ・ｘ＝ｂの両
辺をウェーブレット変換して、新たな行列方程式Ａ_t ・
ｘ_t ＝ｂ_t を得る第５のステップであり、ＳＴ６はその
行列Ａ_t の行列要素中の、絶対値が所定の閾値以下の行
列要素をゼロに置き換え、近似的にこの行列方程式Ａ_t
・ｘ_t ＝ｂ_t を解いて、近似解ｘ_t を求める第６のステ
ップである。ＳＴ７は得られた近似解ｘ_t に対してウェ
ーブレット逆変換を行い、近似的な物理量分布ｘ^* を求
める第７のステップであり、ＳＴ８はその近似的に求め
た物理量分布ｘ^* を繰り返し計算によつて補正する第８
のステップである。

【００１８】次に動作について説明する。ここで、図２
は今回考える対象システムの一例を示す説明図であり、
４つの辺１０、１１、１２、１３で囲まれた解析領域は
正方形の内部である。また、図２に示した対象システム
では、境界領域は当該正方形の４つの辺１０、１１、１
２、１３であり、境界領域（この場合は境界線）上に要
素番号１から３２が付与された境界要素がある。また、
辺１０は電位ｕ＝０Ｖ、辺１２は電位ｕ＝３００Ｖに固
定されているものとする。なお、他の２つの辺１１、１
３は、電界ベクトルが各辺に平行に向き、電界の辺に垂
直な成分ｄｕ／ｄｎはゼロと仮定する。

【００１９】まず、第１のステップＳＴ１において、対
象システムの物理現象を支配する微分方程式と境界条件
式を支配方程式として与える。この支配方程式として
は、ラプラス方程式、ポアソン方程式、ヘルムホルツ方
程式、あるいは波動方程式などが考えられる。例えば、
図２に示した対象システムでは、支配方程式として次式
で示すラプラス方程式が選ばれる。

【００２０】

【数１】

【００２１】ここで、境界条件は、辺１０は電位ｕ＝０
Ｖ、辺１２は電位ｕ＝３００Ｖであり、他の２辺につい
ては電界の各辺に垂直な成分ｄｕ／ｄｎがゼロというも
のである。なお、この明細書においては、ｊをｋで偏微
分することをｄｊ／ｄｋと表記する。

【００２２】次に第２のステップＳＴ２において、上記
第１のステップＳＴ１で与えられた支配方程式を積分定
理（グリーンの定理）を用いて、境界領域に対する次式
の境界積分方程式に変換する。

【００２３】

【数２】

【００２４】ここで、ｕ_i は境界領域上の点ｉにおける
ポテンシャルｕを表わしている。図２に示した対象シス
テムでは、境界領域は正方形の４つの辺１０、１１、１
２、１３である。また、ｑはポテンシャルｕの境界に垂
直な方向の傾きｄｕ／ｄｎを表わしている。なお、積分
は境界領域上の積分であり、図２では正方形の４つの辺
１０、１１、１２、１３に沿った線積分となる。さら
に、ｕ^* は領域内部に点エネルギー源があったときの解
（基本解）を表わし、ｑ^* はｄｕ^* ／ｄｎで定義され
る。

【００２５】次に第３のステップＳＴ３において、上記
第２のステップＳＴ２で得られた境界積分方程式を離散
化方程式に変換する。そのために、境界領域を多くの境
界要素に分割（図２に示した例では、境界領域である４
つの辺１０ないし辺１３は要素番号１ないし３２が付さ
れた３２の境界要素に分割）して、積分演算を加算演算
で近似する。この結果は次式の行列表現で記述できる。Ｈ・ｕ＝Ｇ・ｑここに、ＨとＧは対象システムを記述する行列である。
なお、ウェーブレット変換するために、境界要素の数は
２のべき乗個にすることが好ましい。

【００２６】次に第４のステップＳＴ４において、上記
第３のステップＳＴ３で得られた離散化方程式を行列方
程式に変換する。ここで、上記Ｈ・ｕ＝Ｇ・ｑの形は解
くのに不便である。なぜならば、変数ｕとｑの一部は境
界条件として与えられており、それら以外が未知変数だ
からである。そこで、未知変数の部分だけをまとめ直し
て、未知の物理量分布を表す未知ベクトルｘを考える。
その結果、次式の行列方程式に変換することができる。Ａ・ｘ＝ｂここで、Ａは対象システムに依存するシステム行列であ
り、ｂは定数ベクトルである。なお、この第１のステッ
プＳＴ１から第４のステップＳＴ４までの処理は「境界
要素解析」（榎薗正人著、昭和６１年培風館発行）な
どに示された、従来の物理量分布推定方法における処理
と同様である。

【００２７】次に第５のステップＳＴ５において、上記
第４のステップＳＴ４で得られた行列方程式Ａ・ｘ＝ｂ
において、行列Ａおよびベクトルｘ、ｂに対してそれぞ
れウェーブレット変換を行う。ここで、このウェーブレ
ット変換とは、データの平均値の情報と微分係数の情報
を分離して効率よく保存する変換で、画像の圧縮などに
よく利用されているものであり、ケンブリッジ大学出版
から１９９２年に発刊された、ウイリアム・プレス（Wi
lliam H.Press ）等による「ニューメリカル・レシピー
ズ (Numerical Recipes)」第２版の５８４頁〜５９９頁
などに記載されている公知のものであるため、ここでは
その詳細な説明は割愛する。なお、具体的には、ウェー
ブレット変換行列Ｗおよびその転置行列Ｗ^T に対して、
次の変換式で実行する。Ａ_t ＝Ｗ・Ａ・Ｗ^T ｘ_t ＝Ｗ・ｘｂ_t ＝Ｗ・ｂ

【００２８】なお、このウェーブレット変換をするため
には、行列およびベクトルのサイズは２のべき乗が計算
効率上好ましい。すなわち、３２、６４、１２８、２５
６、５１２、１０２４、２０４８、・・・などが好まし
いサイズである。また、この実施の形態１ではウェーブ
レット変換として、公知のドビッシーの４係数型ウェー
ブレットを用いている。このドビッシーの４係数型ウェ
ーブレットを用いた場合と、より複雑なドビッシーの１
２係数型ウェーブレットを用いた場合について比較して
みたが、その結果は同等であった。他にも種々のウェー
ブレットがあるが、どれを用いても、ここで記載したも
のと同様な結果が得られる。従って、このようにドビッ
シーの４係数型ウェーブレットを用いれば、より簡単な
計算によってウェーブレット変換を行うことができる。

【００２９】上記ウェーブレット変換の結果、上記行列
方程式Ａ・ｘ＝ｂは次式で示す行列方程式に変換され
る。Ａ_t ・ｘ_t ＝ｂ_t

【００３０】次に第６のステップＳＴ６において、上記
第５のステップＳＴ５で得られた行列方程式Ａ_t ・ｘ_t
＝ｂ_t について、その行列Ａ_t の行列要素中の絶対値が
所定の閾値以下の行列要素をゼロに置き換えることによ
り、近似的にそれを解いて近似解ｘ_t を求める。ここ
で、図３および図４はそれぞれ、システム行列Ａおよび
そのウェーブレット変換後の行列Ａ_t を疑似３次元表示
した説明図である。図３からシステム行列Ａは複雑な構
造をしている様子がわかる。また、ゼロ要素はほとんど
なく、非対称行列であることもわかる。一方、図４のウ
ェーブレット変換後の行列Ａ_t は、対角成分が主要にな
り、全体にゼロに近い要素が非常に多くなっていること
がわかる。例えば、図４では絶対値が０．５以上の行列
要素は１２０個（約１０％のデータに相当）しかないの
で、絶対値が０．５以下を０に置き換えると、ゼロ要素
が多い疎行列に近似できる。ここで、このゼロに置き換
える際の閾値を大きくしすぎると、行列方程式Ａ_t ・ｘ
_t ＝ｂ_t が解けなくなることがわかった。経験的には、
行列Ａ_t の各行列要素の絶対値の最大値の０．０１％な
いし０．００００１％程度を閾値とすると、最も高速に
解を求められることがわかった。

【００３１】このような近似を行って、第５のステップ
ＳＴ５でウェーブレット変換された行列方程式Ａ_t ・ｘ
_t ＝ｂ_t を解く方法としては、公知の疎行列用の高速計
算法を用いることができる。この実施の形態１ではその
一つである、前処理付き双共役勾配法（preconditioned
bi-conjugate gradient method ）を用いてこの行列方
程式Ａ_t ・ｘ_t ＝ｂ_t を解いている。ここで、この前処
理付き双共役勾配法とは、行列方程式Ａ・ｘ＝ｂを解く
際に、Ａ・ｘ−ｂがなるべく小さくなる様に、くり返し
反復計算して未知ベクトルｘの値を求めることにより、
この行列方程式Ａ・ｘ＝ｂを簡略的に解くものであり、
例えば、前述した「ニューメリカル・レシピーズ」第２
版（ウイリアム・プレス外ケンブリッジ大学出版１
９９２年）などにも記載されているポピュラーなもので
あるため、ここではその詳細な説明を割愛する。

【００３２】なお、この疎行列用の高速計算法として
は、上記前処理付き双共役勾配法以外にも数多くの手法
が存在しており、いずれの手法を用てもよいが、この前
処理付き双共役勾配法を用いることによって、行列方程
式Ａ_t ・ｘ_t ＝ｂ_t をより高速に解くことができる。

【００３３】次に第７のステップＳＴ７において、上記
第６のステップＳＴ６で行列方程式Ａ_t ・ｘ_t ＝ｂ_t を
解くことによって得られた近似解ｘ_t に対して、次式に
よるウェーブレット逆変換を行い、近似的な物理量分布
ｘ^* を求める。ｘ^* ＝Ｗ^T ・ｘ_t

【００３４】ここで、図５は上述のようにして行列方程
式Ａ_t ・ｘ_t ＝ｂ_t を解き、得られた近似解ｘ_t をウェ
ーブレット逆変換し、近似的な物理量分布（この例では
電位と電界成分の分布）ｘ^* を求めた結果を示す説明図
である。また、図６はそれと対比するために従来の方法
で計算した結果を示す説明図である。なお、図５および
図６において、横軸の座標１ないし３２は図２の要素番
号１ないし３２に対応しており、横軸の座標１ないし
８、および１７ないし２４は各要素中心における電位を
表している。一方、横軸座標９ないし１６、および２５
ないし３２は各要素中心における辺に垂直な電界成分を
表している。この図５と図６とを比較すると、近似的に
は解が求められたが、その精度が悪いことがわかる。こ
れは、ウェーブレット変換後、絶対値が０．５以下の行
列要素を０に置き換えたことによるものである。

【００３５】次に第８のステップＳＴ８において、上記
第７のステップＳＴ７で近似的に求めた物理量分布ｘ^*
の補正を行ってその精度を向上させる。すなわち、近似
的に求めた物理量分布ｘ^* を繰り返し計算によって補正
することにより、その近似的な物理量分布ｘ^* の精度を
向上させるわけである。

【００３６】実施の形態２．次に、上記第８のステップ
ＳＴ８による近似的に求めた物理量分布ｘ^* の補正につ
いての具体的な手法を実施の形態２として説明する。図
７は近似的な物理量分布ｘ^* の繰り返し計算による補正
の手順を示すフローチャートである。図において、ＳＴ
９は近似的な物理量分布ｘ^* とシステム行列Ａおよび定
数ベクトルｂから誤差ｄ_b を計算する第９のステップ、
ＳＴ１０はこの第９のステップＳＴ９で算出された誤差
ｄ_b より誤差行列方程式に基づいて誤差ｄ_x を求める第
１０のステップであり、ＳＴ１１はこの第１０のステッ
プＳＴ１０で得られた誤差ｄ_xを近似的な物理量分布ｘ^*
から差し引いて、真の解ｘ₁ を求める第１１のステッ
プである。ＳＴ１２はこの第１１のステップＳＴ１１で
得られた真の解ｘ₁ を近似的な物理量分布ｘ^* に代入し
て再度第９のステップないし第１１のステップを実行
し、より精度の高い解ｘ₂ を求める第１２のステップで
あり、ＳＴ１３はこの第１２のステップＳＴ１２にて得
られた解ｘ₂ と第１１のステップＳＴ１１で得られた真
の解ｘ₁ とを比較し、それらの相対誤差が所定の設定値
よりも小さければ（例えば、設定値未満であれば）一連
の計算処理を終了し、設定値よりも大きければ（例え
ば、設定値以上であれば）解ｘ₂ を真の解ｘ₁ に代入し
た後、処理を第１２のステップ１２に戻す第１３のステ
ップである。

【００３７】次に動作について説明する。まず第９のス
テップＳＴ９において、図１に示す第７のステップＳＴ
７で近似的に求められた物理量分布ｘ^* と、第４のステ
ップＳＴ４で得られた行列方程式Ａ・ｘ＝ｂにおけるシ
ステム行列Ａおよび定数ベクトルｂから、次式を用いて
誤差ｄ_b の計算を行う。ｄ_b ＝Ａ・ｘ^* −ｂ

【００３８】次に第１０のステップＳＴ１０において、
上記第９のステップＳＴ９で算出された誤差ｄ_b を用い
て、以下に示す誤差行列方程式から誤差ｄ_x を求める。Ａ・ｄ_x ＝ｄ_b ただし、実際の計算は誤差行列方程式Ａ・ｄ_x ＝ｄ_b を
ウェーブレット変換して近似的に解き、得られた近似解
をウェーブレット逆変換することによって誤差ｄ_x をよ
り高速に計算する。すなわち、図１における第５のステ
ップＳＴ５ないし第７のステップＳＴ７と同一の手順に
よってこの誤差行列方程式Ａ・ｄ_x ＝ｄ_b を解いてい
る。

【００３９】次に第１１のステップＳＴ１１において、
上記第１０のステップＳＴ１０で得られた誤差ｄ_x を次
式のように、近似的な物理量分布ｘ^* から差し引いて、
真の解ｘ₁ を求める。ｘ₁ ＝ｘ^* −ｄ_x

【００４０】次に第１２のステップＳＴ１２において、
上記第１１のステップＳＴ１１にて得られた真の解ｘ₁
を近似的な物理量分布ｘ^* に代入して、再度第９のステ
ップＳＴ９ないし第１１のステップＳＴ１１を実行し、
より精度の向上した解ｘ₂ を求める。

【００４１】次に第１３のステップＳＴ１３において、
上記第１１のステップＳＴ１１で得られた真の解ｘ₁ と
第１２のステップＳＴ１２で得られたより精度の向上し
た解ｘ₂ とを比較し、両者の相対誤差があらかじめ定め
られた所定の設定値未満であれば、この一連の計算処理
を終了する。一方、両者の相対誤差が設定値以上であれ
ば、第１２のステップＳＴ１２で得られた解ｘ₂ を真の
解ｘ₁ に代入して、第１２のステップＳＴ１２の処理を
繰り返す。なお、両者の相対誤差の検定に用いられる設
定値としては、０．０１％ないし０．０００１％が現実
的な値であり、この設定値として０．０１％を選んだ場
合には、解の精度は多少ラフなものとなるが解を高速に
求めることができ、０．０００１％を選んだ場合には、
解を得るのに多少時間はかかるが精度の高い正確な解を
求めることができる。

【００４２】ここで、図８は上記繰り返し計算の結果を
示す説明図である。この図８と繰り返し計算による補正
前の結果を示した図５と比較すれば、解の精度は非常に
改善されていることがわかる。また図６と比べれば、時
間がかかるが正確な解が求められる従来の方法で計算し
た場合とほぼ同一精度の解が得られていることも理解で
きる。

【００４３】ここで、上記実施の形態においては境界要
素数が３２、すなわち３２行３２列のシステム行列につ
いて説明したが、実際の推定問題では、さらに行列の要
素数が増大して、推定時間が大幅に増加する。例えば、
１０２４個の要素数の場合を例に取れば、システム行列
Ａの各要素を計算する時間は合計３４秒、ガウスの消去
法で行列方程式Ａ・ｘ＝ｂを解く時間は９７６秒かかっ
た。これに対して、この発明の実施の形態１および２に
よる物理量分布推定方法では、２１秒で行列方程式Ａ・
ｘ＝ｂを解くことができた。なお、解の精度は０．０１
％として収束を判定した。また、計算機にはＤＥＣＡ
ｌｐｈａを用い、計算コードはＦＯＲＴＲＡＮ言語を用
いた。この結果、従来の場合よりも４０倍程度高速に物
理量分布の推定が行えることがわかった。

【００４４】なお、上記各実施の形態では、対象システ
ムの支配方程式が、電位分布に関するラプラス方程式で
あったが、これに限らず、あらゆる境界要素解析の適用
対象が同様に扱える。例えば、電磁界分布、温度分布、
歪み分布、応力分布、変位分布、流速分布など種々の物
理量分布を、この発明の物理量分布推定方法によれば従
来の場合よりも一桁以上も高速に推定することができ
る。その結果、従来は一回の解析しかできなかった時間
で、数十回の計算ができ、例えばシステムの形状を表わ
すパラメータを逐次変化させて、システムの設計者が希
望する物理量分布を発生する最適形状システムを構築す
ることが可能になる。

【００４５】また、図７のフローチャートに示した近似
的な物理量分布ｘ^* を繰り返し計算によって補正する処
理は、最も単純な例を示したに過ぎず、この方法を改良
した種々の方法をとることが可能であり、それぞれ上記
実施の形態と同等の効果を奏する。

【００４６】

【発明の効果】以上のように、請求項１記載の発明によ
れば、ウェーブレット変換によって行列方程式Ａ・ｘ＝
ｂより変換されたＡ_t ・ｘ_t ＝ｂ_t を近似的に解き、得
られた近似解ｘ_t をウェーブレット逆変換して近似的な
物理量分布ｘ^* を求め、それを繰り返し計算によって補
正するように構成したので、密行列に関する行列方程式
を直接解く必要がなくなって、上記行列方程式の計算時
間が大幅に短縮され、物理量分布を高速に推定すること
が可能となるため、各種電磁機器や構造体の最適な形状
設計などが容易になる効果がある。

【００４７】請求項２記載の発明によれば、ウェーブレ
ット変換をドビッシーの４係数型ウェーブレットを用い
て行うように構成したので、より簡単な計算によってウ
ェーブレット変換が可能となる効果がある。

【００４８】請求項３記載の発明によれば、行列方程式
Ａ_t ・ｘ_t ＝ｂ_t を近似的に解く際の閾値を、行列Ａ_t
の各行列要素の絶対値の最大値の０．０１％ないし０．
００００１％の範囲に選定するように構成したので、近
似解ｘ_t をより高速に求めることができる効果がある。

【００４９】請求項４記載の発明によれば、行列方程式
Ａ_t ・ｘ_t ＝ｂ_t を前処理付き双共役勾配法を用いて近
似的に解くように構成したので、近似解ｘ_t の計算速度
を向上させることができる効果がある。

【００５０】請求項５記載の発明によれば、離散化した
ときの境界要素数が２のべき乗となるようにして境界積
分方程式を離散化方程式に変換するように構成したの
で、ウェーブレット変換の際の計算効率を向上させるこ
とができる効果がある。

【００５１】請求項６記載の発明によれば、誤差ｄ_b ＝
Ａ・ｘ^* −ｂにより誤差行列方程式Ａ・ｄ_x ＝ｄ_b に基
づいて誤差ｄ_x を求め、ｘ₁ ＝ｘ^* −ｄ_b より得られた
真の解ｘ₁ を近似的な物理量分布ｘ^* に代入して解ｘ₂
を求め、得られた解ｘ₂ と真の解ｘ₁ の相対誤差が所定
の設定値よりも小さくなるまで、その解ｘ₂ を真の解ｘ
₁ に代入して解ｘ₂ の算出を繰り返すことで、近似的な
物理量分布ｘ^* の補正を行うように構成したので、得ら
れる解の精度を向上させることができる効果がある。

【００５２】請求項７記載の発明によれば、誤差行列方
程式Ａ・ｄ_x ＝ｄ_b から、行列方程式をウェーブレット
変換して近似的に解き、得られた近似解をウェーブレッ
ト逆変換して近似的な物理量分布を求めるという手順に
従って誤差ｄ_x を求めるように構成したので、誤差の補
正を高速に行うことができる効果がある。

【００５３】請求項８記載の発明によれば、真の解ｘ₁
と解ｘ₂ との相対誤差の検定のための所定の設定値を、
０．０１％ないし０．０００１％の範囲で選定するよう
に構成したので、比較的ラフな精度の解を高速で求めた
り、多少の時間をかけて精度の高い解を求めるなど、要
求に応じた選択が可能となる効果がある。

【図面の簡単な説明】

【図１】この発明の実施の形態１による物理量分布推
定方法の処理手順を示すフローチャートである。

【図２】上記実施の形態における対象となるシステム
の一例を示す説明図である。

【図３】上記実施の形態におけるシステム行列Ａを疑
似３次元表示にて示した説明図である。

【図４】上記実施の形態におけるウェーブレット変換
後の行列Ａ_t を疑似３次元表示にて示した説明図であ
る。

【図５】上記実施の形態にて近似的な物理量分布ｘ^*
を求めた結果を示す説明図である。

【図６】上記図５と対比するために従来の方法で計算
した結果を示す説明図である。

【図７】この発明の実施の形態２における近似的な物
理量分布ｘ^* の補正のための繰り返し計算の手順を示す
フローチャートである。

【図８】上記実施の形態における繰り返し計算による
補正結果を示す説明図である。

【図９】従来の境界要素解析方法を示す説明図であ
る。

【符号の説明】

ＳＴ１第１のステップ、ＳＴ２第２のステップ、Ｓ
Ｔ３第３のステップ、ＳＴ４第４のステップ、ＳＴ
５第５のステップ、ＳＴ６第６のステップ、ＳＴ７
第７のステップ、ＳＴ８第８のステップ、ＳＴ９
第９のステップ、ＳＴ１０第１０のステップ、ＳＴ１
１第１１のステップ、ＳＴ１２第１２のステップ、
ＳＴ１３第１３のステップ。

Claims

【特許請求の範囲】

【請求項１】対象システムの支配方程式を与える第１
のステップと、前記支配方程式を境界積分方程式に変換する第２のステ
ップと、前記境界積分方程式を離散化方程式に変換する第３のス
テップと、前記離散化方程式を、前記対象システムに依存したシス
テム行列Ａ、未知の物理量分布を表す未知ベクトルｘ、
および定数ベクトルｂによる行列方程式Ａ・ｘ＝ｂに変
換する第４のステップと、前記行列方程式Ａ・ｘ＝ｂの両辺をウェーブレット変換
して、新たな行列方程式Ａ_t ・ｘ_t ＝ｂ_t を得る第５の
ステップと、前記行列方程式Ａ_t ・ｘ_t ＝ｂ_t における行列Ａ_t の行
列要素中の、絶対値が所定の閾値以下のものをゼロに置
き換え、前記行列方程式Ａ_t ・ｘ_t ＝ｂ_t を近似的に解
いて、近似解ｘ_t を求める第６のステップと、前記近似解ｘ_t をウェーブレット逆変換して、近似的な
物理量分布ｘ^* を求める第７のステップと、前記近似的に求めた物理量分布ｘ^* を繰り返し計算によ
つて補正する第８のステップとを備えた物理量分布推定
方法。
【請求項２】第５のステップにおけるウェーブレット
変換として、ドビッシーの４係数型ウェーブレットを用
いたことを特徴とする請求項１記載の物理量分布推定方
法。
【請求項３】第６のステップにおける所定の閾値を、
行列Ａ_t の各行列要素の絶対値の最大値の０．０１％な
いし０．００００１％に選んだことを特徴とする請求項
１記載の物理量分布推定方法。
【請求項４】第６のステップにおける行列方程式Ａ_t
・ｘ_t ＝ｂ_t を近似的に解く方法として、前処理付き双
共役勾配法を用いたことを特徴とする請求項１記載の物
理量分布推定方法。
【請求項５】境界積分方程式を離散化方程式に変換す
る第３のステップにおいて、離散化したときの境界要素
数を２のべき乗に選んだことを特徴とする請求項１記載
の物理量分布推定方法。
【請求項６】近似的に求めた物理量分布ｘ^* を繰り返
し計算によって補正する第８のステップが、前記近似的な物理量分布ｘ^* とシステム行列Ａおよび定
数ベクトルｂから、誤差ｄ_b ＝Ａ・ｘ^* −ｂを計算する
第９のステップと、誤差行列方程式Ａ・ｄ_x ＝ｄ_b に基づいて、算出された
前記誤差ｄ_b から誤差ｄ_x を求める第１０のステップ
と、得られた誤差ｄ_x を前記近似的な物理量分布ｘ^* から減
算して、真の解ｘ₁ を求める第１１のステップと、得られた真の解ｘ₁ を前記近似的な物理量分布ｘ^* に代
入して再度第９のステップないし第１１のステップを実
行し、解ｘ₂ を求める第１２のステップと、前記真の解ｘ₁ と解ｘ₂ の相対誤差が所定の設定値より
も小さければ計算を終了し、設定値よりも大きければ解
ｘ₂ を真の解ｘ₁ に代入して処理を第１２のステップに
戻す第１３のステップとを備えたことを特徴とする請求
項１記載の物理量分布推定方法。
【請求項７】誤差行列方程式Ａ・ｄ_x ＝ｄ_b から誤差
ｄ_x を求める第１０のステップは、第５のステップない
し第７のステップと同等の手順にて、前記誤差ｄ_x を求
める処理を実行することを特徴とする請求項６記載の物
理量分布推定方法。
【請求項８】第１３のステップにおける所定の設定値
として、０．０１％ないし０．０００１％を選んだこと
を特徴とする請求項６記載の物理量分布推定方法。