JP2005258813A - 情報処理装置及び情報処理方法並びにプログラム - Google Patents

Abstract

【課題】 境界要素法では難しい界面の変化を安定に取り扱うことができるとともに、有限要素法で必要とされる要素分割の数を大幅に減らすことができ、その結果として計算コストを大幅に低減することができる情報処理装置及び情報処理方法並びにプログラムを提供する。
【解決手段】 メモリ６２上に、数値解析の対象となる領域を設定し、当該領域内に同一次元の部分領域を設定する。次に、部分領域を要素分割して離散化された基準点を設定する。次に、外部ポテンシャルが与えられたときの基準点上の近似した外部ポテンシャルと物理量を設定し、それらから定義されるエネルギー関数を最小にするときの物理量を算出する。そして、算出された部分領域内の物理量を領域に対して積分することにより、領域のポテンシャルφを算出する。
【選択図】 図２

Description

本発明は、微分方程式の数値解法において要素分割を用いて近似解を求める情報処理技術に関する。

微分方程式を数値的に解く手法の中心技術として、有限要素法（ＦＥＭ）や境界要素法（ＢＥＭ）が良く知られている。

ＦＥＭは、考察領域全体を有限要素に分割し、重み付け残差法等を用いて、微分方程式を有限要素式に変換する。次に、当該有限要素式の中の未知数を形状関数と離散化された関数で近似し、要素節点での関数値に関する連立一次方程式を得る。そして、この連立一次方程式を解くことによって、領域内節点での値を求めることができる（例えば、非特許文献１参照）。

一方、ＢＥＭは、支配微分方程式をこれと等価な境界積分方程式に変換した後、考察領域の境界表面を要素分割することで離散化し、ＦＥＭと同様に連立一次方程式に帰着させて解く（例えば、非特許文献２参照）。

上記ＦＥＭは、与えられた汎関数に対する「変分＝０」の形の原理である変分原理による汎関数の最小化条件から定式化される。汎関数は、各要素内で定義される要素汎関数の和として与えられる。このため、ＦＥＭは、領域全体の場の平均を扱うことで、非常にロバストで強力な数値解法となっている。また、ＦＥＭは、区分的近似関数を採用するため、近接相互作用のみを考慮すればよい。その結果、連立一次方程式の係数行列はきわめて「疎な状態」となり、その性質を有効に用いた数値解析法が計算機上で実現されている。
O.C.ツィエンキーヴィッツ著、吉識雅夫・山田嘉昭 監訳 "マトリックス有限要素法"、培風館 C.A.ブレビア著、神谷紀生・田中正隆・田中喜久昭 共訳、"境界要素法入門"、培風館

しかしながら、ＦＥＭは、領域全体を有限要素に分割する必要があるため、データ数が膨大となり、入力データを作成するための計算コストがかかるという問題がある。また、取り扱う方程式の元数も大きなものとなるため、特に３次元問題や非線形問題では、より少ない元数の方程式を扱うことが計算機資源の面からみて望ましいと考えられている。

一方、ＢＥＭは、境界表面だけを要素分割し、各節点にデータを入力すればよいので、計算の前処理が容易である。また、ＦＥＭと比較して、次元を一つ下げて取り扱うことができるため、方程式の元数を少なくすることができ、計算コストを下げることができる可能性がある。

しかし、ＢＥＭは、境界という太さ（幅）のない（測度零の）領域での積分を扱うことになるので、特異点の問題等が発生し、数値計算上ロバストではない場合がある。また、境界積分方程式に変換するために遠距離相互作用を取り入れねばならず、離散化して得られる連立一次方程式の係数行列は「密な状態」となる。その結果、ＢＥＭは、元数は少なくて済む一方で、一般にフルマトリックスを取り扱う必要があり、数値計算上のネックと考えられている。

さらに、ＢＥＭには次のような問題点が存在する。例えば、図１１は、界面の変形により要素の大きさが変化する例を説明するための図である。すなわち、図１１に示すように、粒子が成長して界面の形状が変化するような問題である。この場合、境界表面に節点を置くため、粒子が成長すると最初に設定した節点間の長さが変化し、要素の大きさが変化してしまう。その結果、計算の実行に従って、解析精度が変化してしまうことになる。この問題を避けるためには、要素の大きさが一定になるように、節点の数と配置を計算の実行に合わせて調整する必要があり、そのためには新たな節点へのデータの受け渡しや要素数の変化に対応した行列要素の再計算等をする必要が生じてしまう。

また、図１２は、界面の融合により節点の数が変化する例を説明するための図である。すなわち、図１２に示すように、２つの領域がぶつかって一つの領域に合体するような現象を解析しようとすると、融合した近傍で節点を消滅させたり、新たに発生させたりしなければならない。このように、領域の大きさが変化したり、複数の領域が合体したりするような界面の変化を伴う現象を連続して計算しようとする場合、ＢＥＭでは特有の計算手法が必要となり、あまり適当であるとはいえない。

そこで、上述したＦＥＭやＢＥＭのそれぞれの問題を回避するために、ＦＥＭとＢＥＭをハイブリッドに利用するアルゴリズムがこれまでにも提案・研究されている。しかしながら、それらの方法には、特異性の問題等の様々な問題を含んでいる。例えば、Qiya Hu and Dehao Yu, "Solving singularity problems in unbounded domains by coupling of natural BEM and composite grid FEM", Applied Numerical Mathematics 37(2001)127等がある。

