(57)【要約】【課題】数値解析に対する感度解析システムを提供する。【解決手段】数値解析を行なう場合、まず、対象となる現象を数学的なモデル（一般には微分方程式）に変換する（Ｓ１１０）。これをコンピュータにより解くために、モデル化に求めた微分方程式から、離散化手法を用いて離散化した代数方程式（離散化式）を求める（Ｓ１２０）。次に、テイラー解析などの解析的方法により代数方程式から再び連続的な微分方程式（修正方程式）を求める（Ｓ１３０）。その修正方程式に対して、時間刻み幅や空間格子幅などの数値計算上のパラメータに関する感度を求めることにより、数値計算結果の数値パラメータに関する感度を評価する（Ｓ１４０）。感度から数値パラメータにより厳密解を推定する（Ｓ１５０）。厳密解が推定できると、計算解と比較することにより、誤差の評価や、計算解の信頼性が評価できるようになる（Ｓ１６０）。【選択図】  図１
2004246744.tif 000002
【特許請求の範囲】
【請求項１】
対象の数学的モデルである方程式から得られる数値解析結果の感度解析システムであって、
前記方程式を離散化して、代数方程式を求める離散化手段と、
求めた離散化代数方程式から、修正方程式を求める修正手段と、
求めた修正方程式の感度解析を行なう感度解析手段と
を備えることを特徴とする感度解析システム。
【請求項２】
請求項１記載の感度解析システムにおいて、
さらに、前記修正方程式から、厳密解を推定する厳密解推定手段を備えることを特徴とする感度解析システム。
【請求項３】
請求項１又は２に記載の感度解析システムを、コンピュータ・システムに構築するプログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、コンピュータによる数値解析結果の感度解析に関する。
【０００２】
【技術的背景】
近年、計算機性能の向上や低下価格化により、手軽に計算機による数値計算により、自然界や産業・工業分野における現象の予測・評価をすることが可能となっている。その場合、まず、対象となる現象を数学的なモデル（一般には微分方程式＝基礎方程式）に変換する。その際、初期値や基礎方程式が含む物理パラメータなど（たとえば、拡散係数や熱伝導率など）が不確かさを有するため、数学モデルは確率論的性格を有している。従来の感度解析法はこうした数学モデルに含まれる物理パラメータの感度を求める手法であった。
しかし、対象となる現象を記述する、数学モデルである微分方程式は、理論的に解を求めることができないため、多くの場合は数値計算により近似的に解を求めることになる。その際、連続的な微分方程式を、計算機上で処理できる代数方程式に変換する離散化近似が施される。しかし、その離散化に伴い、離散化誤差（数値誤差）は避けられない。この数値誤差を評価するため、古くからテイラー展開に基づく方法や周波数解析による方法など多くの手法が提案されてきたが、それらはいずれも定性的な評価しかできず、数値解析結果が含む誤差を定量的に評価することはできなかった。
【０００３】
【発明が解決しようとする課題】
本発明は、数値解析結果が含む誤差を定量的に評価するために、数値解析結果の感度解析システムを提供することである。
【課題を解決するための手段】
上記目的を達成するために、本発明は、対象の数学的モデルである方程式から得られる数値解析結果の感度解析システムであって、前記方程式を離散化して、代数方程式を求める離散化手段と、求めた離散化代数方程式から、修正方程式を求める修正手段と、求めた修正方程式の感度解析を行なう感度解析手段とを備えることを特徴とする。
さらに、前記修正方程式から、厳密解を推定する厳密解推定手段を備えることもできる。
これらの感度解析システムをコンピュータ・システムに構築するプログラムも本発明である。
【０００４】
【発明の実施の形態】
以下、本発明の実施の形態を、添付図面を参照して説明する。
図１に本実施形態のアルゴリズムの概略を示すフローチャートを示す。
図１に示すように、数値解析（数値計算，コンピュータ・シミュレーション等）を行なう場合、まず、対象となる現象を数学的なモデル（一般には微分方程式＝基礎方程式）に変換する（Ｓ１１０）。その際、初期値や基礎方程式が含む物理パラメータなど（たとえば、拡散係数や熱伝導率など）が不確かさを有するため、数学モデルは確率論的性格を有している。
これをコンピュータにより解くために、モデル化に求めた微分方程式（基礎方程式）から、差分法、有限要素法、境界要素法などの離散化手法を用いて離散化した代数方程式（離散化式）を求める（Ｓ１２０）。次に、テイラー解析やフーリエ解析などの解析的方法により代数方程式から再び連続的な微分方程式（修正方程式）を求める（Ｓ１３０）。その修正方程式に対して、時間刻み幅や空間格子幅などの数値計算上のパラメータに関する感度（パラメータが変化した場合に結果がどの程度変化するか）を求めることにより、数値計算結果の数値パラメータに関する感度を評価する（Ｓ１４０）。
【０００５】
修正方程式は、微分方程式を離散化したものである代数方程式に、テイラー展開式を代入することで求められる。
【数１】
2004246744.tif 000003
ここで、打切り誤差に関して次の式が得られる。
【数２】
2004246744.tif 000004
但し、代数方程式は近似解においてテイラー展開されている。
式（２）は、Ｅ（Ｔ）が空間導関数のみを含むように取り扱われている。ここで、Ｌ（Ｔ）＋Ｅ（Ｔ）＝０が修正方程式である。重要なことは、「代数方程式を解くこと」は、「対象とする微分方程式を解くこと」ではなく、この「修正方程式を解くこと」に相当している点である。打切り誤差の項がメッシュ幅（Δｘ）、時間刻み幅（Δｔ）を含むため、修正方程式はメッシュ幅（Δｘ），時間刻み幅（Δｔ）を含む微分方程式となる。
【０００６】
メッシュ幅（Δｘ），時間刻み幅（Δｔ）は入力データであり、パラメータと見ることができる。感度解析は、パラメータの変化が結果に与える影響を評価するものである。したがって、この修正方程式に対してメッシュ幅（Δｘ）、時間刻み幅（Δｔ）の感度解析を行うことで、メッシュ幅（Δｘ），時間刻み幅（Δｔ）を変化させたときの結果の変化を知ることができる。これにより、ある条件における数値解に、任意のメッシュ幅（Δｘ），時間刻み幅（Δｔ）に変化させたときの結果の変化を加えることで、変化させた後の条件における数値解を予測することができる。
変化させた後のメッシュ幅（Δｘ），時間刻み幅（Δｔ）をともに零としたものが推定された厳密解である（Ｓ１５０）。厳密解が推定できると、計算解と比較することにより、誤差の評価や、計算解の信頼性が評価できるようになる（Ｓ１６０）。
【０００７】
＜感度解析＞
感度解析法は、解析対象とする方程式系が含む種々のパラメータの、指定された応答に対する線形感度を求める手法である。本実施形態ではそのうち、随伴系を用いることにより、系が含む全てのパラメータの感度を一度の計算で求めることを可能にする、随伴感度解析法（一般化感度解析法）を用いる。
随伴系は、対象とする非線型方程式系を、それが含むパラメータに関して微分した系が線形系となることに基づいて構成される。
解析対象とする非線型方程式系の解は通常の差分法などにより得られ（これを前進解という）、随伴系は時間を逆行することにより前進解と同様の計算法で解が得られる（これを後進解あるいは随伴解という）。これら２つの結果を任意の時空間上で用いることにより応答汎関数を計算することができる。
このように、随伴系をとることは、適切な境界条件が与えられた時空間の全域を考慮することにあたる。
以下では、白川，  大田，  梶原，  一般化感度解析手法の過渡時熱水力問題への適用，  ＪＡＥＲＩ−Ｍ  ８５−０１７，  （１９８５）や、堀江，  山本，  白川，  非線形系の隋伴系を用いた一般化感度解析法の応用，  茨城講演会講演論文集，  ＪＳＭＥ  Ｎｏ．０２０−３，  ｐｐ．８７−８８，  （２００２）に基づき、その基本的な定式化方法を説明する。
【０００８】
まず、対象とする方程式系を次のように表す。
【数３】
2004246744.tif 000005
上式に次の境界条件を与える。
【数４】
2004246744.tif 000006
ここで、
ρ＝（ρ１，ρ２，…，ρＪ）  ：位相空間における位置ベクトル
α＝（α１，α２，…，αＩ）  ：系が含むパラメータ
Ｘ（ρ）＝（Ｘ１（ρ），Ｘ２（ρ），…，ＸＫ（ρ））：状態ベクトル
Ｎ＝（Ｎ１，Ｎ２，…，ＮＫ）  ：対象方程式系を定める非線形演算子
ρＳ                  ：位相空間上の境界
【０００９】
さて、応答汎関数を任意の関数Ｆで構成する。
【数５】
2004246744.tif 000007
Ｆは例えば、時空間における位置を表すためのデルタ関数δと対象系が計算する物理量の積とする。たとえば、ある量Ｃのモニタリング・ポイントρｔにおける応答を調べたければ、
【数６】
2004246744.tif 000008
とする。次に、Γ＝（Ｘ，α）およびΓ＋δΓ＝（Ｘ＋δＸ，α＋δα）として、応答Ｒの変分を求める。
【数７】
2004246744.tif 000009
ここで、Ｆ’ＸはＦのＸに関するフレッシェ（Ｆｒｅｃｈｅｔ）微分である。
【００１０】
対象方程式系の微分系を導入するために、（３）式，（４）式のフレッシェ微分をとる。
【数８】
2004246744.tif 000010
Ｌ≡Ｎ’Ｘ，Ｓ≡Ｑ’α−Ｎ’αとおくと、
【数９】
2004246744.tif 000011
（８）式（９）式は、δＸについて線形系である。
この段階でδＸを求める手法を、Ｆｏｒｗａｒｄ  Ｓｅｎｓｉｔｉｖｉｔｙ  Ａｎａｌｙｓｉｓ（前進感度解析）という。しかし、その場合、δＸ＝Ｘ’αδαであるから、一つのベクトルαに対して一回の計算が必要になる。
【００１１】
随伴演算子を導入するために、随伴演算子Ｌ＊を次式で定義する。
【数１０】
2004246744.tif 000012
ここで、φは  δＸの随伴関数（ａｄｊｏｉｎｔ  ｆｕｎｃｔｉｏｎ）、Ｐは境界積分項である。また、内積は次式で定義されている。
【数１１】
2004246744.tif 000013
随伴系のソース項と境界条件の設定を行なう。まず、随伴系を次のように表す。
【数１２】
2004246744.tif 000014
ソース項Ｂ＊はδα，δＸに依存しないように決め、境界積分項に含まれるδＸを消去する。境界条件からφのみの関数となった境界積分項を＾Ｐ［φ］で表すと、（１０）式は次のように書ける。
【数１３】
2004246744.tif 000015
（６）式と比較して、
【数１４】
2004246744.tif 000016
とすると、δＲＸ＝〈Ｆ’Ｘ，δＸ〉となるので、δＲα＝〈Ｆ’α，δα〉とおくと、（６）式は次式となる。
【数１５】
2004246744.tif 000017
上式の（１４）式からわかるように、随伴解φが求められれば任意のδα  について、δＲを求めることができる。なお、随伴系の境界条件の設定から時空反転系が自然に導入される。すなわち、時空反転を行うと随伴系の解は元の対象系と同様のコードで計算することができる。
【００１２】
【実施例】
一次元移流方程式を用いて、図１に示した本実施形態の処理を具体的に説明する。移流方程式は、線形の波動方程式である。１次元移流方程式を次に示す。
【数１６】
2004246744.tif 000018
ここで、ｕは背景速度場、Ｃは濃度などのスカラー量である。簡単のため、背景速度場ｕは一定とする。（１５）式は、濃度などのスカラー量Ｃが速度ｕでｘ方向に伝わる現象を表している。
この１次元移流方程式を、時間に関し一次精度オイラー法で、空間に関し一次精度風上差分法を用いて離散化した式は、
【数１７】
2004246744.tif 000019
である。ただし、Ｃｎｉ≡Ｃ（ｎΔｔ，ｉΔｘ）である。
【００１３】
以下では、テイラー展開を用い、（１６）式から修正方程式を導出する。
（１６）式のＣｎ＋１ｉ，Ｃｎｉ−１にテイラー展開式
【数１８】
2004246744.tif 000020
を代入すると、
【数１９】
2004246744.tif 000021
を得る。ここで、左辺が移流方程式と一致することから、右辺は離散化より生じた打切り誤差である。打切り誤差の項の意味は、時間の導関数を空間の導関数に置き換えることで、より簡単に解釈することができる。Ｃｔｔを空間の導関数に置き換えるために、（１７）式を時間ｔについて偏微分すると、
【数２０】
2004246744.tif 000022
を得る。また、（１７）式を空間ｘについて偏微分し、−ｕを掛けると
【数２１】
2004246744.tif 000023
となる。（１８）式と（１９）式を足すと、以下の式を得る。
【数２２】
2004246744.tif 000024
同様にして、Ｃｔｔｔ，Ｃｔｔｘ，Ｃｘｘｔについて次の式が得られる。
【数２３】
2004246744.tif 000025
（１８）式，（２０）式，（２１）式を組み合わせると、以下の式が求められる。
【数２４】
2004246744.tif 000026
【００１４】
（２２）式のような式は、修正方程式と呼ばれている。この方程式は、偏微分方程式に差分法を適用したときに実際に解いている偏微分方程式となる。修正方程式である（２２）式の右辺は、上で述べたように、それが元の偏微分方程式とその差分近似との差を表していることから、打切り誤差である。よって、修正方程式の右辺の最も低い次数の項がその差分法の次数を与えることになる。今回のケースでは、最も低い次数の項が０［Δｘ，Δｔ］であるから、この手法は１次の精度である。打切り誤差の最低次数が人工粘性項を含むため、ここでは簡単のため、以下に示すような、人工粘性項だけを含む２次の修正方程式を考えることにする。
【数２５】
2004246744.tif 000027
次に、（２３）式で、感度のパラメータ
【数２６】
2004246744.tif 000028
に関して、一般化感度解析法を適用する。そのために、まず、応答汎関数を次のように定義する。
【数２７】
2004246744.tif 000029
ここで、
【数２８】
2004246744.tif 000030
とする。また、（ｘｐ，ｔｐ）はモニタリングポイントに相当する。系のパラメータの変化がδαのとき、応答汎関数の変化（感度）は次式で与えられる。
【数２９】
2004246744.tif 000031
ここで、Ｆ’α，Ｆ’ｃはそれぞれＦのα、Ｃに関するＦｒｅｃｈｅｔ微分であり、本実施例では、
【数３０】
2004246744.tif 000032
となる。
【００１５】
支配方程式をパラメータで微分することにより、Ｄｉｆｆｅｒｅｎｔｉａｌ  Ｆｏｒｗａｒｄ  Ｓｙｓｔｅｍ  を導入する。すなわち、式（２３）をパラメータαで偏微分すると、次式となる。
【数３１】
2004246744.tif 000033
ここで、
【数３２】
2004246744.tif 000034
と定義し、φについて整理すると、
【数３３】
2004246744.tif 000035
となる。この式は、
【数３４】
2004246744.tif 000036
と表せる。ここで、微分系の演算子とソース項はそれぞれ、
【数３５】
2004246744.tif 000037
で与えられる。
（２９）式に適当な関数φ＊をかけて、対象とする時空間領域で積分する。
【数３６】
2004246744.tif 000038
以下、簡単のため次のように簡略化する。
【数３７】
2004246744.tif 000039
この表記では、（３２）式は以下のように書ける。
【数３８】
2004246744.tif 000040
ここで、（３２）式の左辺に部分積分を繰り返し適用することにより、（３２）式は、
【数３９】
2004246744.tif 000041
となる（Ｌ＊は形式的にＬに随伴な作用素である）。すなわち、（３０）式より、
【数４０】
2004246744.tif 000042
右辺のφが時間、または空間の導関数でなくなるように、部分積分を行うと、
【数４１】
2004246744.tif 000043
となる。ここで、Ｐ［φ＊，φ］は次式で与えられる。
【数４２】
2004246744.tif 000044
以上より、Ｌ＊は具体的に次式のように求まる。
【数４３】
2004246744.tif 000045
ここで、適当な関数φ＊として以下の方程式を満たすものを選ぶ。
【数４４】
2004246744.tif 000046
すると、φ＊はφの随伴関数であり、作用素Ｌ＊に関する基本解となる。（３３）式に、（３６）式を代入すると、
【数４５】
2004246744.tif 000047
となる。これより、（２５）式で定義される応答汎関数の感度は、
【数４６】
2004246744.tif 000048
となり、さらに（３１）式より、あるパラメータαに関する感度は、
【数４７】
2004246744.tif 000049
となる。したがって、一旦（３６）式を解いて随伴関数φ＊がわかれば、（３９）式より任意の感度が計算できる。また、系全体の感度δＲは
【数４８】
2004246744.tif 000050
と求めることができる。
【００１６】
最後に、随伴系の時空反転系（ξ，τ）を導入する。すなわち、
【数４９】
2004246744.tif 000051
これより、時間微分、空間微分は次のように与えられる。
【数５０】
2004246744.tif 000052
以上より、随伴演算子に対して、
【数５１】
2004246744.tif 000053
が成り立つ。よって、時空反転した随伴系の解は元の対象系と同様のコードで計算することができる。
【００１７】
具体的な数値解析において、感度解析を適用した例を、（１５）式の１次元の移流方程式に支配された波が移流することで説明する。
（１５）式の１次元の移流方程式で、初期値
【数５２】
2004246744.tif 000054
における初期値問題の厳密解は次のように与えられる。
【数５３】
2004246744.tif 000055
解析時間は０≦ｔ≦１０、解析体系は、図２に示すように、ｘ方向のみの１次元で、−５≦ｘ≦２５とした。初期条件は、ｘ＝０においてピークがあるような次式のｃｏｓ波を与えた。
【数５４】
2004246744.tif 000056
この初期値における濃度の分布が、上で述べたように波形を維持したまま伝わっていくのが厳密解であるが、数値解は離散化に伴う打切り誤差のため、波形が時間とともに拡散していく。
解析条件として、メッシュ幅Δｘ，時間刻み幅Δｔと、境界条件は、
【数５５】
2004246744.tif 000057
と設定した。
微分系、随伴系の初期条件・境界条件は、境界積分項Ｐを消去するように設定するので、以下のようになる。
【数５６】
2004246744.tif 000058
上述のメッシュ幅、時間刻み幅の条件等で解析した数値解と厳密解との比較を図３に、Δｘ，Δｔに対する感度を図４に示す。また、それらの結果により、厳密解を予測した結果を図５に示す。図３で、四角形で示したものが数値解で、直線が厳密解である。図５では、四角形で示したものが推定した厳密解であり、直線が厳密解である。
これらの図からわかるように、１次元移流方程式では、予測の基となる数値結果の数値条件近傍であるとき、予測結果が妥当である。
【００１８】
本発明は、自然界や産業・工業分野における現象の予測・評価に用いられる数値計算において、現象を支配する微分方程式（基礎方程式）を、差分法、有限要素法、境界要素法などの離散化手法を用いて離散化した代数方程式（離散化式）に対し、テイラー解析やフーリエ解析などの解析的方法により代数方程式から再び連続的な微分方程式（修正方程式）を求め、その修正方程式に対して、時間刻み幅や空間格子幅などの数値計算上のパラメータに関する感度（パラメータが変化した場合に結果がどの程度変化するか）を求めることにより、数値計算結果の数値パラメータに関する感度を評価する。その結果、以下のような評価が可能となる。
（１）感度解析の結果を用いた、他の計算パラメータを用いて計算した場合の結果（たとえばより細かい空間格子幅を用いて解析した場合の計算結果）の推定
（２）その結果を用いた厳密解の推定
（３）その結果を用いた数値解析解に含まれている数値誤差の定量的な推定
【００１９】
（１）に関しては、従来の数値計算法では、空間格子幅を変更しようとした場合、メッシュ生成を再実行し、入力データを作成しなおし、もう一度数値計算を実行する必要があった。さらに、空間格子幅を細かくする場合、計算点の数が増え、また、時間刻み幅を細かくする場合、ある時間に到達するまでの計算ステップが増え、計算時間の増大をまねく。本発明では、そうした作業や数値計算の再実行をすることなく、空間格子幅や時間刻み幅を変更した場合の数値計算結果を推定できる。
（２）に関しては、推定した厳密解を、より精度の良い数値解析結果として利用することも可能である。
（３）に関しては、（２）で推定した厳密解と数値解析結果の差をとることにより、誤差の定量的な推定（数値解の信頼性の推定）が可能となる。
なお、本発明は、上記した実施の形態に限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々変更を加え得ることができる。
【００２０】
【発明の効果】
上述したように、本発明では、自然現象を表す数学的モデルの数値解析結果の感度解析を行なうことにより、計算パラメータの依存性を定量的に評価でき、数値解の信頼性の推定や厳密解の推定を行なうことが可能となった。
【図面の簡単な説明】
【図１】実施形態の感度解析処理を示すフローチャートである。
【図２】実施例の１次元移流方程式の例を示した図である。
【図３】数値解と厳密解との比較を示す図である。
【図４】メッシュ幅，時間刻み幅に対する感度を示す図である。
【図５】厳密解と、予測した厳密解との比較を示す図である。
【図１】
2004246744.tif 000059
【図２】
2004246744.tif 000060
【図３】
2004246744.tif 000061
【図４】
2004246744.tif 000062
【図５】
2004246744.tif 000063