(57)【特許請求の範囲】
【請求項１】
空間上の複数の定義点それぞれにおける１種類以上の流体に関する物理的な変数の値を、流体の物理現象を表す偏微分方程式を含む支配方程式と、前記偏微分方程式と局所的に近似する３次多項式であり前記定義点ごとに係数が決定される補間関数と、を用いて、所定の微小時間単位で更新する数値計算装置による数値計算方法であって、
前記数値計算装置は、記憶部と処理部とを有し、
前記記憶部は、前記支配方程式と、前記補間関数と、前記複数の定義点と、前記定義点ごとの前記変数の値と、を記憶しており、
前記定義点のうちの３点を、ｊ−１点、ｊ点、ｊ＋１点とし、基準点をｊ点とし、ｊ−１点側を前記流体が移動してくる元の方向である風上とし、ｊ＋１点側を風上と反対の方向である風下とし、
前記各定義点ｊ−１点、ｊ点、ｊ＋１点について、それらの座標をｘｊ−Δｘ、ｘｊ、ｘｊ＋Δｘとし、それらの前記定義点ごとの変数のうちの１種類の変数の値をｆｊ−１，ｆｊ，ｆｊ＋１とし、
前記定義点ｊ−１点からｊ点までの当該変数の積分値をｘｆｊ−１／２、前記各定義点ｊ点からｊ＋１点までの当該変数の積分値をｘｆｊ＋１／２とし、
前記補間関数を、次の式（６）とした場合、
Ｆ（ｘ−ｘｊ）＝ｃ３（ｘ−ｘｊ）３＋ｃ２（ｘ−ｘｊ）２＋ｃ１（ｘ−ｘｊ）＋ｃ０
・・・式（６）
前記処理部は、
前記記憶部を参照し、前記基準点であるｊ点について、前記式（６）で表した補間関数を用いる際、前記係数ｃ１として、次の式（８）で表す値またはそれに準ずる値を使用する
ことを特徴とする数値計算方法。
【数１】
0005500637.tif 000019
【請求項２】
前記処理部は、
前記記憶部を参照し、前記基準点であるｊ点について、前記式（６）で表した補間関数を用いる際、前記係数ｃ３、ｃ２、ｃ０として、次の式（１０）で表す値またはそれらに準ずる値を使用する
ことを特徴とする請求項１に記載の数値計算方法。
【数２】
0005500637.tif 000020
【請求項３】
風上と風下が入れ替わり、ｊ−１点側が風下となり、ｊ＋１点側を風上となった場合、
前記処理部は、
前記記憶部を参照し、前記基準点であるｊ点について、前記式（６）で表した補間関数を用いる際、前記係数ｃ１として、次の式（１１）で表すｃ１’の値またはそれに準ずる値を使用する
ことを特徴とする請求項１に記載の数値計算方法。
【数３】
0005500637.tif 000021
【請求項４】
前記処理部は、
前記記憶部を参照し、前記基準点であるｊ点について、前記式（６）で表した補間関数を用いる際、前記係数ｃ３、ｃ２、ｃ０として、次の式（１２）で表すｃ３’、ｃ２’、ｃ０’の値またはそれらに準ずる値を使用する
ことを特徴とする請求項３に記載の数値計算方法。
【数４】
0005500637.tif 000022
【請求項５】
空間上の複数の定義点それぞれにおける１種類以上の流体に関する物理的な変数の値を、流体の物理現象を表す偏微分方程式を含む支配方程式と、前記偏微分方程式と局所的に近似する３次多項式であり前記定義点ごとに係数が決定される補間関数と、を用いて、所定の微小時間単位で更新する数値計算装置による数値計算方法であって、
前記数値計算装置は、記憶部と処理部とを有し、
前記記憶部は、前記支配方程式と、前記補間関数と、前記複数の定義点と、前記定義点ごとの前記変数の値と、を記憶しており、
前記定義点のうちの３点を、ｊ−１点、ｊ点、ｊ＋１点とし、基準点をｊ点とし、ｊ−１点側を前記流体が移動してくる元の方向である風上とし、ｊ＋１点側を風上と反対の方向である風下とし、
前記各定義点ｊ−１点、ｊ点、ｊ＋１点について、それらの座標をｘｊ−Δｘ、ｘｊ、ｘｊ＋Δｘとし、それらの前記定義点ごとの変数のうちの１種類の変数の値をｆｊ−１，ｆｊ，ｆｊ＋１とし、
前記定義点ｊ−１点からｊ点までの当該変数の積分値をｘｆｊ−１／２、前記各定義点ｊ点からｊ＋１点までの当該変数の積分値をｘｆｊ＋１／２とし、
前記補間関数を、次の式（６）とした場合、
Ｆ（ｘ−ｘｊ）＝ｃ３（ｘ−ｘｊ）３＋ｃ２（ｘ−ｘｊ）２＋ｃ１（ｘ−ｘｊ）＋ｃ０
・・・式（６）
前記処理部は、
前記記憶部を参照し、前記基準点であるｊ点について、前記式（６）で表した補間関数を用いる際、前記係数ｃ１として、次の式（１３）で表す値またはそれに準ずる値を使用する
ことを特徴とする数値計算方法。
【数５】
0005500637.tif 000023
ただし、（√１７−３)／４＜α≦１／２
【請求項６】
空間上の複数の定義点それぞれにおける１種類以上の流体に関する物理的な変数の値を、流体の物理現象を表す偏微分方程式を含む支配方程式と、前記偏微分方程式と局所的に近似する３次多項式であり前記定義点ごとに係数が決定される補間関数と、を用いて、所定の微小時間単位で更新する数値計算装置による数値計算方法であって、
前記数値計算装置は、記憶部と処理部とを有し、
前記記憶部は、前記支配方程式と、前記補間関数と、前記複数の定義点と、前記定義点ごとの前記変数の値と、を記憶しており、
前記定義点のうちの３点を、ｊ−１点、ｊ点、ｊ＋１点とし、基準点をｊ点とし、ｊ−１点側を前記流体が移動してくる元の方向である風上とし、ｊ＋１点側を風上と反対の方向である風下とし、
前記各定義点ｊ−１点、ｊ点、ｊ＋１点について、それらの座標をｘｊ−Δｘ、ｘｊ、ｘｊ＋Δｘとし、それらの前記定義点ごとの変数のうちの１種類の変数の値をｆｊ−１，ｆｊ，ｆｊ＋１とし、
前記定義点ｊ−１点からｊ点までの当該変数の積分値をｘｆｊ−１／２、前記各定義点ｊ点からｊ＋１点までの当該変数の積分値をｘｆｊ＋１／２とし、
前記補間関数を、次の式（６）とした場合、
Ｆ（ｘ−ｘｊ）＝ｃ３（ｘ−ｘｊ）３＋ｃ２（ｘ−ｘｊ）２＋ｃ１（ｘ−ｘｊ）＋ｃ０
・・・式（６）
前記処理部は、
前記記憶部を参照し、前記基準点であるｊ点について、前記式（６）で表した補間関数を用いる際、前記係数ｃ１として、次の式（１７）で表す値またはそれに準ずる値を使用する
ことを特徴とする数値計算方法。
【数６】
0005500637.tif 000024
ただし、０＜α＜１
【請求項７】
請求項１から請求項６のいずれか一項に記載の数値計算方法をコンピュータに実行させるためのプログラム。
【請求項８】
請求項７に記載のプログラムを記録した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、偏微分方程式の演算に使用する数値計算方法に関する。
【背景技術】
【０００２】
近年、コンピュータを用いた流体の物理現象のシミュレーションに関する研究や開発が、活発になってきている。シミュレーションの対象は、例えば、液体、気体、熱などによる様々な流体運動である。
【０００３】
そのような数値計算を行う場合、流体の変数（速度、圧力、温度、密度など）に関する偏微分方程式を、離散化された空間および時間に関して解く、という手法が用いられるのが一般的である。なぜなら、偏微分方程式の数学的な解析解（厳密解）を求めることは現実的に不可能だからである。したがって、高精度の数値計算方法を開発することが重要となる。
【０００４】
そして、より高精度な計算を行うため、変数の値だけでなく、変数の微分係数（空間勾配）の値または区間積分平均値（以下、単に「積分値」ともいう。）も従属変数とするマルチモーメント手法であるＣＩＰ法（Cubic-Interpolated Pseudo-Particle Scheme）やＩＤＯ法（Interpolated Differential Operator Scheme：局所補間微分オペレータ法）が開発された。
【０００５】
ＩＤＯ法は、物理的に保存しなければならない質量、運動量、エネルギーを厳密に保存する保存形ＩＤＯ法と、高精度数値計算の範囲内で物理的な保存を満足する非保存形ＩＤＯ法とに大別することができる。保存形ＩＤＯ法は、例えば、変数の積分値を従属変数とすることで、その変数に関する物理量を保存することができる。保存形ＩＤＯ法と非保存形ＩＤＯ法は、それぞれ一長一短を有し、場面や状況に応じて選択および使用される。
【０００６】
また、本出願人は、非特許文献１および非特許文献２（非特許文献１とほぼ同内容）において、保存形ＩＤＯ法について、偏微分方程式と局所的に近似する補間関数を用いた計算手法を提案した。保存形ＩＤＯ法では、偏微分方程式中の微分項を補間関数の微分値に置換することで空間の離散化を行う。
【０００７】
非特許文献１および非特許文献２に開示された技術は、保存形ＩＤＯ法において、格子点（定義点）上の変数の値と積分値を用いた補間関数をベースとする高精度数値計算手法である。この手法では、圧縮性流体解析などで用いられる風上補間を行う。風上補間とは、着目する格子点（以下、「基準点」ともいう。）よりも風上（流体が移動してくる元の方向）の情報を用いて行う補間である。つまり、基準点を格子点であるｊとし、風上方向の隣接格子点をｊ−１とした場合、ｊ点上の変数の値、ｊ−１点上の値、および、それらの間の積分値を使って、２次多項式である補間関数の係数を決定していた。
【０００８】
これを、図１２を参照して説明する。保存形ＩＤＯ法において、空間的な間隔がΔｘである３点、ｊ−１点、ｊ点、ｊ＋１点について考える。なお、基準点をｊ点とし、ｊ−１点側を風上とし、ｊ＋１点側を風下（風上の反対側）とする。
【０００９】
ｊ点の座標をｘｊとすると、ｊ−１点とｊ＋１点の座標はそれぞれｘｊ−Δｘ、ｘｊ＋Δｘとなる。また、ｊ−１点、ｊ点、ｊ＋１点における変数の値をｆｊ−１，ｆｊ，ｆｊ＋１とする。また、ｊ−１点からｊ点までの積分値をｘｆｊ−１／２、ｊ点からｊ＋１点までの積分値をｘｆｊ＋１／２とする。
【００１０】
ここで、例えば、ｊ−１点からｊ点までの範囲に関する２次多項式の補間関数
Ｆ（ｘ−ｘｊ）＝ａ（ｘ−ｘｊ）２＋ｂ（ｘ−ｘｊ）＋ｃの係数（ａ，ｂ，ｃ）を決定する場合を想定する。補間関数が２次多項式なので、拘束条件が３つあれば補間関数の係数を決定できる。そこで、ｆｊ−１、ｘｆｊ−１／２、ｆｊの３つの値を使う。拘束条件は、次の式（１）の通りである。
【数１】
0005500637.tif 000002
なお、∫ｘｊｘｊ−ΔｘＦ（ｘ−ｘｊ）ｄｘは、Ｆ（ｘ−ｘｊ）をｘｊ−Δｘからｘｊまでの区間について積分した値であることを示す（以下同様）。
したがって、（１／Δｘ）∫ｘｊｘｊ−ΔｘＦ（ｘ−ｘｊ）ｄｘは、区間積分平均値（その区間における前記積分した値の平均値）を示す。
【００１１】
式（１）により、補間関数Ｆ（ｘ−ｘｊ）の係数（ａ，ｂ，ｃ）の値を決めることができる。これらの係数のうち、特に重要なのはｂである。補間関数Ｆ（ｘ−ｘｊ）を一階微分した∂Ｆ（ｘ−ｘｊ）／∂ｘにｘ＝ｘｊを代入したときの値がｂであり、この値が数値計算全体に与える影響が大きいからである。ｂの値は、具体的には、次の式（２）となる。
【数２】
0005500637.tif 000003
【００１２】
ここで、式（２）の精度（誤差）を調べるために、式（２）の右辺にテイラー展開した離散変数を代入すると、以下の式（３）となる。
【００１３】
【数３】
0005500637.tif 000004
【００１４】
つまり、式（２）と式（３）を比較すればわかるように、式（２）を数値計算に用いた際の誤差は、式（３）の右辺第２項以下である。すなわち、計算誤差はΔｘ２オーダーであり、また、誤差の最上位の項における係数の数字部分は１／１２である。
【先行技術文献】
【非特許文献】
【００１５】
【非特許文献１】Y. Imai, T. Aoki and K. Takizawa, “Conservative form of interpolated differential operator scheme for compressible and incompressible fluid dynamics”, Journal of Computational Physics, the Kingdom of the Netherlands, ELSEVIER, 2008, Vol. 227, Issue 4, p.2263-2285
【非特許文献２】青木尊之、今井陽介、“保存形ＩＤＯ（局所補間微分オペレータ）法”、応用数理学会誌・岩波「応用数理」、岩波書店、２００８年６月、Vol.18、No.2、p.32-45
【発明の概要】
【発明が解決しようとする課題】
【００１６】
前記した従来の手法（式（２）を用いる手法）では、ある程度の計算精度を得ることはできるが、数値粘性（流体力学等の計算において最も重要な誤差であり、空間的に細かいプロファイルを消してしまう数値計算の誤差）をより多く含んでいて、長時間の計算を行うと数値結果がぼやけた状態（空間的に急峻な勾配に関する誤差が大きくなった状態）になるという問題点がある。
【００１７】
この問題の原因としては、例えば、補間関数の係数を決定する際に基準点よりも風下の情報をまったく使っていないことが考えられる。しかし、例えば、補間関数の次数を１つ上げて３次多項式とし、風下の積分値（図１２の例では「ｘｆj+１/2」）も用いて拘束条件を４つとする方法では、高波数（空間的に急峻な勾配）の領域で計算の振動が多く実用的ではない。また、基準点から遠い箇所の情報を使うことは、その境界条件の設定等が困難になったり、計算速度の遅延が起きたりするため、好ましくない。したがって、風下の情報のうち適切な情報を用いて補間関数の係数を決定し、高精度な計算を行う手法が求められていた。
【００１８】
そこで、本発明は、前記状況に鑑みてなされたものであり、偏微分方程式の演算に使用する数値計算方法において、風上の情報に加えて、風下の情報のうち適切な情報も用いて係数を決定した補間関数を使用することで、高精度な計算を実現することを課題とする。
【課題を解決するための手段】
【００１９】
前記課題を解決するために、本発明では、偏微分方程式の補間関数を３次多項式とし、基準点の風下側の近傍に仮想点を想定し、その仮想点における変数の値も拘束条件に使用することで、補間関数の係数を決定する。
【発明の効果】
【００２０】
本発明によれば、偏微分方程式の演算に使用する数値計算方法において、風上の情報に加えて、風下の情報のうち適切な情報も用いて係数を決定した補間関数を使用することで、高精度な計算を実現することができる。
【図面の簡単な説明】
【００２１】
【図１】本実施形態における数値計算装置の構成図である。
【図２】本実施形態における数値計算の全体的な流れを示すフローチャートである。
【図３】第１実施形態における第２の補間関数の決定のプロセスの説明図である。
【図４】第１実施形態における第２の補間関数の決定のための演算の流れを表すフローチャートである。
【図５】実施例１における数値計算結果の例を示すグラフである。
【図６】実施例１における数値計算結果の例を示すグラフであり、（ａ）は振幅比を示すグラフであり、（ｂ）は位相を示すグラフである。
【図７】実施例１における数値計算結果の例を示すグラフであり、（ａ）は振幅比誤差を示すグラフであり、（ｂ）は位相誤差を示すグラフである。
【図８】実施例２における数値計算結果の例を示すグラフであり、（ａ）は従来技術（２次風上）による計算結果の例を示すグラフであり、（ｂ）は本実施形態の数値計算（３次風上）による計算結果の例を示すグラフである。
【図９】実施例２における数値計算結果の例を示すグラフであり、（ａ）は従来技術（２次風上）による数値計算結果の例を示すグラフであり、（ｂ）は本実施形態（３次風上）による数値計算結果の例を示すグラフである。
【図１０】第２実施形態における第２の補間関数の決定のプロセスの説明図である。
【図１１】第３実施形態における第２の補間関数の決定のプロセスの説明図である。
【図１２】従来技術における補間関数の決定のプロセスの説明図である。
【発明を実施するための形態】
【００２２】
以下、本発明を実施するための形態（以下、「実施形態」という。）に関し、第１実施形態、第２実施形態および第３実施形態（それらの総称を「本実施形態」という。）について、適宜図面を参照しながら説明する。本実施形態では、保存形ＩＤＯ法の場合を例にとって説明する。まず、図１を参照しながら、本実施形態における数値計算方法を実行する数値計算装置の構成について説明する。
【００２３】
図１に示すように、数値計算装置１０は、入力部１、記憶部２、メモリ３、演算部４および出力部５を備えて構成されるコンピュータ装置であり、例えば、パソコン、ワークステーション、スーパーコンピュータなどにより実現することができる。
【００２４】
入力部１は、数値計算装置１０の使用者（以下、単に「使用者」という。）が各種情報を入力する手段であり、例えば、キーボード、マウスなどである。
記憶部２は、各種情報を記憶する手段であり、例えば、ＲＯＭ（Read Only Memory）、ハードディスクなどである。記憶部２は、ここでは、流体等の数値計算（シミュレーション）を行うために、数値計算を行う場所を示す定義点である格子点、物理量（変数）やその微分係数の初期条件、支配方程式（流体の物理現象を表わす偏微分方程式など）、補間関数（第１の補間関数、第２の補間関数）、各種プログラムなどを記憶している。
【００２５】
メモリ３は、演算部４が使用する一時記憶手段であり、例えば、ＲＡＭ（Random Access Memory）などである。
演算部４は、記憶部２に記憶された各種情報（プログラム、データなど）などに基づいて演算（処理）を行うものであり、例えば、ＣＰＵ（Central Processing Unit）である。なお、演算部４は、ここでは１つとしているが、複数設けることにより並列計算を行うようにしてもよい。
出力部５は、演算部４による演算の結果を出力するものであり、例えば、ディスプレイなどである。
【００２６】
次に、図２を参照しながら、本実施形態における数値計算の全体的な流れについて説明する（適宜図１参照）。
図２に示すように、まず、使用者が入力部１によって数値計算の開始に関する操作を行うと、演算部４は、記憶部２に記憶された情報に基づいて、演算を行う場所を示す定義点である格子点、各物理量の初期条件などを設定する（ステップＳ１）。格子点は、例えば、１次元の場合は「５１２」、２次元の場合は「１２８×１２８」、３次元の場合は「１６×１６×８」、といったように設定される。また、各物理量の初期条件としては、それぞれの格子点における変数の値、および、それらの値の微分係数（空間勾配）や積分値などが設定される。変数としては、例えば、流体の速度、圧力、温度、密度などが挙げられる。
【００２７】
次に、演算部４は、記憶部２に記憶された情報に基づいて、時間と空間に関する偏微分方程式などからなる支配方程式、各物理量の値を修正するための補間関数などを設定する（ステップＳ２）。なお、支配方程式としては、例えば、圧縮性流体の運動方程式、流体の質量保存方程式などが挙げられる。また、このステップＳ２で補間関数を設定する処理が本発明の特徴であるが、その詳細な説明は後記する。
【００２８】
続いて、演算部４は、ステップＳ１およびステップＳ２で設定した各種情報を用い、時刻「ｔ」における各格子点での各物理量の値に基づいて、時刻「ｔ＋１」における各格子点での各物理量の値を算出する、という演算をｔ＝０,１,２,・・・ｎに関して行う（ステップＳ３〜ステップＳ５）。
【００２９】
次に、演算部４は、ステップＳ３〜ステップＳ５で算出した各物理量の値を、出力部５に出力する（ステップＳ６）。以上で、数値計算が完了する。
なお、ステップＳ１〜Ｓ６の処理のうち、ステップＳ２における補間関数の設定の処理については以下で説明するが、それ以外の処理については、本出願人が先に出願した特許出願に係る特許第４３０４２７７号公報および特開２００８−１９７９４９号公報において記述されているので、ここでは詳細な説明を省略する。
【００３０】
＜第１実施形態＞
次に、第１実施形態として、前記したステップＳ２における補間関数の設定の処理について説明する。まず、図３を参照して、補間関数の設定のイメージについて説明する。ここでは、第１の補間関数を用いて１つの値を算出し、その値と他の３つの値を用いて第２の補間関数を設定（係数を決定）する場合について説明する。なお、第２の補間関数が設定対象の補間関数であり、第１の補間関数は第２の補間関数を作成するための手段である。
【００３１】
図３に示すように、保存形ＩＤＯ法において、空間的な間隔がΔｘである３点、ｊ−１点、ｊ点、ｊ＋１点について考える。なお、基準点をｊ点とし、ｊ−１点側を風上とし、ｊ＋１点側を風下とする。
【００３２】
ｊ点の座標をｘｊとすると、ｊ−１点とｊ＋１点の座標はそれぞれｘｊ−Δｘ、ｘｊ＋Δｘとなる。また、ｊ−１点、ｊ点、ｊ＋１点における変数の値をｆｊ−１，ｆｊ，ｆｊ＋１とする。また、ｊ−１点からｊ点までの積分値をｘｆｊ−１／２、ｊ点からｊ＋１点までの積分値をｘｆｊ＋１／２とする。
【００３３】
このとき、まず、ｊ点からｊ＋１点までの範囲に関する第１の補間関数
Ｈ（ｘ−ｘｊ）＝ａ（ｘ−ｘｊ）２＋ｂ（ｘ−ｘｊ）＋ｃ（・・・式（４ａ））の係数（ａ，ｂ，ｃ）を決定する。第１の補間関数は２次多項式なので、拘束条件が３つあればこの第１の補間関数の係数を決定できる。そこで、ｆｊ、ｘｆｊ＋１／２、ｆｊ＋１の３つの値を使う。拘束条件は、次の式（４ｂ）の通りである。
【００３４】
【数４】
0005500637.tif 000005
【００３５】
式（４ｂ）により、第１の補間関数Ｈ（ｘ）の係数（ａ，ｂ，ｃ）の値が決まる（係数（ａ，ｂ，ｃ）の値の記載は省略）。この第１の補間関数Ｈ（ｘ−ｘｊ）を用いて、ｆｊ＋１／２の値を次の式（５）のように求めることができる。つまり、この第１の補間関数Ｈ（ｘ−ｘｊ）は、ｆｊ＋１／２の値を求めるための手段である。
【００３６】
【数５】
0005500637.tif 000006
【００３７】
つまり、格子点ではないｊ＋１／２点（ｊ点とｊ＋１の中間点）を仮想点とし、その仮想点上のｆの値を、第１の補間関数を用いて算出するのである。
【００３８】
次に、ｊ−１点からｊ＋１／２点までの範囲に関する次の式（６）に示す第２の補間関数の係数（ｃ３、ｃ２、ｃ１、ｃ０）を決定する。
Ｆ（ｘ−ｘｊ）＝ｃ３（ｘ−ｘｊ）３＋ｃ２（ｘ−ｘｊ）２＋ｃ１（ｘ−ｘｊ）＋ｃ０
・・・式（６）
【００３９】
第２の補間関数は３次多項式なので、拘束条件が４つあればこの第２の補間関数の係数を決定できる。そこで、ｆｊ−１、ｘｆｊ−１／２、ｆｊ、ｆｊ＋１／２の４つの値を使う。拘束条件は、次の式（７）の通りである。
【００４０】
【数６】
0005500637.tif 000007
【００４１】
これらにより、第２の補間関数Ｆ（ｘ）の係数（ｃ３、ｃ２、ｃ１、ｃ０）の値が決まる。これらの係数のうち、特に重要なのはｃ１ある。第２の補間関数Ｆ（ｘ−ｘｊ）を一階微分した∂Ｆ（ｘ−ｘｊ）／∂ｘにｘ＝ｘｊを代入したときの値がｃ１であり、この値が数値計算全体に与える影響が大きいからである。ｃ１の値は、具体的には、次の式（８）となる。
【００４２】
【数７】
0005500637.tif 000008
【００４３】
ここで、式（８）の精度（誤差）を調べるために、式（８）の右辺にテイラー展開した離散変数を代入すると、以下の式（９）のようになる。
【数８】
0005500637.tif 000009
【００４４】
つまり、式（８）と式（９）を比較すればわかるように、式（８）を数値計算に用いた際の誤差は、式（９）の右辺第２項以下である。すなわち、計算誤差はΔｘ３オーダーであり、また、誤差の最上位の項における係数の数字部分は１／９０である。前記したように、従来技術における２次多項式の補間関数においては（式（２）、（３）参照）、計算誤差がΔｘ２オーダーであり、また、誤差の最上位の項における係数の数字部分は１／１２であったことと比べると、式（８）を用いた場合の計算精度が従来技術の場合よりも有意に高いことがわかる。換言すれば、式（３）の右辺第２項以下と、式（９）の右辺第２項以下とを比較すると、それぞれの最上位の項における微分の階数が違うので厳密な比較は容易ではないが、それでも前記したようにオーダーと係数の違いから、式（２）を用いた場合の計算誤差のほうが明らかに大きいと言える。
【００４５】
なお、同様にして、その他の係数の値は、次の式（１０）のようになる。
【数９】
0005500637.tif 000010
【００４６】
また、風上と風下が入れ替わった場合、仮想点をｊ−１／２点（ｊ−１点とｊの中間点）として、前記と同様に計算することができる。そのとき、第２の補間関数において、ｃ１に対応する係数であるｃ１’は、次の式（１１）に示す通りとなる。
【数１０】
0005500637.tif 000011
【００４７】
なお、同様にして、ｃ３、ｃ２、ｃ０に対応する係数であるｃ３’、ｃ２’、ｃ０’は、その他の係数の値は、次の式（１２）のようになる。
【数１１】
0005500637.tif 000012
【００４８】
次に、図４を参照しながら、第１実施形態における第２の補間関数の係数の決定のための演算（処理）（図２のステップＳ２における補間関数の設定のための処理）について説明する（適宜図１、図３参照）。数値計算装置１０の演算部４は、記憶部２とメモリ３を使用しながら各演算を行う。
【００４９】
まず、演算部４は、ｆｊ、ｘｆｊ＋１／２、ｆｊ＋１の３つの値を取得する（ステップＳ２１）。
次に、演算部４は、その取得した３つの値を用いて、第１の補間関数の係数を決定する（ステップＳ２２）。
【００５０】
続いて、演算部４は、係数が決定した第１の補間関数を用いて、ｆｊ＋１／２の値を算出する（ステップＳ２３）。
次に、演算部４は、ｆｊ、ｘｆｊ＋１／２、ｆｊ＋１、ｆｊ＋１／２の４つの値を取得する（ステップＳ２４）。
【００５１】
最後に、演算部４は、その取得した４つの値を用いて、第２の補間関数の係数を決定する（ステップＳ２５）。
【００５２】
このように、第１実施形態によれば、数値計算装置１０による偏微分方程式の演算に使用する数値計算方法において、従来から風上補間のときに使用していたｆｊ、ｘｆｊ＋１／２、ｆｊ＋１の３つの値に加えて、第１の補間関数を用いて算出した仮想点ｊ＋１／２点における値ｆｊ＋１／２を用いて第２の補間関数の係数を決定する。したがって、数値計算装置１０は、風上の情報に加えて、風下の適切な情報も用いて係数を決定した第２の補間関数を使用することで、高精度な計算を実現することができる。
【実施例１】
【００５３】
（フーリエ解析）
次に、図５〜図７を参照して、第１実施形態の数値計算方法の実施例に関し、フーリエ解析を用いた精度と安定性の検証について説明する。この実施例１は、線形移流方程式に対してフーリエ変換した離散変数を代入してフーリエ解析を行い、波数と、クーラン数、位相、振幅比とのそれぞれの関係をグラフ化したものである。なお、時間積分には４次精度ルンゲクッタ法を用いた。
【００５４】
この実施例１については、前記した非特許文献１における「3.Fourier analysis(p.2272-2276)」と、前記した非特許文献２における「5 Fourier解析（ｐ．３８−４１）」とに、詳細な記述があるので、ここでは詳細な説明を省略する。それらの非特許文献１，２における記述と比較して、この実施例１における相違点は、補間関数が３次多項式であることと、その３次多項式である補間関数の係数の決定方法とだけである。
【００５５】
図５に示すグラフでは、縦軸をクーラン数とし、横軸を波数とした座標上に、偏微分方程式の時間発展を記述した行列の固有値を表示している。このようなグラフでは、固有値「１．０」を示すラインが上方に位置するほど数値計算が安定することを表す。ここでは、固有値「１．０」を示すラインＬ５１がかなり上方に位置していることで、第１実施形態による数値計算方法が安定していることがわかる。
【００５６】
図６（ａ）に示すグラフでは、縦軸を振幅比とし、横軸を波数として、厳密解と、保存形ＩＤＯ法（２次風上：従来手法である２次多項式の補間関数を使った風上補間による手法）による値と、保存形ＩＤＯ法（３次風上：第１実施形態における３次多項式の補間関数を使った風上補間（ここでは風下の情報も使っているが「風上補間」と称する。）による手法）による値と、を表示している。
【００５７】
このようなグラフでは、当然、厳密解に近いほど数値計算の精度が高いことになる。そして、図６（ａ）を見れば、第１実施形態の手法である保存形ＩＤＯ法（３次風上）は、従来手法である保存形ＩＤＯ法（２次風上）と比べて、厳密解に有意に近く、振幅比についての数値計算の精度が高いことがわかる。
【００５８】
図６（ｂ）に示すグラフでは、縦軸を位相とし、横軸を波数として、厳密解と、保存形ＩＤＯ法（中心：従来手法である四次関数の補間関数を使った中心補間（３つの格子点のうちの真ん中の格子点を基準にして行う補間）による手法）による値と、保存形ＩＤＯ法（２次風上）による値と、保存形ＩＤＯ法（３次風上）による値と、を表示している。
【００５９】
ここでは、保存形ＩＤＯ法（中心）と、保存形ＩＤＯ法（２次風上）と、保存形ＩＤＯ法（３次風上）との間で、厳密解との差に関して大きな差異はない。つまり、第１実施形態の手法である保存形ＩＤＯ法（３次風上）は、位相についての数値計算の精度に関しては、他の従来手法と特に差がないことがわかる。
【００６０】
さらに、図７（ａ）に示すグラフでは、縦軸を振幅比誤差とし、横軸を波数として、保存形ＩＤＯ法（２次風上）による値と、保存形ＩＤＯ法（３次風上）による値と、を表示している。これを見れば、第１実施形態の手法である保存形ＩＤＯ法（３次風上）は、従来手法である保存形ＩＤＯ法（２次風上）と比べて、振幅比誤差が有意に小さく、振幅比についての数値計算の精度が高いことがわかる。
【００６１】
図７（ｂ）に示すグラフでは、縦軸を位相誤差とし、横軸を波数として、保存形ＩＤＯ法（中心）による値と、保存形ＩＤＯ法（２次風上）による値と、保存形ＩＤＯ法（３次風上）による値と、を表示している。
【００６２】
ここでは、保存形ＩＤＯ法（中心）と、保存形ＩＤＯ法（２次風上）と、保存形ＩＤＯ法（３次風上）との間で、大きな差異はない。つまり、第１実施形態の手法である保存形ＩＤＯ法（３次風上）は、位相についての数値計算の精度に関しては、他の従来手法と特に差がないことがわかる。
【００６３】
このように、実施例１からわかるように、第１実施形態の数値計算方法は、位相についての計算精度は従来技術とほぼ同等であるものの、固有値や振幅比についての計算精度は従来技術よりも有意に高く、総合的な計算精度が有意に高い。
【実施例２】
【００６４】
（圧縮性レーリー・テーラー不安定性の成長過程シミュレーション）
次に、図８、図９を参照して、第１実施形態の数値計算方法を圧縮性レーリー・テーラー不安定性の成長過程シミュレーションに用いた場合の実施例について説明する。この実施例２は、圧縮性流体の計算例として、重い（密度の高い）流体が軽い（密度の低い）流体の上に置かれている場合の、２つの流体の界面に与えられた擾乱がレーリー・テーラー不安定性により成長する過程をシミュレーションしたときの２つの流体をグラフ化したものである。なお、そのときの２つの流体について、図８では横から見た様子（水平方向に見た様子）を、また、図９では上から見た様子（鉛直下向きに見た様子）を、それぞれ示している。また、ここでは、圧縮性流体計算における移流項計算部分に第１実施形態の数値計算方法を適用した。
【００６５】
この実施例２については、前記した非特許文献２における「６  計算例  ６．１  圧縮性レーリー・テーラー不安定(p.41-43)」に詳細な記述があるので、ここでは詳細な説明を省略する。非特許文献２における記述と比較して、この実施例２における相違点は、補間関数が３次多項式であることと、その３次多項式である補間関数の係数の決定方法とだけである。
【００６６】
図８（ａ）に示すグラフは、保存形ＩＤＯ法（２次風上）により計算を行った場合の２つの流体の様子を表したものである。図８（ｂ）に示すグラフは、保存形ＩＤＯ法（３次風上）により計算を行った場合の２つの流体の様子を表したものである。いずれのグラフも解像度は２５６×２５６である。
【００６７】
これらを見れば、図８（ａ）に示すグラフでは細部がぼやけているのに比べて、図８（ｂ）に示すグラフでは細部まで明確に表示できていることがわかる。つまり、第１実施形態の手法である保存形ＩＤＯ法（３次風上）は、従来手法である保存形ＩＤＯ法（２次風上）と比べて、数値計算の精度が高いことがわかる。
【００６８】
同様に、図９（ａ）に示すグラフは、保存形ＩＤＯ法（２次風上）により計算を行った場合の２つの流体の様子を表したものである。図９（ｂ）に示すグラフは、保存形ＩＤＯ法（３次風上）により計算を行った場合の２つの流体の様子を表したものである。いずれのグラフも解像度は１２８×１２８である。
【００６９】
これらを見れば、図９（ａ）に示すグラフでは細部がぼやけているのに比べて、図９（ｂ）に示すグラフでは細部まで明確に表示できていることがわかる。つまり、第１実施形態の手法である保存形ＩＤＯ法（３次風上）は、従来手法である保存形ＩＤＯ法（２次風上）と比べて、数値計算の精度が高いことがわかる。
【００７０】
このように、実施例２からわかるように、第１実施形態の数値計算方法は、従来技術よりも計算精度が高い。その理由について説明すると、次のようになる。つまり、このようなシミュレーションでは、初期の線形成長γ＝√ｇｋａ（ｇ：重力加速度、ｋ：擾乱の波数、ａ＝（ρＨ−ρＬ）／（ρＨ＋ρＬ）、ρＨ：重い流体の密度、ρＬ：軽い流体の密度）を過ぎると、ケルビン・ヘルムホルツ不安定性を伴った非線形成長に移り、長波長の波の成長率が大きくなり、非常に複雑な界面を形成しながらFree Fallと呼ばれる過程（ｈ〜0.06ａｇｔ２、ｈ：落下速度、ｔ：時間）へと移る。
【００７１】
そして、数値粘性が大きい従来手法（保存形ＩＤＯ法（２次風上））で計算した場合、波数の大きい波が数値的に減衰してしまい、Free Fallへの遷移や落下加速度が遅くなってしまうことがある。一方、第１実施形態の数値計算方法（保存形ＩＤＯ法（３次風上））では、そういった遅延が少なく、高精度で計算を行うことができる。
【００７２】
第１実施形態の数値計算方法では、第１の補間関数によってｆｊ＋１／２の値を算出し、そのｆｊ＋１／２の値も拘束条件に使用して第２の補間関数の係数を決定した。これにより、補間関数の次数が１つ増えただけでなく、高波数に対して非物理的な数値振動も少なく、安定で高精度な計算を行うことができるようになった。
【００７３】
なお、当業者であれば、第１実施形態に係る数値計算方法をコンピュータに実行させるためのプログラムを作成し、さらに、そのプログラムをＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disk）などの記録媒体に記録することができる。
【００７４】
＜第２実施形態＞
次に、第２実施形態として、前記したステップＳ２（図２参照）における補間関数の設定の処理について説明する。この第２実施形態では、第１実施形態と同様の点については重複説明を適宜省略し、第１実施形態との相違点について重点的に説明する。
【００７５】
図１０を参照して、補間関数の設定のイメージについて説明する。第１実施形態との違いは、第１の補間関数を用いて算出する値がｆtmpであり、そのｆtmpの座標がｘｊ＋α・Δｘ（０＜α＜１）であることである。
【００７６】
そして、第１実施形態と同様にして、第１の補間関数（式（４ａ））と第２の補間関数（式（６））を用いると、ｃ０、ｃ１、ｃ２、ｃ３の値は次のように求まる。
【数１２】
0005500637.tif 000013
【００７７】
ここで、式（８）、式（９）のときと同様に考えると、式（１３）におけるｃ１の元となる厳密解は、次の式（１４）のように表すことができる。
【数１３】
0005500637.tif 000014
【００７８】
また、式（１４）を変形して、次の式（１５）とする。
【数１４】
0005500637.tif 000015
【００７９】
ここで、式（１５）において、３階微分項は波の位相速度に、４階微分項は波の粘性散逸に影響を与える。そして、保存型ＩＤＯ法では波の位相速度が厳密解よりも遅いため、位相速度が上がる方に修正されるとよく、つまり、Ｅ３＜０を満たすのが好ましい。また、数値粘性は安定なものを導入するのがよく、つまり、Ｅ４＞０を満たすのが好ましい。したがって、それらを満たすように計算すると、αは次の式（１６）を満たせばよいことがわかる。
（√１７−３)／４＜α≦１／２  ・・・式（１６）
【００８０】
なお、式（１６）を満たすαのうち、式（１４）、式（１５）からわかるように、α＝１／２のとき、Ｅ３＝０となり、式（１５）におけるΔｘ２の項が消えるので、ｃ１は三次精度（誤差の大きさがΔｘ３に比例）となる。そして、α≠１／２のとき、ｃ１は二次精度（誤差の大きさがΔｘ２に比例）となる。したがって、これは、第１実施形態において、ｊ＋１／２点（ｊ点とｊ＋１の中間点。つまり、第２実施形態におけるα＝１／２のときの点）を仮想点とし、第１の補間関数を用いてその仮想点上のｆの値を算出した（補間した）ことの有効性を示したことにもなっている。また、ｃ２、ｃ３の精度についても同様である。
【００８１】
このように、第２実施形態によれば、数値計算装置１０による偏微分方程式の演算に使用する数値計算方法において、従来から風上補間のときに使用していたｆｊ、ｘｆｊ＋１／２、ｆｊ＋１の３つの値に加えて、第１の補間関数を用いて算出した仮想点ｊ＋α・Δｘにおける値ｆｊ＋α・Δｘを用いて第２の補間関数の係数を決定する。したがって、ｊ点とｊ＋１点の中間点以外に、式（１６）を満たす範囲でαを設定して決定した風下の点についての情報を用いて第２の補間関数の係数を決定し、その第２の補間関数を使用することで、高精度な計算を実現することができる。
【００８２】
＜第３実施形態＞
次に、第３実施形態として、前記したステップＳ２（図２参照）における補間関数の設定の処理について説明する。第３実施形態では、第１実施形態および第２実施形態の少なくとも一方と同様の点については重複説明を適宜省略し、第１実施形態および第２実施形態との相違点について重点的に説明する。
【００８３】
図１１を参照して、補間関数の設定のイメージについて説明する。第２実施形態との違いは、第１の補間関数の範囲がｊ−１点からｊ＋１点までであること、第１の補間関数が式（４ａ）の２次多項式ではなく、それを拡張した４次多項式（詳細な説明は省略）であること、および、それらにともなって、αの範囲が−１＜α＜０、０＜α＜１であること（つまり、仮想点はｊ−１点からｊ＋１点までの範囲内でｊ点以外）である。
【００８４】
第１実施形態と同様に、前記した４次多項式の第１の補間関数と、第２の補間関数（式（６））とを用いると、ｃ０、ｃ１、ｃ２、ｃ３の値は次の式（１７）のように求まる。
【数１５】
0005500637.tif 000016
【００８５】
ここで、式（８）、式（９）のときと同様に考えると、式（１７）におけるｃ１の元となる厳密解は、次の式（１８）のように表すことができる。
【数１６】
0005500637.tif 000017
【００８６】
また、式（１８）を変形して、次の式（１９）とする。
【数１７】
0005500637.tif 000018
【００８７】
ここで、式（１９）におけるＥ４について考えると、まず、分母に（２ａ＋１）の項が含まれるため、分母が０（α＝−１／２）に近づくと、Ｅ４の値が無限（無限大または無限小）に漸近するので計算が安定しない。また、数値粘性の面から考えると、Ｅ４の値はＥ４＞０を満たすことが望ましい。そして、−３／５＜α＜−１／２の範囲では、Ｅ４＞０を満たすが、Ｅ４の値が無限へと漸近することがあるため、計算の安定の面から好ましくない。したがって、０＜α＜１の範囲が、Ｅ４＞０を満たし、かつ、Ｅ４の値が無限へと漸近しないため、計算の安定のために好ましい。また、これは、ｃ２、ｃ３の場合も同様である。
【００８８】
このように、第３実施形態によれば、第１の補間関数の範囲をｊ−１点からｊ＋１点までに拡張し、第１の補間関数を４次多項式としたことで、αの値によらず計算を三次精度とすることができる。また、その際、αの範囲は、０＜α＜１であることが好ましく、これにより、本発明の特徴の１つである風下の適切な情報を用いることの有効性が示された。
【００８９】
なお、前記した３つの実施形態から一般性を導くと、（１）５個以上の物理量を用いて第１の補間関数を作る、あるいは、（２）第１実施形態で行ったようにｆｊ、ｘｆｊ＋１／２、ｆｊ＋１を用いて仮想点ｊ＋１／２点におけるｆｊ＋１／２を算出する、といういずれかの条件を満たすことにより、ｃ１を三次精度とすることができるといえる。
【００９０】
以上で本実施形態の説明を終えるが、本発明の態様はこれらに限定されるものではない。
例えば、本発明は、前記した実施形態に限定されず、偏微分方程式を使った流体の物理現象のシミュレーション全般に広く適用することができる。
【００９１】
また、使用する偏微分方程式における空間の次元数は、１次元以上の何次元でもよい。
また、本実施形態の数値計算方法を適用する流体に関する物理的な変数は、密度、速度、圧力、温度など、いかなる変数であってもよい。
【００９２】
また、ｊ−１点側を風上とし、ｊ＋１点側を風下とした場合について主に説明したが、逆に、ｊ−１点側を風下とし、ｊ＋１点側を風上としてもよい。つまり、連続した計算の中で、風上と風下が入れ替わった場合に、補間関数をそれぞれの場面で使い分けることができる。
また、本発明は、保存形ＩＤＯ法だけでなく、ＣＩＰ法等、他の計算手法にも適用することができる。
【００９３】
また、第１の補間関数の係数を決定する際、ｆｊ、ｘｆｊ＋１／２、ｆｊ＋１の３つの値を使うものとしたが、その他の値を使うものとしてもよい。
また、前記した、第１実施形態および第２実施形態における第１の補間関数、ならびに、第３実施形態における４次多項式の第１の補間関数は、一例であり、他の補間関数を用いてもよい。
【００９４】
また、第２の補間関数の係数としては、式（８）、（１０）〜（１２）に示した値そのものでなくても、それらに準ずる値を使用してもよい。それらに準ずる値とは、例えば、式（８）、（１０）〜（１２）における右辺の数字部分を少し変えた場合の値が考えられる。例えば、式（８）の右辺において、「５ｆｊ−１」を「４．９９９９ｆｊ−１」に置き換えても、計算の精度にあまり影響しないのは明らかである。また、例えば、ｊ点とｊ＋１の間に仮想点を考える際、ｊ点とｊ＋１のちょうど中間の点を仮想点としたが、それ以外に、中間の点から少しずれた点を仮想点としてもよい。その場合、第２の補間関数の係数は、式（８）、（１０）〜（１２）に示した値と少し異なったものとなるが、計算の精度にあまり影響がない範囲であれば、それらの少し異なった値も「準ずる値」に該当する。
【００９５】
その他、数値計算装置や各フローチャートなどの具体的な構成等について、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【符号の説明】
【００９６】
１      入力部
２      記憶部
３      メモリ
４      演算部
５      出力部
１０    数値計算装置
【図１】
0005500637.tif 000025
【図２】
0005500637.tif 000026
【図３】
0005500637.tif 000027
【図４】
0005500637.tif 000028
【図５】
0005500637.tif 000029
【図６】
0005500637.tif 000030
【図７】
0005500637.tif 000031
【図１０】
0005500637.tif 000032
【図１１】
0005500637.tif 000033
【図１２】
0005500637.tif 000034
【図８】
0005500637.tif 000035
【図９】
0005500637.tif 000036