本発明は、このような事情を考慮してなされたものであり、境界要素法では難しい界面の変化を安定に取り扱うことができるとともに、有限要素法で必要とされる要素分割の数を大幅に減らすことができ、その結果として計算コストを大幅に低減することができる情報処理装置及び情報処理方法並びにプログラムを提供することを目的とする。

上記課題を解決するために、本発明は、領域を要素分割して数値解析を行う情報処理装置であって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定手段と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定手段と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定手段と、
前記部分領域内の前記基準点に対して、前記領域に与えられる外部ポテンシャルに基づく物理量Ｑ及び外部ポテンシャルＶを設定する第４の設定手段と、
前記基準点に設定された外部ポテンシャルＶと物理量Ｑと密で正則な行列Ｋとで定義されるエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）を最小にする最小化手段と、
前記最小化手段によって最小のエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）が得られたときの物理量Ｑを算出する第１の算出手段と
前記物理量Ｑを前記領域に対して積分することにより、前記領域内における物理量を算出する第２の算出手段と、
前記第２の算出手段による算出結果を出力する出力手段と
を備えることを特徴とする。

また、本発明は、コンピュータに領域を要素分割して数値解析を行わせる情報処理方法であって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定工程と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定工程と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定工程と、
前記部分領域内の前記基準点に対して、前記領域に与えられる外部ポテンシャルに基づく物理量Ｑ及び外部ポテンシャルＶを設定する第４の設定工程と、
前記基準点に設定された外部ポテンシャルＶと物理量Ｑと密で正則な行列Ｋとで定義されるエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）を最小にする最小化工程と、
前記最小化手段によって最小のエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）が得られたときの物理量Ｑを算出する第１の算出工程と
前記物理量Ｑを前記領域に対して積分することにより、前記領域内における物理量を算出する第２の算出工程と、
前記第２の算出手段による算出結果を出力する出力工程と
を有することを特徴とする。

さらに、本発明は、コンピュータに、領域を要素分割して数値解析を行わせるためのプログラムであって、
前記数値解析の対象となる領域をメモリ上に設定する第１の設定手順と、
前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定手順と、
前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定手順と、
前記部分領域内の前記基準点に対して、前記領域に与えられる外部ポテンシャルに基づく物理量Ｑ及び外部ポテンシャルＶを設定する第４の設定手順と、
前記基準点に設定された外部ポテンシャルＶと物理量Ｑと密で正則な行列Ｋとで定義されるエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）を最小にする最小化手順と、
前記最小化手段によって最小のエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）が得られたときの物理量Ｑを算出する第１の算出手順と
前記物理量Ｑを前記領域に対して積分することにより、前記領域内における物理量を算出する第２の算出手順と、
前記第２の算出手段による算出結果を出力する出力手順と
を実行させることを特徴とする。

本発明によれば、境界要素法では難しい界面の変化を安定に取り扱うことができるとともに、有限要素法で必要とされる要素分割の数を大幅に減らすことができ、その結果として計算コストを大幅に低減することができる。

＜本発明の概要＞
本発明は、幅の無い（すなわち、測度零の）領域での数値解析である境界要素法的手法に、人為的に幅（太さ）を導入することで、問題をエネルギー最小問題として捉えて定式化し、安定で効率的な数値解法を与えるものである。

以下、本発明の概要について具体的に説明する。

図１３は、領域Ｒ中に与えられた外部ポテンシャルｖ（ｘ）の様子を示す図である。図１３に示すように、領域Ｒ内に外部ポテンシャルｖ（ｘ）が与えられる場合、その電荷分布をｑ（ｘ）、その電荷分布によって作られるポテンシャルをφ（ｘ）とすると、エネルギーＥ₁は式（１）に示すようなエネルギー積分の和によって表すことができる。

そして、式（１）のエネルギー積分の和が最小となるところで、正しいポテンシャル分布φ（ｘ）が得られることが知られている。ここで、右辺第１項は静電エネルギー、第２項はポテンシャルエネルギーである（最小作用の原理）。
式（１）に示すエネルギーＥ₁が最小値をとる場合、外部ポテンシャル分布φ（ｘ）や電荷分布ｑ（ｘ）を微小に変化させてもエネルギーＥの変分が零という条件より、次の２式を得る。

ここで、Δは、線形偏微分作用素であるラプラス演算子を表す。また、式（２）はポアソン方程式であり、式（３）はつりあいの式である。上記２式より、外部ポテンシャルｖ（ｘ）を与えることで、ポテンシャル分布φ（ｘ）及び電荷分布ｑ（ｘ）を求めることができる。
一方、領域Ｒ上で定義された線形偏微分作用素Ｌと２点デルタ関数δ（ｘ−ｙ）からなる式（４）で示される方程式は、基本解として２点関数ｕ^*（ｘ，ｙ）を持つことが知られている。

ここで、ｘ及びｙは、領域Ｒ内の座標を表す。
この基本解を用いた場合であって、領域Ｒ上で定義される物理量φ（ｘ）及びｑ（ｘ）を支配する偏微分方程式が式（５）で与えられる場合、物理量φ（ｘ）は、ｕ*（ｘ，ｙ）を積分核として式（６）に示すように表すことができることが知られている。

