第5話 あっと驚くガウスの数値積分法
「応力は積分点で出力されています」と言うと、有限要素法の経験の浅い方にはきょとんとした顔をされることがよくある。ここでいう積分点とはガウスの数値積分点のことである。
有限要素法での重要な局面の1つに各要素の剛性マトリックスを求めることがある。その剛性マトリックスは要素内で仮定した変位分布から誘導された関数を積分することで求まる(ここでは、一番ポピュラーな変位型有限要素法を想定している)。2次元要素や板要素では2変数の面積分となり、ソリッド要素では3変数の体積分となる。
いま、1要素内での応力や歪みが一定とみなしてもいいほどに、メッシュ分割が充分に細かくされていると仮定しよう。すると、要素剛性マトリックスを形成する積分式の被積分関数は一定値となり、積分の外に吐き出され、最終的に積分項は要素の面積(体積)を求める式が残るだけとなる。結局、有限要素法では要素の面積、体積を精度よく求めることが重要課題となる。
面積を求めると言えば、筆者の高校時代に、定積分の応用においてシンプソンの積分公式を使って船体の縦断面の面積を求めたことを今でも覚えている。この公式はニュートン・コーツ系の積分公式の一種で、積分範囲を等間隔に分割しなければならない宿命にある。それで、現在の有限要素法の主流となっている、要素形状の歪みが、ある程度許容されるアイソパラメトリック系要素の面積分、体積分を求めるには非常に不利となってしまう。高次要素では曲線、曲面を対象とする積分が必要となる。ここで登場するのがガウスの数値積分公式である。
筆者が初めてこの積分公式を利用したときは、その抜群の威力に感嘆したものである。シンプソンの公式に比べて、ほんの数点の関数値だけで非常に精度のよい積分値が求まるのである。この旨みを知ってからというも、筆者は何かと愛用し続けている。つい最近も下式のような積分値を求める必要があり、使用したことがある。
この積分を解析的に求めようとするなら、恐らく漸化式を利用していけば解けてしまうのかもしれないが(ひょっとしたら、無理なのかも)、そんなことを考えているより、ガウスの積分公式を適用すれば、10数ステップぐらいのプログラミングで済んでしまうのである。
ところで、筆者がガウスの積分公式を初体験した時、どんな被積分関数であっても数少ない一定の積分点位置だけで、どうしてあんな精度のよい積分値が求まるのか不思議でならなかった。まるで、魔法の玉手箱を使っているような気持ちがしたものである。この公式の誘導には、さぞ、高等数学が駆使されているのだろうと想像していたが、実は、意外にも簡単なことを知り2度目の驚きであった。準備知識として、ルジャンドルの多項式とラグランジュの補間式の2つがあればよいのである。
この公式誘導のエッセンスを言えば、ある定義範囲での任意の関数は直交関数の展開式で表現されるということから始まる。直交関数というのは定義域で下記の性質を持つような多項式のことであり、まず、真っ先に思いつくのはフーリエ級数での正弦関数、余弦関数ではなかろうか。
他にもチェビッシェフの多項式、エルミートの多項式とたくさんある。そして、ここでの主役であるルジャンドル多項式である。
ルジャンドル多項式は、球の静電ポテンシャルを表す変数係数の2階微分方程式の解式の中から見つけ出されている。この多項式の初めの4つを示すと下のとおりである。
さて、ガウスの積分公式(ここでは1変数を扱う)は下の近似式から始まる。
ここで、任意の関数がルジャンドル多項式の重ね合わせで表現できること、また、この多項式が直交性の性質を持つこと(具体的に上の4項を直交性の式に代入するとよく理解できる。ルジャンドル多項式の場合、定義範囲は[-1、1]で重み関数は1である。)を利用すると、x0、x1…が実はルジャンドル多項式をゼロとした方程式の根であることが分かるのである。
いま、2番目、3番目の多項式の場合の根を求めてみると、
ここまでくると、有限要素法の参考書やプログラムリストでよく見かける、何かへんてこな数字列の意味が理解できるというものである。もちろん、冒頭で言った応力の出力位置とはこの位置のことである。
ガウスの積分公式におけるもう1つの定数、重み係数については話のテーマから外れることになるため、積分点を通る関数をラグランジュ補間式で表現することから求められるということだけに留める。
最後に1~3変数の場合のガウスの積分公式を書き下しておく。
2000年11月 記
読者からの寄せられたコメント
馴染みやすい解説ありがとうございます。
連載を終了してから長らく本サイトを覗いていなかったもので、今回久しぶりに覗いてみて初めていただいたコメントに気づきました。返信遅れて申し訳ありませんでした。拙著エッセイをお読みいただきありがとうございます。