JP2005267214A - 有限要素法による制御プログラムおよび記憶媒体 - Google Patents

Abstract

【課題】 従属的節点を有する構造格子に対する有限要素法解析を行うにあたり、基底の台に包含関係が存在する場合において、加法的多重レベル法を利用しても、精度よく、簡便に、かつ高速に解析処理することが可能な情報処理技術を提供することを目的とする。
【解決手段】 有限要素法によるプログラム処理をコンピュータによって実現させるための制御プログラムであって、力学的自由度によって張られる空間をＷ、独立節点と従属的節点を１次独立基底として張られる空間をＶとして、ＶへのＷの埋め込み矩形行列をＳとすると、Ｖにおける共役勾配法の前処理因子ＡをＳ^tＡＳに置き換えることを特徴とする制御プログラム。
【選択図】 図１２

Description

本発明は有限要素法による解析処理に関するものである。

有限要素法の基本的な考え方は、微分作用素の定義域である関数空間Ｆを実有限次元ヒルベルト空間Ｈで近似することにより、境界値問題として定まった偏微分方程式の解を求めるというＲｉｔｚの方法の上に立っている。

当該方法において、有限次元ヒルベルト空間Ｈは有限次元ベクトル空間Ｖと内積＜＞の組(Ｖ，＜＞）によって定められる。このとき、有限次元ベクトル空間Ｖから関数空間Ｆへの線形写像が単射となるように構成されるため、関数空間Ｆへの様々な作用（例えば微分作用）はその自然な制限として有限次元ベクトル空間Ｖ上の作用として定義することができる。

有限要素法の特徴は、有限次元ヒルベルト空間Ｈの基底関数として極めて狭い台（非ゼロ集合の閉包）をもつ関数を利用することにある。そのため、有限次元ベクトル空間Ｖに対する有限次元ヒルベルト空間Ｈ内における(Ｖ，＜＞）の自然な内積＜＞から決まる計量ｍ、即ち、有限次元ベクトル空間Ｖからその双対空間Ｖ*への線形写像として、

は一般に疎な行列として定まる(ｍ^T＝ｍ、転置に対して不変)。具体的には、それぞれの格子点Ｎに張り付いた近接胞内では非零であって、その他ではゼロとなる基底ｅ［Ｎ］（，ｘ）を利用し、該基底に対して実数ｆ［Ｎ］を配分することで、領域全体で定義された関数ｆ（ｘ）が、

という有限の線形和で表現される。ここで、これらの基底に配分された実数ｆ［Ｎ］の集合{ｆ［Ｎ］｜Ｎは全ての格子点}を有限次元ベクトル空間Ｖの要素ｆ＝（ｆ［Ｎ］）とみなすと、計量ｍは、

と表すことができる。
なお、計量ｍは、実数Ｒへの双線形写像である内積＜＞として、

が自然なペアリング［，］：Ｖ×Ｖ* → Ｒを利用して、

と一致するように決定される。つまり、実有限次元ヒルベルト空間Ｈは、（Ｖ，ｍ）で特徴付けられているとみなすことができる。
ここで、有限要素法の基底の台に包含関係が存在する場合について考える。一般に基底の台に包含関係が存在する場合、粗い解像度の基底は細かい解像度の基底で必ず表現することができ、さらに、基底が節点とその解像度によって特徴付けされることとなる。そして、このとき有限要素法で利用する有限次元ヒルベルト空間Ｈは、基底の包含関係より定まるフィルターに対して、フィルター付きヒルベルト空間になっている。

ここで、フィルター付きヒルベルト空間により関数空間Ｆを近似するには、加法的多重レベル法を利用するのが一般的である。当該方法によれば、計量ｍ^-1の算出にあたり計算コストを削減することが可能であるからである。

一方、近年、解析対象に適した最適な格子を構成することで、有限要素法による解析精度を向上させる試みが種々なされている。しかし、最適化格子の構成を行う場合、一般に従属的節点の問題を伴う。従属的節点とは、例えば、図９の（ａ）に示す構造格子において９０１〜９０３に示す点をいい、基底関数では表現できない幾何学的構造を有するという特徴がある。このため、フィルター付きヒルベルト空間により関数空間Ｆを近似するにあたり、加法的多重レベル法をそのまま適用できない。

そこで、例えば、特開平８−８７４９０号公報では、かかる問題を解決すべく、従属的節点のある格子上でのウェーブレットの階層性等を利用した共役勾配法およびその前処理方法を提案し、加法的多重レベル法の適用を実現している。
特開平８−８７４９０号公報

しかしながら、ウエーブレット基底には、一般にそれに自然な双対基底が用意されており、その双対性は有限要素法のものとは異なっている。このため、Ｒｉｔｚ（リッツ）の変分原理を適応する際の双対性と矛盾することなく双対基底を与えることは困難であり、実際、上記従来技術（特開平８−８７４９０号公報の発明）においても１次元においては差分法を利用することで回避している。さらに、２次元以上においては、対象とする行列式が異なることは自明である。この結果、上記従来技術を利用した有限要素法解析によれば、解析結果の精度が低下してしまうという問題があった。

さらに、一般に従属的節点を含む格子上のウェーブレット基底の力学的自由度の取り方は明白ではない。そのため、有限要素法の解析アルゴリズムが複雑になってしまうという問題があった。

本発明は、上記課題に鑑みてなされたものであり、従属的節点を有する構造格子に対する有限要素法解析を行うにあたり、基底の台に包含関係が存在する場合において、加法的多重レベル法を利用しても、精度よく、簡便に、かつ高速に解析処理することが可能な情報処理技術を提供することを目的とする。

上記の目的を達成するために本発明に係る制御プログラムは以下のような構成を備える。即ち、
有限要素法によるプログラム処理をコンピュータによって実現させるための制御プログラムであって、
通常節点に張り付いた基底関数を基底とするベクトル空間をＷ、該通常節点に張り付いた基底関数と従属的節点に張り付いた基底関数とを一次独立な基底とするベクトル空間をＶ、従属的節点に張り付いた基底関数を通常節点に張り付いた基底関数によって表現することで得られるＷとベクトル空間として同型な前記ベクトル空間Ｖの部分ベクトル空間をＶ’、前記ベクトル空間Ｗの前記部分ベクトル空間Ｖ’への矩形型の変換行列をＳ、前記ベクトル空間Ｖから前記ベクトル空間Ｖ’への変換行列として与えられる行列Ａとし、前記部分ベクトル空間Ｖ’のベクトルｂが、
Ａｕ＝ｂ
を満たすベクトルｕを共役勾配法により求める場合において、
前記部分ベクトル空間Ｖ’を前記ベクトル空間Ｖの部分空間とみなす工程と、
前記ベクトル空間Ｖ上の共役勾配法において、前記変換行列Ｓの転置行列Ｓ^tに対して、初期値をＳ（Ｗ）からとる工程と、
前記行列ＡをＳ^tＡＳに置き換える工程とを備える。

従属的節点を有する構造格子に対する有限要素法解析を行うにあたり、基底の台に包含関係が存在する場合において、加法的多重レベル法を利用しても、精度よく、簡便に、かつ高速に解析処理することが可能となる。

以下、必要に応じて添付図面を参照しながら本発明の実施形態を詳細に説明する。

１．本発明の概要
はじめに本発明の概要を説明する。本発明にかかる情報処理方法は、例えば、有限要素法解析プログラムの一機能として、当該有限要素法解析プログラムがコンピュータによって実行されることで実現されうる。

また、本実施形態にかかる有限要素法解析プログラムは、従属的節点が含まれる構造格子であっても精度よく、簡便に、かつ高速に有限要素法解析を行うことを目的とするプログラムであり、同プログラムにおいて計量ｍを算出するにあたり、加法的多重レベル法を利用することを前提としている。そして、本発明の特徴は、このような前提のもとで従属的節点が含まれる構造格子が解析対象となった場合であっても、共役勾配法の適用を可能とした点にある。これにより、有限要素法解析において生じる連立一次方程式の解を求める過程で、精度よく、簡便に、かつ高速な処理を実現することができる。

より具体的に述べると、一般に、フィルター付きヒルベルト空間は、基底関数の分解により簡素化が可能であるところ、従属的節点が含まれる構造格子は、力学的自由度を持たないため、自由度の数え上げが有限次元ベクトル空間Ｖでは困難である。

そこで、力学的自由度によって張られる空間Ｗから出発してＷの情報σによりＶ_wを構成し、ιにより有限次元ベクトル空間Ｖに自然に埋め込んで力学系をつくり、Ｖにおける加法的多重レベル加速を利用した後に、Ｗまたは、Ｖ_wで、自由度の勘定を行うこととした。

換言すると、本発明にかかる制御プログラム（有限要素法解析プログラム）は、有限次元部分ベクトル空間Ｖ’を有限次元ベクトル空間Ｖの部分空間とみなす工程と、有限次元ベクトル空間Ｖ上の共役勾配法において、変換行列Ｓの転置行列Ｓ^tに対して、初期値をＳ（Ｗ）からとる工程と、行列ＡをＳ^tＡＳに置き換える工程とを備えることを特徴とする。

２．装置構成
図１に本発明の一実施形態にかかる情報処理装置の構成を示す。本実施形態では、当該情報処理装置において有限要素法解析プログラムが実行される。

図１において、１０１はＣＲＴ等の表示装置、１０２はキーボードやマウス等の各種入力装置を表す。ハードディスク１０３には各種プログラムが記憶されており、かかるプログラムはバス１０７を介してメモリ１０５に読み込まれ、ＣＰＵ１０６において処理される。

ハードディスク１０３には有限要素法プログラムを実行するうえで必要なプログラムが記憶されている。なお、ハードディスク１０３に記憶されたこれらのプログラムは、フロッピディスク１０８に記憶されていてもよく、フロッピディスクドライブ１０４を介して読み込まれてもよい。なお、参考までに図１３に情報処理装置の外観を示す。図中の参照番号は、それぞれ図１の参照番号と対応付けて付してある。

３．処理フロー
図２は本発明の一実施形態にかかる情報処理装置において実行される有限要素法解析プログラムのうち、有限要素法解析過程において生じる連立一次方程式の解を共役勾配法により求める処理の流れを示すフローチャートである。

図２に示す通り、ステップＳ２０１で、共役勾配法サブルーチンに入り、ステップＳ２０２で初期化を行う。ステップＳ２０３で本体部分の計算を行い、ステップＳ２０４で条件判定を行い、条件が満たされた場合はステップＳ２０５に進み、処理を終了する。そうでない場合は、再びステップＳ２０３の操作を行う。

４．有限要素法の説明
以下に、本実施形態にかかる有限要素法解析プログラムにおいて用いられる有限要素法の詳細について説明することとする。

４−１ 有限要素法の概要
有限要素法の基本的な考え方は、境界値問題として定まった偏微分方程式の解を、微分作用素の定義域である関数空間Ｆを実有限次元ヒルベルト空間Ｈで近似することで求めるというものであるＲｉｔｚの方法に立っている（ここで関数空間Ｆが複素数値の場合も実行列表現があるので、一般性を失わず、実の場合のみを述べている）。ヒルベルト空間Ｈの内積としてはルベーグの関数空間の２乗積分（Ｌ２）ノルム（岩波数学辞典第三版、岩波、1994年発行p.61）から誘導される自然な内積を利用する。また、近似の評価を弱い位相で行うというのもＲｉｔｚの方法と同様である。有限次元ヒルベルト空間Ｈは有限次元ベクトル空間Ｖと内積＜＞の組(Ｖ，＜＞）によって定まるが、この場合ＶからＦへの線形写像が単射となるように構成されており、Ｆへの様々な作用（例えば微分作用）はその自然な制限としてＶ上の作用として定義できる。より正確には、ポアソン方程式に対するエネルギー汎関数は、ベクトルの平行移動等を必要とするので、アファイン空間上にヒルベルト空間の概念を拡張しなければ定式化することもできない。しかしながら、アファイン空間の接空間がベクトル空間である事実から、ディリクレ問題を解く場合は特にその差異を意識することなく定式化できる。従って、概念の簡略化のため、ディリクレ問題を取り扱い、ヒルベルト空間で以下議論をする。

多くの有限要素法の特徴としては、有限次元ヒルベルト空間の基底関数として極めて狭い台（非ゼロ集合の閉包）をもつ関数を利用することである。そのため、有限次元ベクトル空間Ｖに対する上記ヒルベルト空間内Ｈ(Ｖ，＜＞）の自然な内積＜＞から決まる計量ｍ、即ち、Ｖからその双対空間Ｖ*への線形写像として、

が一般に疎な行列として定まる(ｍ^T＝ｍ、転置に対して不変)。具体的には、それぞれの格子点Ｎに張り付いた近接胞内で非零でその他ではゼロとなる基底ｅ［Ｎ］（，ｘ）を利用し、該基底に対して実数ｆ［Ｎ］を配分することで、領域全体で定義された関数ｆ（ｘ）が、

という有限の線形和で表現されている。これらの基底に配分された実数の集合{ｆ［Ｎ］｜Ｎは全ての格子点}を本明細書ではベクトル空間Ｖの要素ｆ＝（ｆ［Ｎ］）と見なしている。

ここでｍは、実数Ｒへの双線形写像である内積＜＞として、

が自然なペアリング［，］：Ｖ×Ｖ* → Ｒを利用して、

と一致するように決定される。これにより、ヒルベルト空間Ｈは（Ｖ，ｍ）で特徴付けされているとみなしてもよいので、以下では、（Ｖ，ｍ）による特徴付けを利用する。また、関数空間Ｆへの微分作用素も上記の自然な制限により、有限次元ヒルベルト空間Ｈ上の有界な（つまり、実数係数有限次元行列表示をもつ）作用素と見なすことができる。上記の制限が作用素環としての準同型になるとは限らないので、作用素環としての自然な包含とはなっていない。しかしながら、ベクトル空間Ｖを実数値ベクトル空間としての関数空間Ｆの自然な包含であると考えることは重要であり、有限要素法は有限次元作用素環上に定義された数値計算技法であると考えることも重要である。
また、上記に述べたように極めて狭い台を持つ関数をＶの基底関数とするために、計量ｍや多くの物理的な微分作用素のＨでの行列表現が、通常非常に疎な行列（非ゼロの要素の数の全自由度に対する比が非常に小さい行列）となる。そのため数値的な計算が非常に簡単に行える。

このように、有限要素法という仕組み自体には汎用性があり、かかる汎用性を利用することで、有限要素法のプログラムの自動生成や有限要素法の基本部分のライブラリとしての提供も可能である。

ここで有限要素法ライブラリとは行列演算ライブラリや関数ライブラリと同様に、ユーザであるプログラマーが自作のプログラムに組み込んで個々の問題に応じて変数を定義し、ライブラリの提供する関数等を組み合わせて、問題に応じた有限要素法プログラムを完成させるための基礎的な変数宣言及びサブルーチンや関数の集合をいう。そして有限要素法による処理がライブラリにおいて数学的に定式化されているため、個別の問題（例えば、材料力学、連続体力学、電磁気学等）に応じた有限要素プログラムを容易に作成することができるという利点がある。

４−２ 基底の台に包含関係がある場合の有限要素法
本実施形態にかかる有限要素法解析プログラムは、基底の台に包含関係がある場合を前提としている。そこで、基底の台に包含関係がある場合の有限要素法について以下に述べる。

簡単のため、図３と図４に示すように、まず１次元の場合について述べる。図３、図４において横軸は空間軸であり、縦軸が関数の値を示している。なお、ここでは、空間軸をｘ軸とする。図３（ａ）と（ｂ）は不図示のある周期的関数を、区間的１次関数で近似したものである。図３（ａ）と（ｂ）とはそれぞれ異なる解像度を持っている。ここで、図３（ａ）の頂点は、{０，ａ，２ａ，３ａ，４ａ，５ａ｝として、５ａを原点０と同一視する周期的境界条件を課している。他方、図３（ｂ）の頂点は{０，ａ／２，ａ，３ａ／２，・・・，５ａ}となるように、点を取っている。これらの区分的１次曲線を、図４に示す各頂点に局在する基底関数の和として眺めることが、有限要素法の慣例的な構成方法である。図４（ａ）に示す基底関数はｊ番目の頂点に対して、

と定義されている。但し、（ｊ＋１）及び（ｊ−１）は、法５で計算されている。他方、図３（ｂ）に示す細かい格子に対する基底もｂ＝ａ／２として、

となっている。これらの基底を利用すると、図３（ａ）と（ｂ）は、図３（ｃ）と（ｄ）に示すように、各基底の線形和として表現できる。ここで、重要なことは、それぞれの基底が節点と、その解像度によって特徴付けられている点である。
更に、図４（ｃ）に示すように、ｅ_j ⁽¹⁾は複数のｅ_j' ⁽²⁾で表現できる。この事実の基本的な背景は、ｅ_j ⁽¹⁾の台（非ゼロの領域の閉包）が、ｅ_j' ⁽²⁾の台の幾つかの連結と等しいということによるものである。

このため、図３（ｅ）に示すように、粗い格子上の区分的一次関数は、細かい格子によっても表現できることとなる。つまり、それぞれの区分的一次関数を基底関数の線形和によって形成されるベクトル空間Ｖ⁽¹⁾及びＶ⁽²⁾として見なした際に、台の包含関係、

により、

が得られる。但し、ｓｕｐｐ（ｅ_j' ^(a))は基底ｅ_j' ^(a)の台（非零値の領域の閉包）を示す。
次に２次元の場合を述べる。図５に示すように、２次元の場合も節点によって特徴付けられた基底関数が定義できる。図５（ａ）が図５（ｂ）で示した基底の台の部分である。中央の節点によって、特徴付けられて、基底関数が図４（ｂ）のように定義できる。但し、それ以外の部分に対してはゼロの値を持つことで、基底関数は、領域全域で定義された関数となる。

今、領域として、図５（ｃ）のような格子空間を考える。但し、１次元の場合と同様に、簡単のために周期境界条件を課している。図示された頂点にはそれぞれ、図５（ｂ）に示したのと同じ局在した台を持ち、関数の形も同様な基底関数が定義されている。そのため、領域の任意の関数は、これらの基底関数の線形和で表現される。

この時、図６に示すように、図６（ａ）の格子で定義された基底関数によって生成されたベクトル空間を、図６（ｂ）あるいは図６（ｃ）の格子上で定義された基底関数で生成されたベクトル空間の内包として眺めることは重要なことである。

つまり、有限要素法の基底の台の包含関係より、

に対応して、ベクトル空間Ｖに次のような包含関係

を持つものとできる場合がある。このような減少列あるいは増加列を備えた（ベクトル）空間のことを数学ではフィルター空間と呼ぶ。これらのフィルターの各階層に付けられた番号を数学では次数と呼ぶが、ここでは階層ｊと呼ぶことにする。
図５（ｂ）の関数形は、ｅ_j ⁽¹⁾のテンソル積として形成されていることに注意すると、３次元への拡張は簡単である。更には、図７のような非構造格子においても同様な基底ベクトルの包含関係が構成できるので、その場合においても、同様のフィルター空間が構成できる。

このような幾何学的内包関係とその上で定義される関数の集合に対する理論は「層」の理論として、岡潔らの貢献もあり発達している。より、一般化して、式（１）と式（２）を眺めることにする。数学の慣習に従うと、基底の台の包含関係として、次のような減少列と、

その上での関数のベクトル空間の増大フィルター列

を考える。つまり、集合としての基底の台の包含関係Ｓ_j+1⊂Ｓ_jに対して、自然な単射（１対１写像）ｉ_j：Ｖ_j→Ｖ_j+1が自然に定まることとなる。これをここでは入射と呼ぶ。標語的に書くと「粗い解像度の基底は細かい解像度の基底で必ず表現できる」ことを意味する。層の理論では、ｉ_jを台を制限する意味で制限写像と呼ぶが、多重格子加速法における制限(restriction)と用語が重なるために、ここではより平易な入射という言葉で置き換える。
ここで、ヒルベルト空間Ｈ＝（Ｖ，ｍ）がこの減少列と無矛盾になるものが存在する場合がある。例えば、各部分ベクトル空間に双対ベクトル空間への線形同型写像ｍ_j：Ｖ_j→Ｖ_j*が存在し（ｍ_n＝ｍ）、ｉ_j：Ｖ_j→Ｖ_j+1から自然に定まる線形写像ｉ_j ^#：Ｖ_j+1*→Ｖ_j*（但し、ここで、ｉ_j ^#はｖ∈Ｖ_j+1*とｘ∈Ｖ_jに対して、ｉ_j ^#ｖ（ｘ）：＝ｖ（ｉ_j（ｘ））で定義されるものである。）に対して、ｍ_j＝ｉ_j ^#ｍ_j+1ｉ_jが成り立っている。このようなヒルベルト空間をフィルター付きヒルベルト空間と呼ぶ。

通常、フィルター付きのベクトル空間の増加列が与えられた際に逆の操作、つまり、π_j：Ｖ_j+1→Ｖ_jは自然には定まらないことは重要なことである。この写像の非対称性が層の理論を深いものにしている。しかしながら、上述のフィルター付きヒルベルト空間においては、

が、各部分空間の同型写像ｍ_jを利用して、自然にπ_j：＝ｍ_j-1 ^-1ｉ_j-1 ^#ｍ_jとして定まる。この時、次の性質が証明される。

ここで、＊は、ヒルベルト空間Ｈの内積から定まる各作用素の共役である。この時、次の系として下記のことが分かる。

ここで、それぞれは次のような式変形で証明される。つまり、

より、弱形式の意味でゼロとなり、有限ヒルベルト空間のために強い意味でも等号の成立が言える。また、ｘ∈Ｖ_j+1、ｙ∈Ｖ_jに対して、

より更に明らかである。また、＜ｘ，ｉ_j-1π_jｘ＞≧０は左辺が＜π_jｘ，π_jｘ＞となることより自明である。一般に、

で、階層性のある空間での非可換性を表現している。
このように有限要素法で利用する有限ヒルベルト空間は、上述の基底の包含関係より定まるフィルターに対して、フィルター付きヒルベルト空間になっている。このことを仮定して、フィルター付きヒルベルト空間の性質を以下に説明する。

４−３ フィルター付きヒルベルト空間の性質
基底の台に包含関係があるフィルター付きヒルベルト空間の性質について説明する。

簡単のために、２次元空間内のポアソン方程式のディリクレ問題として、

を例に取る。但し、Ωは与えられた領域、∂Ωは領域Ωの境界、ｑは電荷分布、Δはラプラス作用素である。
我々はこの微分方程式を近似的に解くために、ｆを有限ヒルベルト空間Ｈから取り、Ｈ⊂Ｆという関係を利用して、次の汎関数、

を定義した後に、Ｈの中で最小化することを行う。Ｈ⊂ＦとしてＨへの微分等は、超関数の意味で定義される。
つまり、問題は「ｆ∈Ｖ_n 但し、Ｊ（ｆ）＝ｍｉｎ_h∈Ｖ_n Ｊ（ｈ）を見つけよ。」という問題を解くことに置き換えられる。そこで、フレッシェの意味での微分∇Ｊ（ｆ）∈Ｖ_nを

として定義する。実際には、フレッシェ微分は、ｍ_n∇Ｊ（ｆ）∈Ｖ_n*として得られるので、これを改めて、ｇｒａｄ（Ｊ（ｆ）)として書き直しておく(ｇｒａｄ（Ｊ（ｆ））＝ｍ_n∇Ｊ（ｆ））。この時、定理として次の事実がわかる。ｖをＶ_nから取ってくる。この時、

を考えると、十分小さいρ₁＞０、ρ₂＞０に対して、

を得る。なぜならば、

であるからである。この事実を利用して計算する方法が、加法的多重レベル法である。同様に

を考えると、十分小さいρ₁＞０、ρ₂＞０、ρ₃＞０に対して、

を得る。証明は同様である。このようにｉ_n-1ｉ_n-2ｉ_n-3π_n-2π_n-1π_nやｉ_n-1ｉ_n-2ｉ_n-3ｉ_n-4π_n-3π_n-2π_n-1π_n等を付加項として作用させても、同様の結果を得ることは自明である。よって、考える階層分作用させることも妥当であると考えられる。
一般に、フィルター空間において、Ｖ_j／Ｖ_j-1を考察することが自然であり、次の直和分解が自然である。

つまり、Ｖのどのベクトルも階層固有の解像度での表現ベクトルＶ_j／Ｖ_j-1の和として構成されることを意味している（これを階層的表現と呼ぶ場合がある）。Ｖ_j／Ｖ_j-1の元はＶ_jに

を作用させることによって生成される。このＶ_j／Ｖ_j-1により上記の収束性を議論してもよい。また、πとｉの非可換性が、フィルター空間の特徴であることは重要であり、以下で述べる加速法もこの非可換性に着目して導かれたものであるといってもよい。
４−４ 加法的多重レベル法
上記４−３に示すように、フィルター付きヒルベルト空間の場合、加法的多重レベル法を利用することが可能となる。そこで、以下に加法的多重レベル法について説明する。

本発明の従来の例としてMarco, Koobus, Dervieuxらの論文を紹介する。（Nathalie Marco, Bruno Koobus and Alain Dervieux, "An Additive Mulitilevel Preconditioning Method", INRIA(Institut National de Recherche en Informatique et en Automatique) Rapport de Recherche N.2310 1994年）論文に従って概略を示す。

論文に従って２次元空間内のポアソン方程式のディリクレ問題を例に取る。但し、Ωは与えられた領域、∂Ωは領域Ωの境界、ｑは電荷分布、Δはラプラス作用素である。

簡単のために、しばらく、二つの階層Ｖ_j⊂Ｖ_j+1との対応について述べる。

有限要素法で現れる関数ｆ（ｘ）は、Ｖ_a(ａ＝ｊ，ｊ＋１)それぞれの格子点Ｎに張り付いた近接胞内で非零でその他ではゼロとなる基底ｅ_a［Ｎ］（，ｘ）を利用し、該基底に対して実数ｆ_a［Ｎ］を配分することで表現される。つまり、該ｆ_a（ｘ）は次のような線形和で

表現されている。これらの基底に配分された実数の集合{ｆ_a［Ｎ］|Ｎは全ての格子点｝を本明細書ではベクトル空間Ｖの要素ｆ_a＝（ｆ_a［Ｎ］)と見なしている。それぞれには、ｍ_j、ｍ_j+1という測度が、

として定義されておりｍ_a：Ｖ_a→Ｖ_a*という自然な対応を与えている。自然な入射として、

が定まる。他方、双対空間の性質より、自然な制限ｉ_j ^#:Ｖ_j+1*→Ｖ_j*が

として定まる。よって、この時、ｍ_j＝ｉ_j ^#ｍ_j+1ｉ_jが成り立つこと、即ち[ｍ_jｗ，ｖ]＝[ｍ_j+1ｉ_jｗ，ｉ_jｖ]が成り立つことは、定義より自明である。つまり、有限要素法の基底の台から決まるフィルター付きベクトル空間上の有限ヒルベルト空間は、フィルター付きヒルベルト空間となる。

を具体的に書き下すとｅ_j［Ｎ］での係数が

となる。但し、ｓｕｐｐ（ｅ_a［Ｎ］)は基底ｅ_a［Ｎ］の台（非零値の領域の閉包）を示す。ここで、

となることを利用した。図３（ｃ）を細かい階層の基底で展開し積分した後に足し合わせても、粗い階層で積分しても積分値は同じなので、ｉ_jｅ_j[Ｎ']をｅ_j[Ｎ']（ｘ）と記している。
しかしながら、ｍ_j ^-1の計算は一般に計算コストが掛かる。正確には、ｍ_n ^-1の計算の代わりに、π_nの類似物で置き換えて計算に利用しようとするものが、加法的多重レベル法である。つまり、π_j+1ｆ_j+1を、次の操作で置き換える。

ここで、Ｖｏｌ（Ｕ）は、領域Ｕの面積を示す。
式（３）で述べたように、π_j+1ｉ_jは恒等写像であった。この恒等性は、体積保存程度のものでなりたっているというアイデアを背景にして、上記のΠ_jは定義されている。従って、

を考えると、十分小さいρ₁＞０と、ρ₂＞０と、π_nをΠ_nとする近似のもとで、

が成り立つことが期待される(上記従来例の文献において、このことを「証明」しているが、特殊な測度においては証明されるが、通常の有限要素法で利用されている測度ｍ_nに対しては一般には成り立たない。しかしながら、π_jをΠ_jで近似する意味において正当化されるものであり、その近似の精度、所謂ランプ近似と呼ばれるものと同等であり、求める問題や、格子の出来具合によって、その精度は変動すると考えられる。従来例や、本実施例で扱う物理的な作用素に対しては子収束を早めることから、上記の近似はよい近似であると思われる)。
更に、Π_nをπ_nの代わりとして利用し、分子の操作をＶ_n*内の通常の有限和の操作として理解することで、Π_n:Ｖ_n*→Ｖ_n-1と理解する。つまり、Π_nのみｍ_n-1 ^-1ｉ_n-1 ^#の近似であると見なす。この考え方により、Π_nを、共役勾配法においてｍ_n-1：Ｖ*→Ｖを計算しなければならない工程において、ｍ_n-1に代わる別の行列を施すことで、最小化を加速するという前処理という部分に関わらせるものとなる。しかしながら、π_jをΠ_jに置き換えたための補いのために、ρを最適化しようとするのが従来例の手法であった。本発明においても、この考え方を利用する。

実際には、有限要素法において、フレッシェ微分は、∇Ｊ（ｖ）よりもｖ∈Ｖ_n＝Ｖに対してｇｒａｄ（Ｊ（ｖ））∈Ｖ_n*＝Ｖ*として得られる。また同様に、ｆ∈Ｖ_n＝Ｖに対して、Ｌａｐ（ｆ）をＬａｐ（ｆ）:＝ｍ_nΔ（ｆ）∈Ｖ_n*＝Ｖ*として定義しておく。更にＰｒｅＣＣとして、

と定義し、これらは、ＰｒｅＣＣ:Ｖ_n*→Ｖ_nと見なすこととする。
このとき、十分小さいρ₁＞０、ρ₂＞０、ρ₃＞０、ρ₄＞０、ρ₅＞０に対して、ｖ'＝ＰｒｅＣＣ（ｖ）は

を得ることが期待される。ある初期値ｘをＶから取って来て、

となる。
４−５ 従属的節点
次に、本発明の適用対象となる、従属的節点について説明する。近年、問題に適した最適化格子を構成しようとする流れが生じている。このような考え方は従来にも存在したが、近年、図８に示すような、構造格子を利用したものが現れてきた。この場合、図９（ａ）の黒点で示すように、従属的節点（slaved node あるいはhanged node）の問題が生じる。それは、図４（ｂ）の基底関数では表現できない幾何学構造をもっている点である。

ここで、従属的節点の線形代数的な観点からの性質について述べる。まず、従属的節点を表現するために、物理的な自由度を表現するベクトル空間Ｗ

に対して、Ｖ⊃Ｗとなる従属的節点を一次独立基底として含む空間Ｖを

として用意する。Ｖをｅ_lで張られる空間ＷとＳ_aで張られる空間Ｖ_sとで直和分解する。つまり、

とする。この時、自然な射影として、

が存在する。
他方、節点上の基底には、Ｓ_a＝Σ_lＳ_alｅ_lとなる拘束条件があるとする。つまり、(Ｓ_a＝Σ_lＳ_alｅ_l)の関係式をゼロと同一視する意味で、即ち、

と定義すると、σ：Ｗ→Ｖ_wは同型対応がある。この意味で従属的節点を含む空間は、商空間として眺めることができる。通常、商空間の位相は元のものと異なるが、ここでは、Ｖ_w⊂Ｖとして、Ｖからの自然な位相を導入する。このことにより、次の包含写像、

が定義でき、同時に、その双対射

が定義される。
また、上記のことより、σ・ｐ:Ｖ→Ｖ_wが定義できる。

この時、双対空間の性質より、

が定義できることとなる。行列表現した際には、(σ・ｐ）^#は(σ・ｐ）の転置行列となる。よって、

となり、(σ・ｐ）^#:は実質的には、Ｖ*内の部分空間であるＷ*への射影する射影作用素ｐ*が作用していると考えてもよい。
従属的節点は力学的自由度を持たないので、系の関数ｆ（ｘ）はＷでのみ構成されている。しかしながら、Ｖ自身の方が幾何学との対応が明白である場合が多い。つまり、Ｖの方が、基底関数の台によるフィルター関係が簡単である。

つまり、従来の例で述べたように、Ｖの方には有限要素法の基底の減少列が与えられているとする。

その上での関数のベクトル空間の増大フィルター列

が存在したとする。
このことをより明確に述べるために、ポアンカレ双対性の意味での、節点の双対である胞を基本としたフィルター構造を考えると自然である。図４（ａ）の節点に張り付いた基底に対して、図１０（ｃ）、（ｄ）に示すように、胞Ｃを基底とする２つの基底関数に分解して考えると、台のフィルター構造がより簡単となる。２次元の場合は図５（ｂ）の基底を図１１に示したような胞Ｃに対する４つの基底関数に分解して考えればよい。このような基底の取り替えを行うと、図８のような階層構造と従属的節点を持つ場合も、自由度の勘定さえ正確であれば、フィルター構造は簡単になる。実際、胞による表現を利用すると、式（５）の計算は非常に簡単になる。

しかしながら、従属的節点は力学的自由度を持たないので、系の自由度の数え上げは、Ｖでは一般に困難である。

そこで、構成の簡単なＷから出発して、Ｗの情報をσにより、Ｖ_wで構成し、ιによりＶに自然に埋め込んで力学系をつくり、Ｖにおける加法的多重レベル加速を利用した後に、Ｗまたは、Ｖ_wで、自由度の勘定を行うこととした。

５．具体的な実施例
５−１ 実施例１
Ｗのみ一次独立な基底が存在するので、有限要素法をＷ上で形成し、Ｗ上の線形一次方程式を解くことにより問題を（近似的に）解くということは実に明快である。

そこで、Ｗから決まるヒルベルト空間をＶ上のヒルベルト空間Ｈと無矛盾に、Ｈ_wと定義することとする。つまり、Ｈ_wの測度行列（ｍ_w:Ｗ→Ｗ*)をｍ_w:＝σ^#ι^#mισによって定義する。

つまり、行列としてはσ^# ι^#＝（ισ）^Tとして実現されている。

自然なペアリングが、

として定義されている。これとｍ_wを通して、Ｈ_wの内積＜,＞ｗが定義されている。
これらを利用して、従属的節点の存在する場合の、前処理付共役勾配法のアルゴリズムを示す。

簡単のために、２次元空間内のポアソン方程式のディリクレ問題として、

を例に取る。但し、Ωは与えられた領域、∂Ωは領域Ωの境界、ｑは電荷分布、Δはラプラス作用素である。
我々はこの微分方程式を近似的に解くために、ｆを有限ヒルベルト空間Ｈ_w＝(Ｗ，ｍ_w)から取り、Ｗ上で定義される汎関数として、

を定義した後に、Ｈ_wの中で最小化することを行う。Ｖ⊂ＦとしてＶへの微分等は、超関数の意味で定義される。ここで、ｑ'はｑを近似したＷの元である。この時、問題は「ｆ∈Ｗ 但し、Ｊ_w（ｆ）＝ｍｉｎ_h∈W Ｊ_w（ｈ）を見つけよ。」という問題を解くことに置き換えられる。
しかしながら、そこで、Ｖ上には、フィルター構造等が自然に定義されているので、従来例と同様に、次のような定式化に変更することは有効である。つまり、ισ:Ｗ→Ｖ⊂Ｆという関係を利用する。まず、Ｖの元ｆに対して汎関数が定義されているとする。

但し、ｑは式（６）のｑを模擬したもので、ι(Ｖ Ｗ)内の要素である。この時、問題を「ｖ∈Ｗ 但し、Ｊ(ισ（ｖ）)＝ｍｉｎ_h∈W Ｊ(ισ（ｈ）)を見つけよ。」という問題を解くことに置き換える。
フレッシェの意味での微分∇Ｊ（ｆ）∈Ｖ_nを

をとして定義する。実際には、フレッシェ微分は、ｍ_n∇Ｊ（ｆ）∈Ｖ_n*として得られるので、これを改めて、ｇｒａｄ（Ｊ（ｆ））として書き直し、ｇｒａｄ（Ｊ（ｆ））＝ｍ_n∇Ｊ（ｆ）∈Ｖ_n*＝Ｖ*を定義する。また同様に、Ｌａｐ（ｆ）をＬａｐ（ｆ）:＝ｍ_nΔ（ｆ）として定義しておく。更にＰｒｅＣとして、

と定義しておく。但し、ρ_jは正の実数である。
このとき、π_nｍ_n ^-1：Ｖ*はＶに写像するものであるので、

とみなす。
ある初期値ｘをＷから取って来て、

となる。
図１２に従うと、ステップＳ１２０１で、共役勾配法サブルーチンに入り、ステップＳ１２０２で初期化を行う。ステップＳ１２０３で本体部分の計算を行い、ステップＳ１２０４で条件判定を行い、条件が満たされた場合はステップＳ１２０５に向かい終了する。そうでない場合は、再びステップＳ１２０３の操作を行う。

ここで、ｒ（ｋ）はＷ*の元である。ｇ（ｋ）、ｈ（ｋ）、ｘ（ｋ）はＷの元であるので、自由度は常に保たれている。ι及びι^#は型変換のようなもので、実際はなにもしない。ｐも射影であるので、ｐ^#もｐも部分的にゼロを代入する作用素である。σは行列であり、従属的節点上の関数の値を形成する。σ^#はσの行列表示の転置行列として構成される。計算を行うと、前者の方が、後者より収束が早くなった。

本実施例においてはρ_jの値としてすべて１の場合と、ρ₀のみ１で他をゼロとするものを利用した。後者は、フィルター構造に沿う加法的多重レベル加速を行わない場合である。もちろん、さまざまなρ_jに対して計算を行ってもよい。

ここで、特開平８−８７４９０号公報との違いについて述べておく。特開平８−８７４９０号公報において、σの役割をするのは、不完全離散ウエーブレット変換行列τ（特開平８−８７４９０号公報においてはＷと記している。）であるが、τが正方行列であるのに対して、σは矩形行列である。なぜならば、Ｗの次元とＶの次元は従属的節点の存在のために異なるためである。しかしながら、本発明において、力学的自由度の勘定はＷで行えばよく、実に明快であるのに対して、特開平８−８７４９０号公報においては力学的自由度を勘定することが難しいことが分かる。また、有限要素法のような弱形式によって、弱解を求めるというバックグランドがないために、差分作用素の定式化に自由度があり、実際に問題を解く際に困難が生じるのに対して、本発明は、リッツの原理から出発しており、力学的に安定な解法であり、かつ、アルゴリズムが明確である。そのために、ラプラス方程式以外の方程式に対しても、困難なく定式を行うことができる。また、弱形式という背景を持たないために、特開平８−８７４９０号公報の前処理行列は、本実施例のものとは異なり、式（４）の近似としての意味もない。

本発明提案のアルゴリズムに従うと、従属的節点をもった系においても、共役勾配法が合理性を伴って構成できる。また、パラメータρ_jの値により、加法的多重レベル加速を取り入れた前処理付き共役勾配法から、通常の共役勾配法まで、従属的節点をもった系においてもスケーラブルに取り扱うことが可能であり、有限要素法ソルバーの収束パラメータとして、簡便に利用することができる。

また、加法的多重格子加速を実行的に行った場合は、従来のものより収束が早くなる。

５−２ 実施例２
実施例１と同様の作用素と問題に対して以下のようなアルゴリズムが提案できる。ある初期値ある初期値ｘ₀をＷから取って来て、ｘ＝σ（ｘ₀)∈Ｖ_wと定義する。

となる。
図１２に従うと、ステップＳ１２０１で、共役勾配法サブルーチンに入り、ステップＳ１２０２で初期化を行う。ステップＳ１２０３で本体部分の計算を行い、ステップＳ１２０４で条件判定を行い、条件が満たされた場合はステップＳ１２０５に向かい終了する。そうでない場合は、再びＳ１２０３の操作を行う。

ここで、ｒ（ｋ）はＶ_w*の元であり、ｇ（ｋ）、ｈ（ｋ）、ｘ（ｋ）はＶ_wの元であるので、ペアリングはＶ上のものを利用している。

本実施例においてもρ_jの値としてすべて１の場合と、ρ₀のみ１で他をゼロとするものを利用した。後者は、フィルター構造に沿う加法的多重レベル加速を行わない場合である。この場合は、(σ・ｐ)^#を、前処理行列と思う共役勾配法になっている。

計算を行うと、前者の方が、後者より収束が早くなった。
もちろん、従来例の場合と同様にさまざまなρ_jに対して計算を行ってもよい。

本実施例のアルゴリズムは、実施例１と比較して、σの作用が少なくて済んでいることは重要なことである。よって演算回数が実施例１より少ない。特開平８−８７４９０号公報との差異については実施例１と同様である。

５−３ 実施例３
実施例１と２の計算は、図１４のような３次元格子にも簡単に拡張できる。更には、図６に示したような従属的節点を持つ非構造格子においても、同様に適応できることは自明である。なお、本発明は、複数の機器（例えばホストコンピュータ、インタフェイス機器、リーダ、プリンタなど）から構成されるシステムに適用しても、一つの機器からなる装置（例えば、複写機、ファクシミリ装置など）に適用してもよい。

また、本発明の目的は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記憶媒体を、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ（またはＣＰＵやＭＰＵ）が記憶媒体に格納されたプログラムコードを読出し実行することによっても、達成されることは言うまでもない。

この場合、記憶媒体から読出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記憶した記憶媒体は本発明を構成することになる。

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

また、コンピュータが読出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているＯＳ（オペレーティングシステム）などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。

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

従属的節点を有する構造格子にも適用可能な有限要素法解析プログラムとして利用することができる。

本発明の一実施形態にかかる情報処理装置の構成を示す図である。 本発明の一実施形態にかかる情報処理装置における有限要素法解析プログラムのうち、連立一次方程式の解を求める処理の流れを示すフローチャートである。 横軸に空間軸、縦軸に関数の値を記載した様子を示す図である。 横軸に空間軸、縦軸に関数の値を記載した様子を示す図である。 ２次元空間において、節点によって特徴付けられた基底関数を示す図である。 ２次元空間における基底の台を示す図である。 非構造格子を示す図である。 構造格子の一例を示す図である。 従属的節点を示す図である。 図４（ａ）の基底に対して、胞Ｃを基底とする２つの基底関数に分解した様子を示す図である。 図４（ｂ）の基底に対して、胞Ｃに対する４つの基底関数に分解した様子を示す図である。 本発明の一実施形態にかかる情報処理装置における有限要素法解析プログラムのうち、連立一次方程式の解を求める処理の流れを示すフローチャートである。 本発明の一実施形態にかかる情報処理装置の外観を示す図である。 ３次元格子の一例を示す図である。

Claims (2)

1. 有限要素法によるプログラム処理をコンピュータによって実現させるための制御プログラムであって、
   通常節点に張り付いた基底関数を基底とするベクトル空間をＷ、該通常節点に張り付いた基底関数と従属的節点に張り付いた基底関数とを一次独立な基底とするベクトル空間をＶ、従属的節点に張り付いた基底関数を通常節点に張り付いた基底関数によって表現することで得られるＷとベクトル空間として同型な前記ベクトル空間Ｖの部分ベクトル空間をＶ’、前記ベクトル空間Ｗの前記部分ベクトル空間Ｖ’への矩形型の変換行列をＳ、前記ベクトル空間Ｖから前記ベクトル空間Ｖ’への変換行列として与えられる行列Ａとし、前記部分ベクトル空間Ｖ’のベクトルｂが、
   Ａｕ＝ｂ
   を満たすベクトルｕを共役勾配法により求める場合において、
   前記部分ベクトル空間Ｖ’を前記ベクトル空間Ｖの部分空間とみなす工程と、
   前記ベクトル空間Ｖ上の共役勾配法において、前記変換行列Ｓの転置行列Ｓ^tに対して、初期値をＳ（Ｗ）からとる工程と、
   前記行列ＡをＳ^tＡＳに置き換える工程と
   を備えることを特徴とする制御プログラム。
2. 請求項１に記載の制御プログラムを格納した記憶媒体。

Images