尚、線形偏微分作用素Ｌとしてはラプラス型、ヘルムホルツ型がよく知られているが、本実施形態ではこれらに限られるものではない。ラプラス型及びヘルムホルツ型には、それぞれに対応して基本解ｕ^*（ｘ，ｙ）が求められている。尚、基本解については前出の非特許文献１にその詳細が記されている。
線形偏微分作用素Ｌをラプラス演算子としたとき、式（１）で与えられるエネルギーＥ₁は、上記基本解ｕ*（ｘ，ｙ）を積分核として、式（７）で示すようなエネルギーＥ₂のように別表現されることが知られている（例えば、V.N.Popov著、"Functional integrals and collective excitations" 、ケンブリッジ大学出版、５９頁参照）。

上記エネルギーＥ₂の電荷分布ｑ（ｘ）の微小変化に対する変分が零という条件（最小作用の原理）から式（８）が得られる。

ここで、式（６）を用いると、式（８）は式（９）に示すようになる。

また、式（５）の線形偏微分作用素Ｌがラプラス演算子Δであることから、式（１０）を得る。

すなわち、式（１０）及び式（９）は、それぞれ式（２）及び式（３）と同じものとなる。
以上説明したように、式（１）で表されるエネルギーＥ₁の最小値を求める問題は、式（７）で表されるエネルギーＥ₂の最小化問題と等価であることが示された。そこで、本発明では、以下に示す実施形態においてより具体的に、式（７）で表されるエネルギーＥ₂を離散化して、計算機上で安定で効率的な数値解析手法を与えることとする。

＜第１の実施形態＞
［情報処理装置の構成］
以下、図面を参照して、本発明の第１の実施形態に係る情報処理装置及び当該装置を用いた情報処理方法についてより具体的に説明する。図６は、本発明の第１の実施形態に係る微分方程式の数値解法において要素分割を用いて近似解を求める情報処理装置の構成を示すブロック図である。図６に示すように、本実施形態に係る情報処理装置は、各種演算を実行するＣＰＵ６１やメモリ６２等から構成される演算処理部１、ハードディスク６３やフレキシブルディスク６４からなる記憶部２、キーボード６５やマウス６６等からなる入力部３、そしてディスプレイ６７等からなる出力表示部４からで構成されている。

［情報処理装置の基本的動作］
図２は、本発明の第１の実施形態に係る情報処理装置の動作手順を説明するためのフローチャートである。本実施形態に係る情報処理装置では、メモリ６２上に図２に示したフローチャートに沿って設計、コード化されたプログラムをロードし、必要な計算領域を確保し、適当な方法で入力された領域形状や格子点座標及び物理量に基づいて所定の演算処理を行う。そして、その結果得られた領域上での物理量の値をハードディスク６３等に書き込んで記憶させたり、ディスプレイ６７上に表示させる。

まず、図１に示すようにポテンシャルφを演算する対象となる領域Ａを設定する（ステップＳ１）。図１は、本発明の第１の実施形態を説明するための領域Ａとその部分領域ａとを説明するための概念図である。尚、本実施形態に係る領域Ａの次元は２次元としているが、本発明の適用は２次元に限られるものではなく、３次元であってもよい。

次に、領域Ａの部分領域ａを設定する（ステップＳ２）。ここで、部分領域ａは、領域Ａと同一次元の領域である。前述したように、ＦＥＭを用いて領域Ａの物理量を計算する場合は、通常、領域Ａ内の全格子点を未知数として計算しなければならない。また、ＢＥＭを用いて領域Ａの物理量を計算する場合は、部分領域ａに相当するものは次元が一つ少ないものを用いる。しかし、本実施形態では、領域Ａが２次元の場合は部分領域ａは２次元（面）であり、領域Ａが３次元（（立体）の場合は部分領域ａも立体であるとする。また、部分領域ａの形状は任意である。

次に、部分領域ａ上で、式（７）に示すエネルギー積分を離散化するため、部分領域ａ上に離散化された格子点を設定する（ステップＳ３）。また、各格子点には番号を付け、それらをｉ＝１，２，…，Ｎとする。ここで、Ｎは、領域Ａ上に設定された格子点の総数である。

