(57)【要約】【課題】  精度が向上した移流方程式の数値解法を提供すること。【解決手段】  空間上で定義された物理量及び速度場とから、この物理量の移流を求める移流方程式を数値計算によって解く。初期値記憶手段１２は、初期時刻ｔ０における物理量ａ（ｘ、ｔ０）を記憶している。流跡線算出手段は、前記空間上の速度場ｕから、前記空間の速度場を表す流跡線を表すＧ（ｘ、ｔ）を求める。流跡線記憶手段は、流跡線ＧにおけるＧ（ｘ、ｔ）を記憶する。移流計算手段は、初期時刻ｔ０における物理量である前記ａ（ｘ、ｔ０）に前記流跡線を表すＧ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める。このようにして、移流方程式を解くことができる。【選択図】      図２
2006331351.tif 000002
【特許請求の範囲】
【請求項１】
所定の空間と、その空間上で定義された物理量と、前記空間上で定義された速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く方法において、
初期時刻ｔ０における物理量ａ（ｘ、ｔ０）（ここで、ｘは空間上の位置ベクトル）を初期値記憶手段に記憶する初期値記憶ステップと、
流跡線算出手段が、前記空間上の速度場ｕから、前記空間の速度場を表す流跡線を表すＧ（ｘ、ｔ）（ここで、ｔは時間を表す）を求める流跡線算出ステップと、
流跡線記憶手段が、前記空間上の速度場における流跡線を表すＧ（ｘ、ｔ）を記憶するステップと、
移流計算手段が、初期時刻ｔ０における物理量である前記ａ（ｘ、ｔ０）に前記流跡線を表すＧ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める移流計算ステップと、
を含むことを特徴とする移流方程式を解く方法。
【請求項２】
請求項１記載の移流方程式を解く方法であって、
前記流跡線算出ステップにおいて、
前記流跡線算出手段は、前記速度場に沿って、前記空間上の粒子を追跡することによって流跡線を表す前記Ｇ（ｘ、ｔ）求めることを特徴とする移流方程式を解く方法。
【請求項３】
請求項１記載の移流方程式を解く方法であって、
前記流跡線算出ステップにおいて、
前記流跡線算出手段は、前記速度場を表すベクトル場を時間微分によって離散化し、得られた差分方程式によって、流跡線を表す前記Ｇ（ｘ、ｔ）を算出することを特徴とする移流方程式を解く方法。
【請求項４】
請求項１記載の移流方程式を解く方法であって、
前記流跡線算出ステップにおいて、
前記流跡線算出手段は、以下の式によって、流跡線を表すＧ（ｘ、ｔ）を算出することを特徴とする移流方程式を解く方法。
【数１】
2006331351.tif 000090
ここで、Δｔは、所定の時間ステップを表し、ベクトルｕは速度場を表す。
【請求項５】
所定の空間と、その空間上で定義された物理量と、前記空間上で定義された速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く方法において、
初期値記憶手段が、初期時刻ｔ０における物理量ａ（ｘ、ｔ０）（ここで、ｘは空間上の位置ベクトル）を記憶する初期値記憶ステップと、
写像関数算出手段が、初期時刻ｔ０における移流前の空間Ω０から、時刻ｔｎにおける移流後の空間Ωｎへの写像関数Ｆ（ｘ、ｔ）を求める写像関数算出ステップと、
逆写像関数算出手段が、時刻ｔｎにおける移流後の空間Ωｎから、時刻ｔ０における移流前の空間Ω０への逆写像関数Ｇ（ｘ、ｔ）を求める逆写像算出ステップと、
写像関数・逆写像関数記憶手段が、前記写像関数Ｆ（ｘ、ｔ）及び逆写像関数Ｇ（ｘ、ｔ）を記憶する写像関数・逆写像関数記憶ステップと、
非移流項加算手段が、前記写像関数Ｆ（ｘ、ｔ）を用いて、移流後の時刻ｔｎにおける空間Ωｎで計算された非移流項ｄ（ｘ、ｔｎ）を、時刻ｔ０における物理量ａ（ｘ、ｔ０）に加える非移流項加算ステップと、
移流計算手段が、前記加算後の物理量ａ（ｘ、ｔ０）に、前記逆写像関数Ｇ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める移流計算ステップと、
を含むことを特徴とする移流方程式を解く方法。
【請求項６】
請求項５記載の移流方程式を解く方法であって、
前記写像関数算出ステップにおいて、
前記写像関数算出手段が、前記速度場に沿って、前記空間上の粒子を追跡することによって写像関数Ｆ（ｘ、ｔ）を求め、
前記逆写像関数算出ステップにおいて、
前記逆写像関数算出手段が、前記速度場に沿って、前記空間上の粒子を追跡することによって逆写像関数Ｇ（ｘ、ｔ）を求めることを特徴とする移流方程式を解く方法。
【請求項７】
請求項５記載の移流方程式を解く方法であって、
前記写像関数算出ステップにおいて、
前記写像関数算出手段が、前記速度場を表すベクトル場を時間微分によって離散化し、得られた差分方程式によって、前記写像関数Ｆ（ｘ、ｔ）を算出し、
前記逆写像算出ステップにおいて、
前記逆写像関数算出手段が、前記速度場を表すベクトル場を時間微分によって離散化し、得られた差分方程式によって、前記逆写像関数Ｇ（ｘ、ｔ）を算出することを特徴とする移流方程式を解く方法。
【請求項８】
請求項５記載の移流方程式を解く方法であって、
前記写像関数算出ステップにおいて、
前記写像関数算出手段が、以下の式によって、前記写像関数Ｆ（ｘ、ｔ）を算出し、
【数２】
2006331351.tif 000091
前記逆写像算出ステップにおいて、
前記逆写像関数算出手段が、以下の式によって、前記逆写像関数Ｇ（ｘ、ｔ）を算出することを特徴とする移流方程式を解く方法。
【数３】
2006331351.tif 000092
ここで、Δｔは、所定の時間ステップを表し、ベクトルｕは速度場を表す。
【請求項９】
コンピュータを、所定の空間と、その空間上で定義された物理量と、前記空間上で定義された速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く装置として動作させるプログラムにおいて、前記コンピュータに、
初期時刻ｔ０における物理量ａ（ｘ、ｔ０）（ここで、ｘは空間上の位置ベクトル）を記憶する初期値記憶手順と、
前記空間上の速度場ｕから、前記空間の速度場を表す流跡線を表すＧ（ｘ、ｔ）（ここで、ｔは時間を表す）を求める流跡線算出手順と、
前記空間上の速度場における流跡線を表すＧ（ｘ、ｔ）を記憶する流跡線記憶手段と、
初期時刻ｔ０における物理量である前記ａ（ｘ、ｔ０）に前記流跡線を表すＧ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める移流計算手順と、
を実行させることを特徴とするプログラム。
【請求項１０】
請求項９記載のプログラムであって、
前記流跡線算出手順においては、前記コンピュータに、前記速度場に沿って、前記空間上の粒子を追跡することによって流跡線を表す前記Ｇ（ｘ、ｔ）を求める動作を実行させることを特徴とするプログラム。
【請求項１１】
請求項９記載のプログラムであって、
前記流跡線算出手順においては、前記コンピュータに、前記速度場を表すベクトル場を時間微分で離散化し、得られた差分方程式によって、流跡線を表す前記Ｇ（ｘ、ｔ）を算出する動作を実行させることを特徴とする移流方程式を解く方法。
【請求項１２】
請求項９記載のプログラムであって、
前記流跡線算出手順において、前記コンピュータに、以下の式によって、流跡線を表すＧ（ｘ、ｔ）を算出する動作を実行させることを特徴とする移流方程式を解く方法。
【数４】
2006331351.tif 000093
ここで、Δｔは、所定の時間ステップを表し、ベクトルｕは速度場を表す。
【請求項１３】
コンピュータを、所定の空間と、その空間上で定義された物理量と、前記空間上で定義された速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く装置として動作させるプログラムにおいて、
初期時刻ｔ０における物理量ａ（ｘ、ｔ０）（ここで、ｘは空間上の位置ベクトル）を記憶する初期値記憶手順と、
初期時刻ｔ０における移流前の空間Ω０から、時刻ｔｎにおける移流後の空間Ωｎへの写像関数Ｆ（ｘ、ｔ）を求める写像関数算出手順と、
時刻ｔｎにおける移流後の空間Ωｎから、時刻ｔ０における移流前の空間Ω０への逆写像関数Ｇ（ｘ、ｔ）を求める逆写像算出手順と、
前記写像関数Ｆ（ｘ、ｔ）及び逆写像関数Ｇ（ｘ、ｔ）を記憶する写像関数・逆写像関数記憶手順と、
前記写像関数Ｆ（ｘ、ｔ）を用いて、移流後の時刻ｔｎにおける空間Ωｎで計算された非移流項ｄ（ｘ、ｔｎ）を、時刻ｔ０における物理量ａ（ｘ、ｔ０）に加える非移流項加算手順と、
前記加算後の物理量ａ（ｘ、ｔ０）に、前記逆写像関数Ｇ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める移流計算手順と、
を実行させることを特徴とするプログラム。
【請求項１４】
請求項１３記載のプログラムであって、
前記写像関数算出手順において、前記コンピュータに、前記速度場に沿って、前記空間上の粒子を追跡することによって写像関数Ｆ（ｘ、ｔ）を求める動作を実行させ、
前記逆写像算出ステップにおいて、前記コンピュータに、前記速度場に沿って、前記空間上の粒子を追跡することによって逆写像関数Ｇ（ｘ、ｔ）を求める動作を実行させることを特徴とするプログラム。
【請求項１５】
請求項１３記載のプログラムであって、
前記写像関数算出手順において、前記コンピュータに、前記速度場を表すベクトル場を時間微分で離散化し、得られた差分方程式によって、前記写像関数Ｆ（ｘ、ｔ）を算出する動作を実行させ、
前記逆写像関数算出手順において、前記コンピュータに、前記速度場を表すベクトル場を時間微分によって離散化し、得られた差分方程式によって、前記逆写像関数Ｇ（ｘ、ｔ）を算出する動作を実行させることを特徴とするプログラム。
【請求項１６】
請求項１３記載のプログラムであって、
前記写像関数算出手順において、前記コンピュータに、以下の式によって、前記写像関数Ｆ（ｘ、ｔ）を算出する動作を実行させ、
【数５】
2006331351.tif 000094
前記逆写像関数算出手順において、前記コンピュータに、以下の式によって、前記逆写像関数Ｇ（ｘ、ｔ）を算出する動作を実行させることを特徴とするプログラム。
【数６】
2006331351.tif 000095
ここで、Δｔは、所定の時間ステップを表し、ベクトルｕは速度場を表す。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物理現象を解析するための方程式である移流方程式を解く装置及び方法に関する。また、コンピュータを移流方程式を解く装置として動作させるプログラムに関する。
【背景技術】
【０００２】
移流方程式とは、ある物理量が与えられた速度場によって伝搬する現象を記述した方程式である。移流は、流体の運動を初めとして、多くの物理現象に関わるため、移流方程式を高精度に解くことは、物理現象の数値計算において重要な課題となっている。
【０００３】
移流方程式
自然界には、波に関する現象が数多く存在する。海面を伝わる波や、空気中を伝わる波、種々の電磁波など、多くの波が存在する。この波を記述する方程式は移流方程式と呼ばれる。このように、移流方程式は波が伝播する方程式であり、種々の物理現象を解析するのに重要な役割を果たす方程式である。
【０００４】
したがって、この移流方程式を効率よく解くことは、単に学術的な面に留まらず、産業上重要な意味を有する。
【０００５】
さて、方程式をいわゆる数学的な手法を用いて代数的に厳密解を求めるのは一般に困難である場合が多い。そのため、従来から知られている移流方程式を「解く」手法は、いわゆる数値解法（数値計算）の手法を用いて、近似的な解を求めようとするアプローチがほとんどである。本発明も、このような数値解法の一つである。
【０００６】
従来の解法とＣＩＰ法
図８には、この「移流」の概念を表すグラフが示されている。この図においては、局所的な関数（物質）が時空間を所定の速度ｕで移動していく様子が示されている。このグラフでは、横軸が空間上のｘ座標を表し、奥行きが時間ｔを表す。縦軸は、ｆ（ｘ−ｕｔ）であり、所定の物理量の値である。図６のような局所的に存在する「物質」の場合は、密度のような物理量が好ましい。図８では、物質の例を示したが、波動でも同様である。
【０００７】
移流方程式の数値解法は、様々な方法が提案されている。例えば、下記非特許文献１には、ＣＩＰ（Constrained Interpolation Profile）法と呼ばれる数値解法が記載されている。
【０００８】
このＣＩＰ法は、端的に言えば、物理量を移流させる際に、各点の物理量を単に移動させただけでなく、各点間の傾きまで移流の対象とした点を特徴とする方法である。
【０００９】
このＣＩＰ法の原理を簡単に表した図が図７（下記非特許文献１の２６ページから引用）である。この図９の（１）に実線で示されているようなプロファイルを移動する際に、各点（●で示されている）の値だけでなく、その傾き（矢印）も移流されている。その結果、図９（２）に示すように、元の形に近い形で、移流後の結果が得られている。
【００１０】
その他の従来技術
方程式そのものを解く手法ではないが、明確な境界を持つような物理量の移流に関して、境界の補足に特化した手法が知られている。例えば、このような手法が下記非特許文献２に記載されている。このような手法によれば、計算誤差を抑えることができると非特許文献２に記載されている。
【００１１】
下記非特許文献３には、 テクスチャマッピングによる流れ場の可視化に関する技術が記載されている。この技術は、後述するように、本件発明の処理の制御の判断基準として用いることができる可能性がある。この事項は後に再び詳述する。
【００１２】
また、下記特許文献１には、雨域のレーダー画像に基づき、雨域の輪郭線から速度ベクトルを算出し、この速度場を雨域の濃淡情報と共に移流方程式に代入する技術が開示されている。この移流方程式を解き、雨域の変化予測をする発明が開示されている。
【００１３】
また、下記特許文献２には、移流方程式を利用した気液自由界面シミュレーション装置が開示されている。ここに記載されている発明によれば、数値計算で求めた各格子の流体率や密度が数値拡散によって誤差を生じた場合に、正しい値との差を気液自由界面法線方向に移動させることによって、数値拡散を補正する技術が開示されている。
【００１４】
また、下記特許文献３には、入力画像から、連続した動画を生成する方法が開示されている。ここに記載の発明によれば、入力画像である静止画像に対して速度場を設定し、その速度場に基づいて移流方程式を適用し、出力画像を生成する技術が開示されている。
【００１５】
【非特許文献１】矢部孝、内海隆行、尾形陽一、「ＣＩＰ法−原子から宇宙までを解くマルチスケール解法」、森北出版、２００３年
【非特許文献２】Stanley J. Osher, Ronald P. Fedkiw :"Level Set Methods and Dynamic Implicit Surfaces", Springer, 2002
【非特許文献３】Fablic Neyret :Advected Textures, Proceedings of the 2003 ACM SIGRAPH/Eurographics Symposium on Computer animation, pp147-153, 2003
【特許文献１】特開平１０−２６０２６７号公報（特許登録：３３７７０７６号）
【特許文献２】特開２００４−１１８３９４号公報
【特許文献３】特開２００４−３３４６９４号公報
【発明の開示】
【発明が解決しようとする課題】
【００１６】
さて、移流方程式を離散化して計算（数値計算）しようとする場合、２つの記述法（見方・立場）が存在する。
【００１７】
その一つは、ラグランジアン（Lagrangian）（ラグランジュ式記述）と呼ばれる立場・記法であり、流れと共に視点を移動させていく見方である。他の一つは、オイラリアン（Eulaerian）（オイラー式記述）と呼ばれる立場・記法であり、空間に視点を固定し、その視点に対して移流を計算するものである。
【００１８】
さて、オイラー的な空間に固定の格子を用いて移流方程式を解いた場合、計算手法に起因する数値拡散や振動などが問題となる場合が多い。
【００１９】
明確な境界を持つような物理量の移流に関しては、上記非特許文献２に記載されているような手法によって、計算誤差を抑えることができる場合も存在するが、特に、明確な境界を持たないような物理量の移流については、上記の手法が使えず計算誤差を抑えることができない。
【００２０】
したがって、移流方程式の解法以外の補助的な手法ではなく、高精度な移流方程式の解法そのものが必要とされている。
【００２１】
本発明は、移流方程式に関するこのような状況に鑑みなされたものであり、従来より精度が向上した移流方程式の数値解法を提供することを目的とする。
【課題を解決するための手段】
【００２２】
Ａ．移流方程式を解く方法
以下、上記移流方程式を解く方法の発明である。
【００２３】
（１）本発明は、上記課題を解決するために、所定の空間と、その空間上で定義された物理量と、前記空間上で定義された速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く方法において、初期時刻ｔ０における物理量ａ（ｘ、ｔ０）（ここで、ｘは空間上の位置ベクトル）を初期値記憶手段に記憶する初期値記憶ステップと、流跡線算出手段が、前記空間上の速度場ｕから、前記空間の速度場を表す流跡線を表すＧ（ｘ、ｔ）（ここで、ｔは時間を表す）を求める流跡線算出ステップと、流跡線記憶手段が、前記空間上の速度場における流跡線を表すＧ（ｘ、ｔ）を記憶するステップと、移流計算手段が、初期時刻ｔ０における物理量である前記ａ（ｘ、ｔ０）に前記流跡線を表すＧ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める移流計算ステップと、を含むことを特徴とする移流方程式を解く方法である。
【００２４】
（２）また、本発明は、上記（１）の移流方程式を解く方法において、前記流跡線算出ステップにおいて、前記流跡線算出手段は、前記速度場に沿って、前記空間上の粒子を追跡することによって流跡線を表す前記Ｇ（ｘ、ｔ）求めることを特徴とする移流方程式を解く方法である。
【００２５】
このような構成によって、従来技術の固定格子におけるＧの数値拡散の問題は生じなくなる。それは、本件発明が、Full-LagrangianでＧを求めているからである。また、本発明で、ａ（ｘ、ｔ）において数値拡散が存在しないことの本質は、時間発展による誤差が蓄積しないことである。
【００２６】
（３）また、本発明は、上記（１）記載の移流方程式を解く方法において、前記流跡線算出ステップにおいて、前記流跡線算出手段は、前記速度場を表すベクトル場を時間微分で離散化し、得られた差分方程式によって、流跡線を表す前記Ｇ（ｘ、ｔ）を算出することを特徴とする移流方程式を解く方法である。
【００２７】
このような構成によって、従来技術の固定格子におけるＧの数値拡散の問題は生じなくなる。それは、本件発明が、Full-LagrangianでＧを求めているからである。また、本発明で、ａ（ｘ、ｔ）において数値拡散が存在しないことの本質は、時間発展による誤差が蓄積しないことである。
【００２８】
（４）また、本発明は、上記（１）記載の移流方程式を解く方法において、前記流跡線算出ステップにおいて、前記流跡線算出手段は、以下の式によって、流跡線を表すＧ（ｘ、ｔ）を算出することを特徴とする移流方程式を解く方法である。
【００２９】
【数１】
2006331351.tif 000003
ここで、Δｔは、所定の時間ステップを表し、ベクトルｕは速度場を表す。
【００３０】
このような構成によって、従来技術の固定格子におけるＧの数値拡散の問題は生じなくなる。それは、本件発明が、Full-LagrangianでＧを求めているからである。また、本発明で、ａ（ｘ、ｔ）において数値拡散が存在しないことの本質は、時間発展による誤差が蓄積しないことである。
【００３１】
（５）また、本発明は、所定の空間と、その空間上で定義された物理量と、前記空間上で定義された速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く方法において、初期値記憶手段が、初期時刻ｔ０における物理量ａ（ｘ、ｔ０）（ここで、ｘは空間上の位置ベクトル）を記憶する初期値記憶ステップと、写像関数算出手段が、初期時刻ｔ０における移流前の空間Ω０から、時刻ｔｎにおける移流後の空間Ωｎへの写像関数Ｆ（ｘ、ｔ）を求める写像関数算出ステップと、逆写像関数算出手段が、時刻ｔｎにおける移流後の空間Ωｎから、時刻ｔ０における移流前の空間Ω０への逆写像関数Ｇ（ｘ、ｔ）を求める逆写像算出ステップと、写像関数・逆写像関数記憶手段が、前記写像関数Ｆ（ｘ、ｔ）及び逆写像関数Ｇ（ｘ、ｔ）を記憶する写像関数・逆写像関数記憶ステップと、非移流項加算手段が、前記写像関数Ｆ（ｘ、ｔ）を用いて、移流後の時刻ｔｎにおける空間Ωｎで計算された非移流項ｄ（ｘ、ｔｎ）を、時刻ｔ０における物理量ａ（ｘ、ｔ０）に加える非移流項加算ステップと、移流計算手段が、前記加算後の物理量ａ（ｘ、ｔ０）に、前記逆写像関数Ｇ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める移流計算ステップと、を含むことを特徴とする移流方程式を解く方法である。
【００３２】
なお、時刻ｔｎとは、その時刻における移流、すなわち物理量を求めたい時刻である。
【００３３】
（６）また、本発明は、上記（５）記載の移流方程式を解く方法であって、前記写像関数算出ステップにおいて、前記写像関数算出手段が、前記速度場に沿って、前記空間上の粒子を追跡することによって写像関数Ｆ（ｘ、ｔ）を求め、前記逆写像関数算出ステップにおいて、前記逆写像関数算出手段が、前記速度場に沿って、前記空間上の粒子を追跡することによって逆写像関数Ｇ（ｘ、ｔ）を求めることを特徴とする移流方程式を解く方法である。
【００３４】
このような構成によって、従来技術の固定格子におけるＧの数値拡散の問題は生じなくなる。それは、本件発明が、Full-LagrangianでＦ及びＧを求めているからである。また、本発明で、ａ（ｘ、ｔ）において数値拡散が存在しないことの本質は、時間発展による誤差が蓄積しないことである。
【００３５】
（７）また、本発明は、上記（５）記載の移流方程式を解く方法であって、前記写像関数算出ステップにおいて、前記写像関数算出手段が、前記速度場を表すベクトル場を時間微分で離散化し、得られた差分方程式によって、前記写像関数Ｆ（ｘ、ｔ）を算出し、前記逆写像算出ステップにおいて、前記逆写像関数算出手段が、前記速度場を表すベクトル場を時間微分で離散化し、得られた差分方程式によって、前記逆写像関数Ｇ（ｘ、ｔ）を算出することを特徴とする移流方程式を解く方法である。
【００３６】
このような構成によって、従来技術の固定格子におけるＧの数値拡散の問題は生じなくなる。それは、本件発明が、Full-LagrangianでＦ及びＧを求めているからである。また、本発明で、ａ（ｘ、ｔ）において数値拡散が存在しないことの本質は、時間発展による誤差が蓄積しないことである。
【００３７】
（８）また、本発明は、上記（５）記載の移流方程式を解く方法であって、前記写像関数算出ステップにおいて、前記写像関数算出手段が、以下の式によって、前記写像関数Ｆ（ｘ、ｔ）を算出し、
【数２】
2006331351.tif 000004
前記逆写像算出ステップにおいて、前記逆写像関数算出手段が、以下の式によって、前記逆写像関数Ｇ（ｘ、ｔ）を算出することを特徴とする移流方程式を解く方法である。
【００３８】
【数３】
2006331351.tif 000005
ここで、Δｔは、所定の時間ステップを表し、ベクトルｕは速度場を表す。
【００３９】
このような構成によって、従来技術の固定格子におけるＧの数値拡散の問題は生じなくなる。それは、本件発明が、Full-LagrangianでＦ及びＧを求めているからである。また、本発明で、ａ（ｘ、ｔ）において数値拡散が存在しないことの本質は、時間発展による誤差が蓄積しないことである。
【００４０】
Ｂ．コンピュータを、移流方程式を解く装置として動作させるプログラム
以下、コンピュータを、上記移流方程式を解く装置として動作させるプログラムの発明であり、その作用効果は、実質的に上記方法の発明と同様である。
【００４１】
（９）本発明は、上記課題を解決するために、コンピュータを、所定の空間と、その空間上で定義された物理量と、前記空間上で定義された速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く装置として動作させるプログラムにおいて、前記コンピュータに、初期時刻ｔ０における物理量ａ（ｘ、ｔ０）（ここで、ｘは空間上の位置ベクトル）を記憶する初期値記憶手順と、前記空間上の速度場ｕから、前記空間の速度場を表す流跡線を表すＧ（ｘ、ｔ）（ここで、ｔは時間を表す）を求める流跡線算出手順と、前記空間上の速度場における流跡線を表すＧ（ｘ、ｔ）を記憶する流跡線記憶手段と、初期時刻ｔ０における物理量である前記ａ（ｘ、ｔ０）に前記流跡線を表すＧ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める移流計算手順と、を実行させることを特徴とするプログラムである。
【００４２】
（１０）また、本発明は、上記（９）記載のプログラムであって、前記流跡線算出手順においては、前記コンピュータに、前記速度場に沿って、前記空間上の粒子を追跡することによって流跡線を表す前記Ｇ（ｘ、ｔ）を求める動作を実行させることを特徴とするプログラム。
【００４３】
（１１）また、本発明は、上記（９）記載のプログラムであって、前記流跡線算出手順においては、前記コンピュータに、前記速度場を表すベクトル場を時間微分で離散化し、得られた差分方程式によって、流跡線を表す前記Ｇ（ｘ、ｔ）を算出する動作を実行させることを特徴とする移流方程式を解く方法である。
【００４４】
（１２）また、本発明は、上記（９）記載のプログラムであって、前記流跡線算出手順において、前記コンピュータに、以下の式によって、流跡線を表すＧ（ｘ、ｔ）を算出する動作を実行させることを特徴とする移流方程式を解く方法である。
【００４５】
【数４】
2006331351.tif 000006
ここで、Δｔは、所定の時間ステップを表し、ベクトルｕは速度場を表す。
【００４６】
（１３）また、本発明は、コンピュータを、所定の空間と、その空間上で定義された物理量と、前記空間上で定義された速度場と、から前記物理量の移流を求める移流方程式を数値計算によって解く装置として動作させるプログラムにおいて、初期時刻ｔ０における物理量ａ（ｘ、ｔ０）（ここで、ｘは空間上の位置ベクトル）を記憶する初期値記憶手順と、初期時刻ｔ０における移流前の空間Ω０から、時刻ｔｎにおける移流後の空間Ωｎへの写像関数Ｆ（ｘ、ｔ）を求める写像関数算出手順と、時刻ｔｎにおける移流後の空間Ωｎから、時刻ｔ０における移流前の空間Ω０への逆写像関数Ｇ（ｘ、ｔ）を求める逆写像算出手順と、前記写像関数Ｆ（ｘ、ｔ）及び逆写像関数Ｇ（ｘ、ｔ）を記憶する写像関数・逆写像関数記憶手順と、前記写像関数Ｆ（ｘ、ｔ）を用いて、移流後の時刻ｔｎにおける空間Ωｎで計算された非移流項ｄ（ｘ、ｔｎ）を、時刻ｔ０における物理量ａ（ｘ、ｔ０）に加える非移流項加算手順と、前記加算後の物理量ａ（ｘ、ｔ０）に、前記逆写像関数Ｇ（ｘ、ｔ）を適用し、ａ（ｘ、ｔｎ）＝ａ（Ｇ（ｘ、ｔｎ）、０）を算出することによって、時刻ｔｎにおける物理量ａ（ｘ、ｔｎ）を、求める移流計算手順と、を実行させることを特徴とするプログラムである。
【００４７】
（１４）また、本発明は、上記（１３）記載のプログラムであって、前記写像関数算出手順において、前記コンピュータに、前記速度場に沿って、前記空間上の粒子を追跡することによって写像関数Ｆ（ｘ、ｔ）を求める動作を実行させ、前記逆写像算出ステップにおいて、前記コンピュータに、前記速度場に沿って、前記空間上の粒子を追跡することによって逆写像関数Ｇ（ｘ、ｔ）を求める動作を実行させることを特徴とするプログラムである。
【００４８】
（１５）また、本発明は、上記（１３）記載のプログラムであって、前記写像関数算出手順において、前記コンピュータに、前記速度場を表すベクトル場を時間微分で離散化し、得られた差分方程式によって、前記写像関数Ｆ（ｘ、ｔ）を算出する動作を実行させ、前記逆写像関数算出手順において、前記コンピュータに、前記速度場を表すベクトル場を時間微分によって離散化し、得られた差分方程式によって、前記逆写像関数Ｇ（ｘ、ｔ）を算出する動作を実行させることを特徴とするプログラムである。
【００４９】
（１６）また、本発明は、上記（１３）記載のプログラムにおいて、前記写像関数算出手順において、前記コンピュータに、以下の式によって、前記写像関数Ｆ（ｘ、ｔ）を算出する動作を実行させ、
【数５】
2006331351.tif 000007
前記逆写像関数算出手順において、前記コンピュータに、以下の式によって、前記逆写像関数Ｇ（ｘ、ｔ）を算出する動作を実行させることを特徴とするプログラムである。
【００５０】
【数６】
2006331351.tif 000008
ここで、Δｔは、所定の時間ステップを表し、ベクトルｕは速度場を表す。
【発明の効果】
【００５１】
以上述べたように、本発明によれば、精度や計算速度が向上した移流方程式の解法及び装置が得られる。
【発明を実施するための最良の形態】
【００５２】
以下、本発明の好適な実施の形態を図面に基づき説明する。以下、順次、第１章定式化及び実装を述べ、第２章計算例、第３章本実施の形態の効果とまとめを述べる。
【００５３】
第１章  定式化及び実装
（１−１）流跡線による移流方程式の定式化
本実施の形態で提案する手法では、非移流項を含む一般化された移流方程式を扱うことが可能であるが、説明の便宜のため、下記（１）式のような非移流項を含まない線形移流方程式について説明する。
【００５４】
【数７】
2006331351.tif 000009
ここで、
【数８】
2006331351.tif 000010
は、速度場を表し、
【数９】
2006331351.tif 000011
は、移流の対象であるスカラー値を表す。
【００５５】
このような移流方程式は、Semi-Lagrangian法（非特許文献１及び非特許文献２参照）によれば、下記（２）式のような形式で表すことができる。
【００５６】
【数１０】
2006331351.tif 000012
ここで、Δｔは、ある微小時間（時間刻み幅）を表し、ｔｎ＝ｎΔｔである。また、
【数１１】
2006331351.tif 000013
は、
【数１２】
2006331351.tif 000014
となる流跡線を表している。
【００５７】
なお、本特許では左向き矢印は、右辺から左辺への代入を表している。Semi-Lagrangian法では、上記（２）式によって、
【数１３】
2006331351.tif 000015
の時間発展を求めており、
【数１４】
2006331351.tif 000016
については、一時的に表れるだけである。
【００５８】
これに対して、本発明では、この一時的に表れる「数１４」についての情報を陽に保持する。まず、上記（２）式を各計算ステップに連続して適用すれば、
【数１５】
2006331351.tif 000017
となる。したがって、この（３）式で
【数１６】
2006331351.tif 000018
を整理して、両辺を時間軸について−Δｔだけずらせば、
【数１７】
2006331351.tif 000019
となることが理解されよう（（４）式）。本発明では、上記（４）式のような定式化において、以下のような関数
【数１８】
2006331351.tif 000020
を導入する。これによって、以下の（５）式が得られる。
【００５９】
【数１９】
2006331351.tif 000021
ここで
【数２０】
2006331351.tif 000022
とすれば、
【数２１】
2006331351.tif 000023
は、ｔ＝０で
【数２２】
2006331351.tif 000024
に存在する粒子の流跡線であることに注意を向けられたい。このように本実施の形態では、初期値の
【数２３】
2006331351.tif 000025
と、
【数２４】
2006331351.tif 000026
とを保持することによって、各計算ステップで必要な
【数２５】
2006331351.tif 000027
を、
【数２６】
2006331351.tif 000028
を用いて
【数２７】
2006331351.tif 000029
から毎回構築している。
【００６０】
したがって、本実施の形態の手法は、時間発展による誤差の蓄積がないという特徴を有する。
【００６１】
【数２８】
2006331351.tif 000030
の更新については、既に述べたように、この関数は流跡線となるので、粒子を速度場によって追跡すればよい。したがって、時間１次精度のオイラー法によって離散化すれば、
【数２９】
2006331351.tif 000031
となる（（６）式）。これは移流に伴う数値拡散などの誤差を含まないFull-Lagrangianに基づいた更新であるので、上記流跡線の関数（「数２１」等）の更新においても、従来の固定格子における数値拡散などの問題を回避することができる。
【００６２】
なお、ここでは、時間１次精度のオイラー法によって離散化する手法を述べたが、時間微分を離散化する手法ならどのようなものでも利用可能である。例えば、ルンゲクッタ法、リープフロッグ法などが挙げられる。以下、本文で、時間１次精度のオイラー法と述べている部分においては、すべて、時間微分を離散化する手法ならどのようなものでも利用可能である。
【００６３】
既に述べたように、Semi-Lagrangian法では、このような流跡線の情報を持たずに、移流の対象となる値を直接、次のステップでの値に書き換えている。数学的にはこの方法でも（理論上は）誤差は生じないはずであるが、「数値計算上」は、数値拡散などの誤差が生じることが知られている。これは、移流の対象となる値が離散的にしか記述されていないため、更新に必要な上流での値は何らかの補間によって求めざるを得ないためである。補間方法には、様々な方法が考えられるが、「補間」であるから少なからず誤差は存在する。そのため、時間発展による誤差の蓄積が生じ、数値拡散や振動の原因となっている。
【００６４】
本実施の形態の手法と、従来のSemi-Lagrangian法の相違を、概念的に表した説明図が図１に示されている。
【００６５】
この図１（１）に示すように、従来のSemi-Lagrangian法においては、各時点での物理量
【数３０】
2006331351.tif 000032
は、
【数３１】
2006331351.tif 000033
によって、逐次求められていく。ここで、
【数３２】
2006331351.tif 000034
は、離散的に定義されているので、位置
【数３３】
2006331351.tif 000035
における物理量
【数３４】
2006331351.tif 000036
を得るために補間を用いる必要がある。この補間による誤差は、図１（１）に示すように更新（補間）を行うごとに累積的に蓄積して行くので、誤差が拡大し、その結果、数値拡散が生じてしまうのである。
【００６６】
これに対して、本実施の形態の手法は、
【数３５】
2006331351.tif 000037
を逐次、更新していくのではなく、この値を、初期値
【数３６】
2006331351.tif 000038
から関数
【数３７】
2006331351.tif 000039
を用いて求めている。この初期値は、移流の計算においては常に一定である。
【００６７】
本実施の形態の手法は、直感的に理解しやすいものであると考える。それは各時点における物理量は、初期値及び速度場によって全て完全に記述されるはずだからである。したがって、初期値に基づき各時点の値を算出することは理にかなったことである。
【００６８】
図１（２）には、本実施の形態の手法の原理を示す概念図が示されている。この図に示すように、移流の計算においては、常に初期値が用いられ、各時点の値を算出する際に必要な上記関数が逐次更新されていく。この図１（２）にも示すように、各時点での値を求める際には、本実施の形態の手法でも「補間」は必要であり、そこで誤差が生じるかもしれない。しかし、図１（２）に示すように、この誤差は各時点の計算における誤差であり、誤差が累積していくことはないのである。本実施の形態の手法では、常に初期値に基づき計算が行われており、この初期値には誤差が累積していないからである。
【００６９】
なお、様々な補間に関する説明が上記非特許文献１に記載されている。
【００７０】
本実施の形態においては、初期値を保持して後の計算でも用いるという特徴を有している。このような特徴を有しているので、本実施の形態の手法は、構造解析で用いられるTotal-Lagrangianに基づいた手法と位置づけることができよう。これに対して、移流方程式で従来提案されてきたSemi-Lagrangian法は、直前の値を用いて計算を進めるという特徴を有している。このような特徴は、このSemi-Lagrangian法がUpdate-Lagrangianに基づく手法であることを意味している。
【００７１】
また、本実施の形態の手法では、空間に固定された計算点を用いているため、純粋な意味でのFull-Lagrangianではないが、数値拡散を起こさない点などで、Full-Lagrangianにおける利点を備えたEulerian的な手法であると言えよう。
【００７２】
（１−２）実装
上記（１−１）で述べた原理に基づく移流方程式を解く装置の構成ブロック図が図２に示されている。
【００７３】
この図２に示すように、移流方程式を解く装置１０は、初期値記憶手段１２と、流跡線算出手段１４と、流跡線記憶手段１６と、移流計算手段１８と、から構成されている。
【００７４】
初期値記憶手段１２は、外部から供給されてきた初期値を記憶する。この初期値は、既に説明したように、種々の物理量を採用することができる。水溶液中の物質の濃度でもよいし、音波の音圧でもよい。また、固体の物質でもよい。この初期値は、上で説明した例によれば、ａ（ｘ、０）等の表記法を採用することが好ましい。また、初期値記憶手段は、所定のデータが記憶できればよく、種々の半導体記憶装置、種々の光学的記録装置、種々の磁気的記憶装置、を利用することが可能である。典型的には、磁気的記憶装置の一種であるハードディスク等を利用することが好ましい。
【００７５】
流跡線算出手段１４は、外部から与えられる速度場から、流跡線を算出する。流跡線の算出動作自体は、上記式（５）（「数１０」）や式（６）（「数２０」）等の式に基づき行う。特に、この流跡線算出手段１４は、式（６）に基づき流跡線を更新し、この更新された流跡線が流跡線記憶手段１６に逐次格納される。
【００７６】
また、この流跡線算出手段は、所定の演算プロセッサを利用することが好ましい。また、いわゆるディジタルシグナルプロセッサ（ＤＳＰ）を利用することも好ましい。
【００７７】
流跡線記憶手段１６は、上記算出した流跡線を記憶しておく手段である。具体的には上記の式（５）に示す関数、「数１７」の関数を記憶する。この流跡線算出手段１４は、上述した初期値記憶手段１２と同様にハードディスク等種々の記憶手段を採用可能である。但し、流跡線を表す上記関数は、式（６）のように、計算の各ステップにおいて更新されるので、書き換えが容易な手段であることが望ましい。例えば、半導体記憶手段を採用することが好適である。
【００７８】
移流計算手段１８は、流跡線記憶手段１６が保持する流跡線を表す関数と、初期値記憶手段１２に格納されている初期値と、に基づき移流を計算し、最終的に、時刻ｔｎにおける物理量の値を求めるのである。この計算は、上述した式（５）によって計算される。
【００７９】
この移流計算手段１８は、所定の演算プロセッサを利用することが好ましい。また、いわゆるディジタルシグナルプロセッサ（ＤＳＰ）を利用することも好ましい。
【００８０】
このような構成によって、本移流方程式を解く装置１０によれば、数値拡散を起こさずに移流方程式を（数値計算的に）解くことができる。
【００８１】
なお、本特許で、「解く」とは、数値計算上で解くことを意味し、具体的には時刻ｔｎにおける物理量を求めることを言う。
【００８２】
（１−３）変形例
上述した図２においては、各手段はハードウェアを用いて構成している。初期値記憶手段１２や、流跡線記憶手段１６は、データを保持する各種の手段であればどのようなものでもかまわない。また、流跡線算出手段１４や移流計算手段１８は、ＤＳＰ等、所定の演算ハードウェアを利用することが計算速度の点から望ましい。
【００８３】
しかし、対象となる空間の大きさがごく小さい場合や、高速な計算速度が要求されない場合は、汎用的なプロセッサと、このプロセッサが実行するプログラムと、から上記流跡線記憶手段１６や、移流計算手段１８を構成することも好ましい。このような構成を採用する場合には、プログラムを格納しておく記憶手段が必要であるが、初期値記憶手段１２や流跡線記憶手段１６などの記憶手段に共用することも好ましい。
【００８４】
また、図２においては、初期値記憶手段１２・流跡線記憶手段１６は別体に構成されているが、同一の記憶装置を利用することも好ましい。
【００８５】
（２−１）逆方向写像を用いた非移流段階の計算
次に、本手法において、非移流項を加える方法を説明する。
【００８６】
これまで説明したように、本実施の形態の手法では、各計算ステップにおいて得られる値
【数３８】
2006331351.tif 000040
は、常に初期値である以下の値
【数３９】
2006331351.tif 000041
を元にして構築されているため、単に上記各ステップにおいて得られる値（上記「数３８」）に非移流項を加えても、その次の計算ステップには影響がない。したがって、何らかの手法で、上記初期値である値（上記「数３８」）に非移流項による変化を加えることが必要である。
【００８７】
本実施の形態で採用する手法
第１に考えられる手法（本実施の形態では採用していない）は、非移流項を加えた
【数４０】
2006331351.tif 000042
を改めて、
【数４１】
2006331351.tif 000043
と置き、そのまま計算を続行する手法である。これを第１の手法と呼ぶ。この第１の手法は、前述した意味において、Update-Lagrangianとなり、時間発展による誤差が再び生じてしまうことが明らかである。
【００８８】
したがって、上記「数４０」での非移流項を、上記「数３６」に直接加える方法が望ましいことは言うまでもない。このような手法を第２の手法と呼び、本実施の形態ではこの第２の手法を採用している。以下、具体的にこの第２の手法を説明する。
【００８９】
第２の手法のこのような直接加える処理は、上記（５）式（「数１９」）を、移流前から、移流後の逆方向写像であると考えることによって定式化が可能である。逆方向写像とは、逆写像関数を用いて写像を求める方法であり、本実施の形態での記号（上記（５）式で採用した記法）を利用すれば、下記（７）式のように定義することができる。
【００９０】
【数４２】
2006331351.tif 000044
ここで、
【数４３】
2006331351.tif 000045
は、移流前の空間（ｔ＝０）Ω０から、移流後の空間（ｔ＝ｔｎ）Ωｎへの写像についての逆写像関数となる。また、
【数４４】
2006331351.tif 000046
は、それぞれ空間Ωｎ、Ω０で定義されているスカラー値であり、これを用いて次に述べる式（８）を求めることができる。式（８）を求めるにあたり、
【数４５】
2006331351.tif 000047
を、空間Ω０からΩｎへの写像関数を表す記号と定義する。以上のようにすれば、以下の式（８）が成り立つ。
【００９１】
【数４６】
2006331351.tif 000048
このように、空間Ωｎから空間Ω０への逆方向写像が定義可能である。つまり、
【数４７】
2006331351.tif 000049
を逆写像関数として考えた場合に、それに対応する写像関数
【数４８】
2006331351.tif 000050
を求めることができれば、移流後の空間Ωｎから移流前の空間Ω０への逆方向写像が定義できる。したがって、
【数４９】
2006331351.tif 000051
における計算で得られた非移流項を
【数５０】
2006331351.tif 000052
に加えることが可能となる。このような
【数５１】
2006331351.tif 000053
は、上述した（２）式（「数１０」）において、流跡線上では
【数５２】
2006331351.tif 000054
が一定であることを考慮すれば、時間軸についてΔｔシフトすることによって、下記（９）式が得られる。
【００９２】
【数５３】
2006331351.tif 000055
この（９）式を用いて、
【数５４】
2006331351.tif 000056
を定義した手順と同様の手順で
【数５５】
2006331351.tif 000057
を定義すれば、下記（１０）式のように定義することができる。
【００９３】
【数５６】
2006331351.tif 000058
この（１０）式では、
【数５７】
2006331351.tif 000059
は、移流前の空間Ω０から移流後の空間Ωｎへの写像関数となっていることが明らかである。時間発展についても、
【数５８】
2006331351.tif 000060
と同様に、上記写像関数（「数５７」）は、流跡線であるので、オイラー法によって時間軸で離散化すれば、上記（６）式（「数２６」）と同様に、下記式（１１）のように記述することができる。
【００９４】
（※なお、ここでは、オイラー法によって時間微分による離散化を利用したが、世の中で知られている種々の時間微分の手法を利用可能なことは言うまでもない）
【数５９】
2006331351.tif 000061
したがって、
【数６０】
2006331351.tif 000062
は、
【数６１】
2006331351.tif 000063
とは時間軸に逆方向に更新を行えばよい。このように定義された上記写像関数（「数６０」）を用いれば、後述する（１２）式によって、
【数６２】
2006331351.tif 000064
が定義されている時空間上で計算された非移流項
【数６３】
2006331351.tif 000065
を、
【数６４】
2006331351.tif 000066
に加えることができる。なお、先に述べた（１２）式は以下の通りである。
【００９５】
【数６５】
2006331351.tif 000067
この（１２）式において、
【数６６】
2006331351.tif 000068
及び
【数６７】
2006331351.tif 000069
は、非移流項を加えた結果である。なお、
【数６８】
2006331351.tif 000070
については、写像後の
【数６９】
2006331351.tif 000071
のみを用いて計算することができるので、非移流項自体の計算において写像関数を意識する必要はない。
【００９６】
非移流項を加えるさらに他の手法（第３の手法：本実施の形態では採用しない）として、以下のような手法も考えられるかもしれない。この第３の手法は別途
【数７０】
2006331351.tif 000072
を用意せずに、
【数７１】
2006331351.tif 000073
のみを用いた順方向写像を考え、これを利用する手法であり、理論的には可能である。この手法を採用した場合の順方向写像とは、下記の（１３）式
【数７２】
2006331351.tif 000074
で表される処理（によって得られる写像）である。
【００９７】
ところで、一般に、順方向写像と、逆方向写像は、数学的には（数学理論上は）等価であるが、数値計算上では、一般に等価ではない。これは、値域である
【数７３】
2006331351.tif 000075
が離散的であるので、その関数域
【数７４】
2006331351.tif 000076
が写像先を全て埋め尽くさない可能性があり、順方向写像を用いると写像後の分布に穴があくことがあるからである。したがって、この第３の手法は、理論的には可能であるかもしれないが、計算工学的的には問題が発生する。
【００９８】
これに対して、本実施の形態で提案する上記第２の手法においては、その全ての写像で採用している逆方向写像が、写像先の全ての点で処理を行っているため、写像後の分布に穴があくことはない。したがって、計算工学的には問題が発生する。この点が本実施の形態において特徴的な事項である。
【００９９】
（２−２）移流前後の空間と、関数との関係
本実施の形態において説明した上記関数Ｆ（ｘ）とＧ（ｘ）及び移流前の空間と移流後の空間との関係を示した説明図が図３に示されている。この図に示すように、関数Ｆ（ｘ）の写像の方向は移流前の空間２０から移流後の空間２２に向かう方向である。一方、関数Ｇ（ｘ）の写像の方向は、移流後の空間２２から移流前の空間２０に向かう方向である。したがって、関数Ｆ（ｘ）による正方向の写像の方向及び関数Ｇ（ｘ）による逆方向写像の方向は、ともに移流前の空間２０から、移流後の空間２２へ向かう方向である。また、関数Ｇ（ｘ）による正方向の写像の方向及び関数Ｆ（ｘ）による逆方向写像の方向は、ともに移流後の空間２２から、移流前の空間２０へ向かう方向である。
【０１００】
（２−２）実装
上記（２−１）で述べた原理に基づく移流方程式を解く装置３０の構成ブロック図が図４に示されている。
【０１０１】
この図４に示すように、移流方程式を解く装置３０は、初期値記憶手段３２と、写像関数算出手段３４と、逆写像関数算出手段３６と、写像関数・逆写像関数記憶手段３８と、非移流項加算手段４０と、移流計算手段４２と、から構成されている。
【０１０２】
初期値記憶手段３２は、図２の初期値記憶手段１２と実質的に同様の構成である。
【０１０３】
写像関数算出手段３４は、外部から与えられる速度場から、上で述べた写像関数Ｆを算出する。この関数の算出は既に説明した通り、上記式（１０）（「数５６」）式（１１）（「数５９」）に基づき行われる。なお、写像関数算出手段３４には、移流を求めたい時刻ｔｎも与えられており、写像関数算出手段３４はこの時刻に至るまで、式（１１）による更新動作をも実行する。このようにして毎回更新された写像関数Ｆが移流の計算において非移流項の加算に用いられることは今まで説明した通りである。
【０１０４】
なお、写像関数算出手段３４は、所定の演算プロセッサを利用することが好ましい。また、いわゆるディジタルシグナルプロセッサ（ＤＳＰ）を利用することも好ましい。
【０１０５】
逆写像関数算出手段３６は、外部から与えられる速度場から、上で述べた逆写像関数Ｇを算出する。この関数の算出は既に説明した通り、上記式（５）（「数１９」）式（６）（「数２９」）に基づき行われる。なお、逆写像関数算出手段３６には、移流を求めたい時刻ｔｎも与えられており、逆写像関数算出手段３６はこの時刻に至るまで、式（６）による更新動作をも実行する。このようにして毎回更新された逆写像関数Ｇが移流の計算に用いられることは今まで説明した通りである。
【０１０６】
なお、逆写像関数算出手段３４は、所定の演算プロセッサを利用することが好ましい。また、いわゆるディジタルシグナルプロセッサ（ＤＳＰ）を利用することも好ましい。
【０１０７】
写像関数・逆写像関数記憶手段３８は、上記算出した写像関数、逆写像関数を記憶しておく手段である。この写像関数・逆写像関数記憶手段３８は、図２の流跡線記憶手段１６と基本的に同様の構成である。
【０１０８】
非移流項加算手段４０は、写像関数Ｆを用いて非移流項を加算する。この動作は、上記式（１２）（「数６５」）による計算動作である。この非移流項加算手段４０も、所定の演算プロセッサを利用することが好ましい。また、いわゆるディジタルシグナルプロセッサ（ＤＳＰ）を利用することも好ましい。
【０１０９】
移流計算手段４２は、基本的に図２における移流計算手段１８と同様の構成である。異なる点は、初期値に非移流項が加算された初期値を用いている点である。この相違点以外は、上記移流計算手段１８と同様の構成である。
【０１１０】
この移流計算手段４２は、所定の演算プロセッサを利用することが好ましい。また、いわゆるディジタルシグナルプロセッサ（ＤＳＰ）を利用することも好ましい。
【０１１１】
このような構成によって、本移流方程式を解く装置３０によれば、非移流項を含むような移流方程式であっても、数値拡散を起こさずに移流方程式を（数値計算的に）解くことができる。
【０１１２】
（２−３）変形例
上述した図４においては、各手段はハードウェアを用いて構成している。初期値記憶手段３２や、写像関数・逆写像関数記憶手段３８は、データを保持する各種の記憶手段であればどのようなものでもかまわない。また、写像関数算出手段３４や逆写像関数算出手段３６や、移流計算手段４２は、ＤＳＰ等、所定の演算ハードウェアを利用することが計算速度の点から望ましい。
【０１１３】
しかし、対象となる空間の大きさがごく小さい場合や、高速な計算速度が要求されない場合は、汎用的なプロセッサと、このプロセッサが実行するプログラムと、から上記写像関数算出手段３４や逆写像関数算出手段３６と、さらに移流計算手段４２を構成することも好ましい。このような構成を採用する場合には、プログラムを格納しておく記憶手段が必要であるが、初期値記憶手段３２や写像関数・逆写像関数記憶手段３８などの記憶手段と共用することも好ましい。
【０１１４】
また、図４においては、初期値記憶手段３２と、写像関数・逆写像関数記憶手段３８とは、別体に構成されているが、同一の記憶装置を利用することも好ましい。
【０１１５】
（３）写像関数の再初期化
これまで述べた本実施の形態の手法（上記第２の手法）は、理論的には初期値を保ったまま、計算を無限に進めることが可能であり、その場合でも誤差は累積しない。
【０１１６】
しかしながら、実用上は、写像関数にある程度のゆがみが生じてしまった時点で、その時点の場を初期値として置き換える処理を行った方が好ましい結果が得られる場合が多い。この処理（再初期化と呼ぶ）は、以下の（１４）式のように定式化される。
【０１１７】
【数７５】
2006331351.tif 000077
つまり、現在の
【数７６】
2006331351.tif 000078
を、
【数７７】
2006331351.tif 000079
として置き換えて、
【数７８】
2006331351.tif 000080
と、
【数７９】
2006331351.tif 000081
とを
【数８０】
2006331351.tif 000082
に初期化するという処理を行う。既に述べたように、
【数８１】
2006331351.tif 000083
を求めるための
【数８２】
2006331351.tif 000084
は、
【数８３】
2006331351.tif 000085
の格子点上には、存在しないことがほとんどであるので、何らかの補間が必要となる。この補間による誤差は、その後の計算には残存しないので、誤差が累積してしまうことがないことは既に説明した通りである。
【０１１８】
しかし、上述したような再初期化を行うと、補間による誤差が計算結果に加わってしまう。つまり、初期値に誤差が累積してしまう。そのため、上記初期化はなるべく行わない方が望ましい。そして初期値をなるべく長期間保ったまま計算を進めることが望ましい。このようなことを実現するための手法としては、テクスチャマッピングによる流れ場の可視化において研究されている手法が応用できると考えられる。このような手法は、上記非特許文献３に記載されている。
【０１１９】
第３章  計算例
本実施の形態で提案する手法を用いた結果を説明する。以下説明する例においては、全て、空間１次精度の補間と、時間２次精度の写像関数の更新を用いている。
【０１２０】
まず、図５には、本実施の形態で提案する手法を用いて２次元の移流方程式を解いた結果（ｔ＝３１４）と、その初期値が示されている。ここで示す「値」は所定の物理量であればどのようなものでもよく、また既に説明したように波だけでなく物でもよい。水中の絵の具の濃度や、温度等、種々の物理量に適用することができる。
【０１２１】
まず、図５（１）には、この初期値が示されており、２次元上のグラフで初期値が示されている。初期値は各グラフの「高さ」で示されており、その結果、図５（１）のグラフはいわゆる３次元グラフとして描かれている。図５（１）の「平面」は、移流方程式が適用される２次元空間をそのまま表している。
【０１２２】
図５（２）には、図５（１）の初期値から一定時間（ｔ＝３１４）経過した後の平面上の値を示すグラフが示されている。このグラフも、２次元の空間（平面）とその空間上の値を「高さ」で示す３次元グラフである。
【０１２３】
この例において、速度場は、
【数８４】
2006331351.tif 000086
と与えられる。格子の解像度は５０＾２である。すなわち、上記平面には５０×５０の格子が設けられている。２次元空間（平面）のこの格子は、物理量の値に応じた高さを持つように変形した結果、図５（１）（２）に示すようなグラフが形成される。
【０１２４】
また、Δｘ＝Δｔ＝１としている。移流対象の分布の初期値は、
【数８５】
2006331351.tif 000087
である。この場合、初期値の分布が、その初期の分布形状を保ったまま、（ｔ＝１００π）で初期値に一致するのが厳密解である。
【０１２５】
ここで示した計算例（図５）では写像関数の再初期化は必要ないが、比較のため、同様の例に関して３０ステップごとに再初期化を行った場合の結果を図６（２）に示す。また、従来の３次精度Semi-Lagrangian法での結果を図６（１）に示す。
【０１２６】
次に、解析解と従来手法であるSemi-Lagrangian法との誤差評価と、解析解と本手法（１次精度）との誤差評価と、を比較する表が表１に示されている。この表１に示すように、本実施の形態の手法によれば、誤差評価の値が従来より小さな値となる。このことは、図５からも明らかであり、本実施の形態の手法によれば再初期化を頻繁に行った場合でも、従来の３次精度の計算結果に比べて、視覚上、分布形状の「なまり」が少ないという結果が得られた。
【０１２７】
【表１】
2006331351.tif 000088
次に、図７は、非移流項を含む一般化された移流方程式を従来の手法と、本実施の形態の手法で解いた結果を示すグラフである。
【０１２８】
この図に示す例においては、速度場や、空間、時間刻みについては図５や図６と同様である。但し、分布の初期値は０とした。また、非移流項は、
【数８６】
2006331351.tif 000089
というものを採用した。このような非移流項の場合、厳密解は非移流項による分布が、１格子幅の線による円を描くことになる。図７の（１）は、３次精度の従来の手法（Semi-Lagrangian法）によるグラフであり、横軸がｘ軸、縦軸がｙ軸の平面を表し、グラフ上、値の大きさはグラフ中の色の濃さで表されている。図７の（２）は同様の計算を本実施の形態で提案する手法で計算した結果を示すグラフであり、縦軸横軸は図５の（１）と同様の意味である。
【０１２９】
第３章  本実施の形態の効果とまとめ
図５に関して
まず、図５に関しては、１次精度による本実施の形態の手法は、オーバーシュートやアンダーシュートを生じさせることがなく、厳密解に非常に近い結果が得られることが理解されよう。初期値として、明確な境界を有する分布と、滑らかな分布と、の２種類の分布を含む場合を採用したが、本実施の形態の手法によれば、特段の取り扱いをすることなく、高精度に「移流」を求めることが可能である。
【０１３０】
図６及び表１に関して
次に、図６及び表１に関しては、以下のことが判明した。
【０１３１】
すなわち、本実施の形態の手法で、再初期化を行った場合においても、３次精度の従来手法より高い精度を実現することができる。特に、本実施の形態の手法が１次精度であるにもかかわらず、３次精度の従来手法より高い精度を実現することができた点に注目すべきである。
【０１３２】
図７に関して
次に、図７に関しては、本実施の形態の手法では、非移流項がある場合でも非常に精度よく移流が計算できることが理解されよう。
【０１３３】
なお、本実施の形態においては、いずれの場合も、境界に対して特別な取り扱いをせず、単純に本実施の形態で提案する移流方程式の解法を適用して移流を計算した。
【０１３４】
従来の手法では、境界に対して特段の処理や取り扱いをしなければ良好な結果が得られなかったが、本実施の形態の手法によれば、そのような特段の取り扱いをしなくとも、良好な結果が得られた。例えば、図４（２）における円の図形も、暈（ぼ）けずに描かれていることは、本実施の形態の手法によれば高精度に移流方程式を解くことができるということを如実に示していると言えよう。
【０１３５】
まとめ
以上述べたように、本実施の形態によれば、移流方程式を移流前から移流後への写像として定式化し、各計算ステップごとに、移流後の分布を再構築する手法を提案した。この手法によれば、時間発展によって蓄積する誤差を抑制することができる。また、非移流項については、順方向写像と逆方向写像の考え方を導入することによって、写像に起因する誤差の発生を防止することができた。
【０１３６】
また、本発明の特徴の一つは、上記Ｇの計算については，速度場を反転させたもので粒子を追跡し，上記Ｆの計算については，速度場をそのまま用いて粒子を追跡するという点であることに留意されたい。
【０１３７】
また、本実施の形態で述べた定式化は、従来の移流方程式の解法とは全く異なるものであり、構造解析におけるTotal-Lagrangianの概念を、移流方程式で実現したものである。
【０１３８】
したがって、本発明は、移流方程式の新たな解法の方向性を世の中に対して示すものであり、極めて斬新なものであると言えよう。
【０１３９】
発明性について
なお、本実施の形態で説明したように、本件発明は、空間上の所定の物理量が「流れ」によって移動していくという物理法則、又は、波動が空間を伝播していくという自然法則に基づいて方程式を解く技術的なアイデアであり、特許の対象となる技術思想であることは明らかである。
【０１４０】
また、本実施の形態で説明したように、本件発明は、純粋な机上の数学理論ではなく、実際に値を求める場合に誤差や数値拡散を防止し、精度の高い計算を行おうとする工学的なアプローチであり、エンジニアリング上の工夫であるので、自然法則を利用した技術的な思想であることは明らかである。
【図面の簡単な説明】
【０１４１】
【図１】従来手法と本実施の形態の手法との原理の相違を説明する説明図である。
【図２】本実施の形態における移流方程式を解く装置の構成ブロック図である。
【図３】本実施の形態における関数Ｆ（ｘ）とＧ（ｘ）の写像の方向に関する説明図である。
【図４】本実施の形態における移流方程式を解く装置の構成ブロック図である。
【図５】本実施の形態ので提案する手法を用いて２次元の移流方程式を解いた結果（ｔ＝３１４）と、その初期値を示すグラフである。
【図６】図５の例において、３０ステップごとに再初期化を行った場合の結果と、従来の３次精度Semi-Lagrangian法での結果を示すグラフである。
【図７】非移流項を含む一般化された移流方程式を従来の手法と、本実施の形態の手法で解いた結果を示すグラフである。
【図８】「移流」の概念を表すグラフである。
【図９】従来のＣＩＰ法の原理を示す説明図である。
【符号の説明】
【０１４２】
１０  移流方程式を解く装置
１２  初期値記憶手段
１４  流跡線算出手段
１６  流跡線記憶手段
１８  移流計算手段
２０  移流前の空間
２２  移流後の空間
３０  移流方程式を解く装置
３２  初期値記憶手段
３４  写像関数算出手段
３６  逆写像関数算出手段
３８  写像関数・逆写像関数記憶手段
４０  非移流項加算手段
４２  移流計算手段
【図８】
2006331351.tif 000096
【図１】
2006331351.tif 000097
【図２】
2006331351.tif 000098
【図３】
2006331351.tif 000099
【図４】
2006331351.tif 000100
【図５】
2006331351.tif 000101
【図６】
2006331351.tif 000102
【図７】
2006331351.tif 000103
【図９】
2006331351.tif 000104