次に、物理量ｑを近似表現した格子点上の値Ｑ_i (ｉ＝１，２，…，Ｎ）及び外部ポテンシャルｖを近似表現した格子点上の値Ｖ_i（ｉ＝１，２，…，Ｎ）を設定する（ステップＳ４）。尚、本実施形態は、便宜上格子点を用いるが、離散点は格子点だけに限られるものではなく、格子辺上で定義してもよいし、また格子胞上等で定義してもよい。

上記の定義を用いることで、式（７）に示されるエネルギー積分は、式（１１）に示すような格子上で定義される物理量の和の形に書き換えることができる（ステップＳ５）。

ここで、密で正則なマトリックスＫ_ij （ｉ，ｊ＝１，２，…，Ｎ）を基本解ｕ^*（ｘ，ｙ）を積分核として、式（１２）に示すように定義する。

また、同様に、物理量ｑ及び外部ポテンシャルｖの格子点上の値Ｑ_i及びＶ_iを、それぞれ式（１３）、（１４）に示すように定義する。

ここで、ｆ_i（ｘ）及びｇ_j（ｙ）は、それぞれ部分領域ａ上に定義された非零の領域（すなわち、境界）が、格子点ｉ又はｊの周りに局在した関数であり、積分は部分領域ａ上で計算される。また、上記の積分は、解析的又は数値的に行ってもよく、数値積分を実行する際は、ｆ_i（ｘ）及びｇ_j（ｙ）を適当な内挿関数とする。尚、このマトリックスＫ_ijは、その大部分の要素（例えば、９０パーセント以上の要素）が非零の値を持つ行列要素からなることを特徴とする。
以上の定式化より、部分領域ａ内のＮ個の格子上での値｛Ｖ_i｝を与えることにより、式（１１）で与えられたエネルギーＥ₂を最小にするような電荷分布ｑの近似値の組｛Ｑ_i｝を求めることができる（ステップＳ６）。

このようにして求められた部分領域ａ上の格子上での値Ｑ_iから、部分領域ａの領域Ａでの補集合上で、座標ｘ^Pを有する任意の点Ｐ上でのポテンシャルφ（ｘ^P）は、式（１５）に示すように近似される。

尚、和は部分領域ａのＮ個の格子点にわたってとるものとする。尚、式（１６）において、ｇ_j（ｙ）は前述した関数であり、積分は部分領域ａ上で行うものとする。
以上の手順により、領域Ａの部分領域ａの離散化された格子点に設定された外部ポテンシャル値Ｖ_iから、部分領域ａの領域Ａでの補集合上の任意の点におけるポテンシャルφを求めることができる（ステップＳ７）。

＜第２の実施形態＞
次に、上述した第１の実施形態に係る情報処理装置を用いて、例えば、電磁場のような物理場のポテンシャルを求めるプログラムの流れについて説明する。図７は、本発明の第２の実施形態において計算対象となる領域Ｂ及びその部分領域ｂを説明するための図である。本実施形態では、図７に示すように、領域Ｂを２次元で正方な領域とし、部分領域ｂを同心円で囲まれた領域とした。尚、領域Ｂ及び部分領域ｂの形状は、これだけに限定されることはなく、任意の次元の領域を設定することが可能である。尚、本実施形態における情報処理装置は、前述した第１の実施形態における情報処理装置と同一であり、その処理手順を示すフローチャートは、前述した第１の実施形態で用いられた図２のフローチャートと同一である。すなわち、ステップＳ１、Ｓ２において、図７に示す領域Ｂ及びその部分領域Ｂが設定される。

まず、正方領域ＢをＭ×Ｍの等間隔に分割し、メモリ６２上に（Ｍ＋１）×（Ｍ＋１）点分の格納スペースを作り、そこに各格子点の座標を収納する。また、ｖやｑの離散版であるＵやＱを（Ｍ＋１）×（Ｍ＋１）ベクトルとしてメモリ６２上にとる。次に、各格子点についてその座標が部分領域ｂ上に存在するかどうかの判定を行い、存在する場合のみ選択して、新たなメモリ領域に登録する。その結果、部分領域ｂ上に存在するＮ個の格子点座標ｘ_iからなる、メモリ空間で連続したベクトルＸを得ることができる。上記処理は、前述のステップＳ３に相当する。尚、正方領域ではない領域を分割する場合は、当該領域をＭ×Ｎに分割し、ＶやＱを（Ｍ＋１）×（Ｎ＋１）ベクトルとしてメモリ上にとればよい。

次に、メモリ空間にＮ個のデータを格納するリストベクトルＬ（ｉ）を用意し、リストベクトルＬ（ｉ）が（Ｍ＋１）×（Ｍ＋１）格子内における部分領域ｂのＮ個のデータを指すようにする。以下、簡単のために、Ｖ（Ｌ（ｉ））をＵ_i、Ｑ（Ｌ（ｉ））をＱ_iと記す。但し、ｉ＝１，２，…，Ｎである。また同様に、行列要素Ｋ_ijを格納するためのマトリックスＫをメモリ６２上に用意する。

図８は、第２の実施形態において部分領域ｂの格子点上に設定した値Ｖ_iを示す鳥瞰図である。図８では、図を見やすくするために、部分領域ｂ外に対しては一定の値が表示されているが、実際に値は設定されていない。すなわち、この物理量の設定は前述のステップＳ４に相当する。

次に、線形微分作用素Ｌとして２次元ラプラス演算子を選択すると、式（４）から基本解として、次式が得られる。

次に、マトリックＫの成分を求めるため、式（１２）における関数ｆ_i（ｘ）及びｇ_j（ｙ）を次式のようにおく。

尚、本実施形態では、関数ｆ_i（ｘ）及びｇ_j（ｙ）をδ関数としているが、これに限られるものではなく、部分領域ｂで定義された非零の関数であれば良い。上式において、ｘⁱ及びｙ^jは、それぞれ部分領域ｂ上のｉ番目及びｊ番目の格子点座標である。また、式（１７）、式（１８）及び式（１９）より式（１２）を用いて、Ｋ_ijを次のように計算する。

すなわち、このＫ_ijの計算は、前述のステップＳ５に相当する。上記計算をするために、ベクトルＸのｉ番目とｊ番目を参照して格子点座標を取り出し、計算された行列要素をマトリックスＫの（ｉ，ｊ）番地に格納する。
次に、ベクトルＶ、マトリックスＫに対して、式（１１）のエネルギーＥ2を最小化するようなＱiの組を最急勾配法により求めた。エネルギー最小化問題は、良く知られているように、初期値に依存して局所最小値（極小値）に落ち着く場合もある。本実施形態における最小もこの意味で用いられる。図９は、第２の実施形態で計算された部分領域ｂ上で定義された値Ｑ_iを示す鳥瞰図である。図９においても、図８と同様に、部分領域ｂ外に対しては一定の値を入れ、図を見やすくしている。すなわち、このベクトルＱの演算は、前述したステップＳ６に相当する。

次に、部分領域ｂ上の各格子点に対して求められたＱ_iを用いて、式（１５）及び式（１６）から、部分領域ｂの領域Ｂでの補集合上の座標ｘ^Pを有する任意の点Ｐに対してポテンシャルφ（ｘ^P）を求め、それを最初に設定した領域ｂ上の格子点上での値Ｖ_iと合わせて表示する。図１０は、第２の実施形態において領域Ｂ上で計算されるポテンシャルφを表す鳥瞰図である。すなわち、ポテンシャルφの取得が、前述したステップＳ７の処理に相当する。

以上のように、計算機上のメモリ空間に設定したベクトルやマトリックスを演算することで、部分領域ｂ上の格子点上の値Ｖ_iから領域Ｂ上のポテンシャルφを近似的に求めることができる。

＜第３の実施形態＞
図３は、本発明の第３の実施形態に係る情報処理装置による処理内容を説明するための図である。本実施形態における情報処理装置の構成及びその動作は前述したものと同様である。図３において、領域Ｃの部分領域ｃは２つの分離された領域からなる。

本実施形態においても、前述した第２の実施形態と同様に部分領域ｃ上で離散化された格子点を設定する。そして、設定された格子点数をＮとして、各格子点上に物理量を設定することにより、第２の実施形態と同様の処理手順に従って、領域Ｃ上の任意の点におけるポテンシャルφを求めることができる。

尚、部分領域ｃが３つ以上の分離された領域からなる場合であっても、上記手順が同様に成り立つことは明らかである。

＜第４の実施形態＞
図４は、本発明の第４の実施形態に係る情報処理装置による処理内容を説明するための図である。本実施形態における情報処理装置の構成及びその動作は前述したものと同様である。図４に示すように、本実施形態では、領域Ｄは、その部分領域ｄの境界の一部が一致している。

本実施形態においても、前述した第２の実施形態と同様に、まず、部分領域ｄ上で離散化された格子点を設定する。そして、設定された格子点数をＮとして、各格子点上に物理量を設定することにより、第２の実施形態と同様の処理手順に従って、領域Ｄ上の任意の点におけるポテンシャルφを求めることができる。

尚、部分領域ｄの境界が領域Ｄの境界全部と一致している場合でも、上記手順が同様に成り立つことは明らかである。

＜第５の実施形態＞
図５は、本発明の第５の実施形態に係る情報処理装置による処理内容を説明するための図である。本実施形態における情報処理装置の構成及びその動作は前述したものと同様である。図５に示すように、本実施形態では、部分領域ｅ上に非正方格子状に格子点を設定した場合を示す。この場合も、前述した第２の実施形態と同様に、まず、部分領域ｅ上の適当な位置にＮ個の離散化された格子点を設定する。そして、このときの各格子点上に物理量を設定することにより、第２の実施形態と同様の処理手順に従って、領域Ｅ上の任意の点におけるポテンシャルφを求めることができる。

また、格子点の離散化を正方格子ではなく、三角格子や任意の非正方規則格子を用いて行っても同様に処理することができる。さらに、適当な擬似乱数を用いて格子点位置を設定した場合であっても、上記手順が同様に成り立つことは明らかである。

＜第６の実施形態＞
前述した第１の実施形態では、図２に示すように、まず領域Ａを設定し（ステップＳ１）、次に領域Ａと同一次元の部分領域ａを設定した（ステップＳ２）。本実施形態では、幅のない（すなわち、測度零の）境界Γでの数値解解析である境界要素法的手法に対して、人為的に幅（太さ）を導入することで、問題をエネルギー最小問題として捉えて定式化する手法に関した本発明に対し、さらに、導入する外部ポテンシャルｖ（ｘ）を決めるために、境界Γからの距離関数を用いることで、安定で効率的な数値解法を与えるものである。

以下、本実施形態について具体的に説明する。

図１４は、本発明の第６の実施形態を説明するための領域Ｇとその部分領域ｇ及び幅のない境界Γを説明するための概念図である。本実施形態では、図１４に示すような幅のない（測度零の）境界Γからの距離関数φをレベルセット法により定義する。レベルセット法は、ある曲面を空間の座標（ｘ）の関数がゼロを持つ点の集合と見なし、そこからの距離が一定の等高面を計算する方法である（例えば、Level Set Methods and First Marching Methods, J.A.Sethian, Cambridge, 1999を参照）。レベルセット法では、上記距離関数φ（ｘ）を、次式で示す微分方程式の定常解として与える。

ここで、ｆ₀は、境界Γで零の値、境界Γの内側で負の値、境界Γの外側で正の値をとる関数である。また、φ（ｘ） が与えられた時、領域Ｇ内での外部ポテンシャルｖ（ｘ）を次式に示すように与える。

このように、領域Ｆ内で外部ポテンシャルｖ（ｘ）を与えることにより、測度零の領域Γから、等距離の位置に同電位の外部ポテンシャルが形成されることになる。
そして、以下では、＜本発明の概要＞において説明した手順に従って、上記外部ポテンシャルが与えられた時の電荷分布をｑ（ｘ）、その電荷分布によって作られるポテンシャルをφ（ｘ）としたときの式（３）で与えられるエネルギー積分Ｅ₁に基づいて、最小値を求める問題を式（９）で表されるエネルギー積分Ｅ₂の最小化問題と等価であるとして、エネルギー積分Ｅ₂を離散化して、計算機上で安定で効率的な数値解析手法を与えることができる。

尚、本実施形態における情報処理装置の構成及びその動作は第１の実施形態で用いられた情報処理装置と同様である。また、図１５は、本発明の第６の実施形態に係る情報処理装置の動作手順を説明するためのフローチャートである。ここで、領域Ｇの設定（ステップＳ２１）は第１の実施形態でのステップＳ１の処理と同様である。

本実施形態では、領域Ｇの部分領域中に測度零の領域、つまり、領域Ｇの次元より一次元低い境界Γを設定する（ステップＳ２２）。一般に、境界Γは、解く問題が複数の領域によって構成されている場合、これらの領域の境界に一致させることが望ましい。すなわち、例えば、真空中に置かれた金属におけるポテンシャル問題の場合、金属と真空の境界に境界Γを一致させることが望ましいが、これに限るものではない。

次に、領域Ｇの部分領域であるｇを設定する（ステップＳ２３）。尚、部分領域ｇは、第１の実施形態と同様に領域Ｇと同一次元の領域である。

また、第１の実施形態と同様に、部分領域ｇ上でのエネルギー式を離散化するため、部分領域ｇ上に離散化された格子点を設定する（ステップＳ４）。そして、各格子点には番号を付ける（ｉ＝１，２，…，Ｎ）。ここで、Ｎは部分領域ｆ上に設定された格子点の総数である。

次に、境界Γを、空間の座標（ｘ）の関数が、ゼロを持つ点の集合と見なし、そこからの距離が一定の等高面を表す距離関数φを前記のレベルセット法により計算する。このためには、領域Ｆ内で、境界Γの内側の格子点では負、境界Γの外側の格子点では正となるように、前述のｆ₀を近似表現した格子点上の値ｆ_0i 、（但し、ｉ＝１，２，…，Ｍ。尚、Ｍは，領域Ｆ内の格子点の総数。）、及び、φ_i（ｉ＝１，２，…，Ｍ）を設定し（ステップＳ２５）、式（２１）を数値的に解く。

次に、上記数値解φ_iを用いて、部分領域ｆ内の格子点に対して、外部ポテンシャルｖを与え、式（２２）によって、Ｖ_i（ｉ＝１，２，…，Ｎ）を設定する（ステップＳ２６）。式（２２）におけるｆの具体形としては、問題によって種々の関数を用いることができる。

次に、物理量ｑを近似表現した部分領域ｆ内の格子点上の値Ｑ_i（ｉ＝１，２，…，Ｎ）を設定する。以上の説明では、格子点上に値を設定した例を示したが、これらは、格子点に限られるものではなく、格子辺上で定義してもよいし、また格子胞上等で定義してもよい。

上記の定義を用いることで、エネルギーの積分は、前述した式（１１）で示されるような格子上で定義される物理量の和の形に書き換えることができる（ステップＳ２７）。以下、前述した手順に従って、式（１１）で示されるエネルギーを最小にするような電荷分布の近似値の組｛Ｑ_i｝を求める（ステップＳ２８）。そして、同様にして、部分領域ｇの領域Ｇでの補集合上の任意の点におけるポテンシャルφを求める（ステップＳ２９）。

＜第７の実施形態＞
図１６は、本発明の第７の実施形態で用いられる領域Ｈと境界Γ及び部分領域ｈを説明するための図である。本実施形態では、領域Ｈは２次元で正方の領域であり、部分領域ｈは同心円で囲まれた領域とする。尚、前述と同様に、領域等の形状は図１６に示す形状の限りではなく、任意の次元の領域を設定することが可能である。

まず、正方領域ＨをＭ×Ｍの等間隔に分割し、計算機メモリ上に（Ｍ＋１）×（Ｍ＋１）点分の格納スペースを作り、そこに各格子点の座標を収納する。次に、Ｆ₀及びφを（Ｍ＋１）×（Ｍ＋１）ベクトルとしてメモリ６２上に取る。その後、境界Γよりも内側で負、外側で正となるようにＦ₀を設定し、式（２１）をこの離散化したＦ₀及びφについて解くことにより、各格子点でのレベルセット関数φの値を数値的に求めることができる。尚、本実施形態では、式（２１）は風上差分法を用いて解くことにより定常解を求める。

以下、前述した第２の実施形態と同様に、領域Ｈ上の物理量と部分領域ｈ上の物理量との関係を定義する偏微分方程式を基本解を用いて離散化する。但し、本実施形態では、さらに、レベルセット関数φから式（２２）によって、部分領域ｈ内の格子点上にＶ_iを決定する。尚、本実施形態では、式（２２）における関数ｆとして次式を用いる。

尚、ａは定数であり、格子間隔や、問題に応じて任意に決定することができる。前述した図８は、上記Ｖ_iを３次元表示したものでもある。
そして、以下では同様に、ベクトルＶ、マトリックスＫに対して、式（１１）で示すエネルギーＥ₂を最小化するようなＱ_iの組を最急降下法を用いて求め、図９に示すような結果を得る。さらに、部分領域ｇ上の各格子点に対して求められたＱ_iを用いて、式（１５）及び式（１６）から、部分領域ｈの領域Ｈでの補集合上の座標ｘ^Pを有する任意の点Ｐに対してポテンシャルφ（ｘ^P）を求め、それを最初に設定した領域ｇ上の格子点上での値Ｖ_iと合わせて表示し、図１０に示すポテンシャルφを得ることができる。

以上のように、界面からの距離関数を求めて、当該距離関数よりポテンシャルを定義してエネルギー関数の最小化問題を解くことによって、界面が複雑な形状を持つ場合においても、部分領域ｂ上の格子点上の値Ｖ_iから領域Ｂ上のポテンシャルφを近似的に、効率的に求めることができる。

＜その他の実施形態＞
以上、実施形態例を詳述したが、本発明は、例えば、システム、装置、方法、プログラム若しくは記憶媒体等としての実施態様をとることが可能であり、具体的には、複数の機器から構成されるシステムに適用しても良いし、また、一つの機器からなる装置に適用しても良い。

尚、本発明は、前述した実施形態の機能を実現するソフトウェアのプログラム（実施形態では図に示すフローチャートに対応したプログラム）を、システムあるいは装置に直接あるいは遠隔から供給し、そのシステムあるいは装置のコンピュータが該供給されたプログラムコードを読み出して実行することによっても達成される場合を含む。

従って、本発明の機能処理をコンピュータで実現するために、該コンピュータにインストールされるプログラムコード自体も本発明を実現するものである。つまり、本発明は、本発明の機能処理を実現するためのコンピュータプログラム自体も含まれる。

その場合、プログラムの機能を有していれば、オブジェクトコード、インタプリタにより実行されるプログラム、ＯＳに供給するスクリプトデータ等の形態であっても良い。

プログラムを供給するための記録媒体としては、例えば、フロッピー（登録商標）ディスク、ハードディスク、光ディスク、光磁気ディスク、ＭＯ、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、磁気テープ、不揮発性のメモリカード、ＲＯＭ、ＤＶＤ（ＤＶＤ−ＲＯＭ，ＤＶＤ−Ｒ）などがある。

その他、プログラムの供給方法としては、クライアントコンピュータのブラウザを用いてインターネットのホームページに接続し、該ホームページから本発明のコンピュータプログラムそのもの、もしくは圧縮され自動インストール機能を含むファイルをハードディスク等の記録媒体にダウンロードすることによっても供給できる。また、本発明のプログラムを構成するプログラムコードを複数のファイルに分割し、それぞれのファイルを異なるホームページからダウンロードすることによっても実現可能である。つまり、本発明の機能処理をコンピュータで実現するためのプログラムファイルを複数のユーザに対してダウンロードさせるＷＷＷサーバも、本発明に含まれるものである。

また、本発明のプログラムを暗号化してＣＤ−ＲＯＭ等の記憶媒体に格納してユーザに配布し、所定の条件をクリアしたユーザに対し、インターネットを介してホームページから暗号化を解く鍵情報をダウンロードさせ、その鍵情報を使用することにより暗号化されたプログラムを実行してコンピュータにインストールさせて実現することも可能である。

また、コンピュータが、読み出したプログラムを実行することによって、前述した実施形態の機能が実現される他、そのプログラムの指示に基づき、コンピュータ上で稼動しているＯＳなどが、実際の処理の一部または全部を行い、その処理によっても前述した実施形態の機能が実現され得る。

さらに、記録媒体から読み出されたプログラムが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれた後、そのプログラムの指示に基づき、その機能拡張ボードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によっても前述した実施形態の機能が実現される。

本発明の第１の実施形態を説明するための領域Ａとその部分領域ａとを説明するための概念図である。 本発明の第１の実施形態に係る情報処理装置の動作手順を説明するためのフローチャートである。 本発明の第３の実施形態に係る情報処理装置による処理内容を説明するための図である。 本発明の第４の実施形態に係る情報処理装置による処理内容を説明するための図である。 本発明の第５の実施形態に係る情報処理装置による処理内容を説明するための図である。 本発明の第１の実施形態に係る微分方程式の数値解法において要素分割を用いて近似解を求める情報処理装置の構成を示すブロック図である。 本発明の第２の実施形態において計算対象となる領域Ｂ及びその部分領域ｂを説明するための図である。 第２の実施形態において部分領域ｂの格子点上に設定した値Ｖ_iを示す鳥瞰図である。 第２の実施形態で計算された部分領域ｂ上で定義された値Ｑ_iを示す鳥瞰図である。 第２の実施形態において領域Ｂ上で計算されるポテンシャルφを表す鳥瞰図である。 界面の変形により要素の大きさが変化する例を説明するための図である。 界面の融合により節点の数が変化する例を説明するための図である。 領域Ｒ中に与えられた外部ポテンシャルｖ（ｘ）の様子を示す図である。 本発明の第６の実施形態を説明するための領域Ｇとその部分領域ｇ及び幅のない境界Γを説明するための概念図である。 本発明の第６の実施形態に係る情報処理装置の動作手順を説明するためのフローチャートである。 本発明の第７の実施形態で用いられる領域Ｈと境界Γ及び部分領域ｈを説明するための図である。

符号の説明

１ 演算処理部
２ 記憶部
３ 入力部
４ 出力表示部
６１ ＣＰＵ
６２ メモリ
６３ ハードディスク
６４ フレキシブルディスク
６５ キーボード
６６ マウス
６７ ディスプレイ

Claims (15)

1. 領域を要素分割して数値解析を行う情報処理装置であって、
   前記数値解析の対象となる領域をメモリ上に設定する第１の設定手段と、
   前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定手段と、
   前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定手段と、
   前記部分領域内の前記基準点に対して、前記領域に与えられる外部ポテンシャルに基づく物理量Ｑ及び外部ポテンシャルＶを設定する第４の設定手段と、
   前記基準点に設定された外部ポテンシャルＶと物理量Ｑと密で正則な行列Ｋとで定義されるエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）を最小にする最小化手段と、
   前記最小化手段によって最小のエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）が得られたときの物理量Ｑを算出する第１の算出手段と
   前記物理量Ｑを前記領域に対して積分することにより、前記領域内における物理量を算出する第２の算出手段と、
   前記第２の算出手段による算出結果を出力する出力手段と
   を備えることを特徴とする情報処理装置。
2. 前記物理量Ｑが、２次元又は３次元のベクトル量であることを特徴とする請求項１に記載の情報処理装置。
3. 前記第３の設定手段が、
   前記領域をＭ×Ｎに分割する分割手段と、
   前記領域内における（Ｍ＋１）×（Ｎ＋１）点分の格子点のうち、前記部分領域に属する格子点を前記基準点として選択する選択手段と
   を備えることを特徴とする請求項１又は２に記載の情報処理装置。
4. 前記密で正則な行列Ｋは、線形偏微分作用素Ｌと、前記部分領域内の２点であるｘ及びｙから成るデルタ関数δ（ｘ−ｙ）とを用いて、
   Ｌ・ｕ^*＋δ（ｘ−ｙ）＝０
   を満足する２点関数ｕ^*を積分核として前記部分領域上で積分することにより、該部分領域内の前記基準点上で定められることを特徴とする請求項１から３までに記載のいずれか１項に記載の情報処理装置。
5. 前記エネルギー関数Ｅは、
   Ｅ＝Ｑ^TＫＱ＋Ｑ^TＶ
   で与えられることを特徴とする請求項１から４までのいずれか１項に記載の情報処理装置。
6. 前記密で正則な行列Ｋが、前記部分領域上の関数の集合ｆ_i（ｘ）及びｇ_j（ｙ）を用いて、
   Ｋ_ij＝∫ｕ^*（ｘ，ｙ）ｆ_i（ｘ）ｇ_j（ｙ）ｄｘｄｙ
   で定義されることを特徴とする請求項１から５までのいずれか１項に記載の情報処理装置。
7. 前記密で正則な行列Ｋは、大部分の要素が零でない値を有する行列要素からなることを特徴とする請求項１から６までのいずれか１項に記載の情報処理装置。
8. 前記線形偏微分作用素Ｌが、ラプラス演算子Δであることを特徴とする請求項４に記載の情報処理装置。
9. 前記部分領域内に属する前記格子点に通し番号を付与する番号付与手段をさらに備え、
   前記第１の算出手段は、前記通し番号が付与された前記部分領域内に属する格子点上の物理量の近似解を算出する
   ことを特徴とする請求項３に記載の情報処理装置。
10. 前記分割手段が、前記領域をＭ×Ｎの正方格子に分割することを特徴とする請求項３に記載の情報処理装置。
11. 前記領域内に測度零の境界Γを設定する境界設定手段をさらに備え、
    前記第２の設定手段は、前記境界Γを内包するように前記部分領域を設定し、
    前記第４の設定手段は、前記測度零の境界Γからの距離関数を用いて前記外部ポテンシャルＶを設定する
    ことを特徴とする請求項１から１０までのいずれか１項に記載の情報処理装置。
12. 前記距離関数は、前記測度零の境界Γ上で零、該境界Γよりも内側で負、該境界Γよりも外側で正の値をとる関数ｆ₀を用いて、
    ∂φ／∂τ＝ｓｉｇｎ（ｆ0）（１−｜∇φ（τ）｜）
    の定常解として与えられるφに基づいて計算されることを特徴とする請求項１１に記載の情報処理装置。
13. コンピュータに領域を要素分割して数値解析を行わせる情報処理方法であって、
    前記数値解析の対象となる領域をメモリ上に設定する第１の設定工程と、
    前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定工程と、
    前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定工程と、
    前記部分領域内の前記基準点に対して、前記領域に与えられる外部ポテンシャルに基づく物理量Ｑ及び外部ポテンシャルＶを設定する第４の設定工程と、
    前記基準点に設定された外部ポテンシャルＶと物理量Ｑと密で正則な行列Ｋとで定義されるエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）を最小にする最小化工程と、
    前記最小化手段によって最小のエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）が得られたときの物理量Ｑを算出する第１の算出工程と
    前記物理量Ｑを前記領域に対して積分することにより、前記領域内における物理量を算出する第２の算出工程と、
    前記第２の算出手段による算出結果を出力する出力工程と
    を有することを特徴とする情報処理方法。
14. コンピュータに、領域を要素分割して数値解析を行わせるためのプログラムであって、
    前記数値解析の対象となる領域をメモリ上に設定する第１の設定手順と、
    前記領域内に該領域と同一次元の部分領域を前記メモリ上に設定する第２の設定手順と、
    前記部分領域を要素分割して離散化された基準点を前記メモリ上に設定する第３の設定手順と、
    前記部分領域内の前記基準点に対して、前記領域に与えられる外部ポテンシャルに基づく物理量Ｑ及び外部ポテンシャルＶを設定する第４の設定手順と、
    前記基準点に設定された外部ポテンシャルＶと物理量Ｑと密で正則な行列Ｋとで定義されるエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）を最小にする最小化手順と、
    前記最小化手段によって最小のエネルギー関数Ｅ（Ｖ，Ｑ，Ｋ）が得られたときの物理量Ｑを算出する第１の算出手順と
    前記物理量Ｑを前記領域に対して積分することにより、前記領域内における物理量を算出する第２の算出手順と、
    前記第２の算出手段による算出結果を出力する出力手順と
    を実行させるためのプログラム。
15. 請求項１４に記載のプログラムを格納したことを特徴とするコンピュータ読み取り可能な記録媒体。

Images