極座標を用いて数値計算する問題について

極座標は境界形状が丸いものに対して適用するということになりますが、最も簡単な事例は円だと思います。そこで円という境界に囲まれた内部領域の物理現象（流体とか固体とか）を考えます（２次元空間）。その内部を極座標(r, θ)という座標で領域分割（r_n=ndr, θ_m=mdθ)して差分法でも有限要素法でも適用できると思います。ここでどうしても気になるのは円の中心です。これは極座標の特異点ではないでしょうか（座標軸が集中するとか、物理方程式で1/rが出てきたときのr=0だとか）。しかし、実際は円の中心は境界（円周）からも遠いし、いちばん平凡なところともいえると思います。つまり、極座標を導入することによって特異点が導入されてしまうのではないかと思うのですが、どのように対処するものでしょうか。このような問題で一番いい方法は通常の直交デカルト座標で形状に合致するように非構造格子をつくって（これもやり方がいろいろあるのでしょうが）普通の有限要素法を適用することだろうと思います。しかし、特異点のある極座標の場合では計算や理論展開においてどのようにしてこの問題（円の中心の特異性）に対処するのでしょうか（もともと問題などない（特異点ではなかった）となるでしょうか。）

地球科学では球座標を使う（この場合はほぼこの方法しかない？）と思いますが、両極は特異点ではないかと思うのですが。あるいは球面（直交デカルト座標のデータとしては３次元）を任意の三角形で被覆したときの座標データはないでしょうか。球なので誰が作っても同じものになるし、球の半径の違いによる値の変換は簡単だと思うのですが。

量子力学、分子物理のようなところでも球座標は出てくるはずですが。

よろしくお願いします。

No.1 回答者： yhr2 回答日時： 2019/07/10 14:10

＞最も簡単な事例は円だと思います。

2次元であればそうです。

＞地球科学では球座標を使う（この場合はほぼこの方法しかない？）と思いますが、両極は特異点ではないかと思うのですが。

それは3次元であって、「両極（北極、南極）」は球の中心ではありませんよ？
万有引力は、地球の中心では無限大になるので特異点です。

この回答へのお礼

回答ありがとうございます。両極は経線が集中してくるので特異点になるのではないでしょうか（経度方向に値が決まらない？）。球面という平面上のことですが。万有引力は距離の２乗分の１という式になるので距離がゼロのときはどうするかということはありますね。平面の極座標ではあちこちに1/r（ｒで割る）というのが出てくるので何らかの措置が必要だと思うのですが。球面の場合はルジャンドル関数を使うとかスペクトル的な方法（関数展開してその和で解を構成する）で問題を回避しているのではないかと思います。式を満たすことが分かっている関数の重ね合わせで解を構成するので式に関する問題は回避されるように思います。とにかく式に示されたとおりにプログラムしようと思ってもできない気がするのですが。

お礼日時：2019/07/10 16:58

No.2 回答者： tknakamuri 回答日時： 2019/07/11 08:14

なんか論理がメチャクチャ

式が1/rを含んでいて原点で発散するなら
デカルト座標でもそれが1/√(x^2+y^2)になるだけで
発散することには変わりはない。
特異点を持つということは座標系とは無関係。

極座標で少し厄介なのは原点で角度が定まらない
ことだけど、実用上問題になることは少ない。

この回答へのお礼

ご指摘ありがとうございます。確かに、そういうことになりますが、偏微分方程式において∇、ラプラス演算子、発散演算子、回転演算子というものは1/rを含んでしまうことがあり、それでいて状況としてr=0が特別な点でない（円筒というなら缶詰の真ん中）ということとの乖離はあると思います（解釈に困る）。また、座標軸が1点に収れんするという性質は円筒、球面でも起こっており、難しいように思うのですが（これは式系ではなく座標系としての特異点とは言わないのでしょうか）。今回、そういうのをいっしょくたにして特異点とか特別な点と言う風に書いています。
実際の問題ですが、プログラミングする場合、どう処理すればいいでしょうか（プログラミングという実用上の処理）。
例えば、2次元の(r,θ)空間の極座標は0<r<R, 0<θ<2πという領域を考えるのでr=0は境界（端っこ）ともいえるように思います。でも極座標（円筒形の缶詰）では境界は円周上であり、r=Rのところのはずですね。これも何となく理解しにくいところがあります。缶詰の真ん中(r=0)は境界ではないですよね。

お礼日時：2019/07/12 15:38

No.3 回答者： eatern27 回答日時： 2019/07/12 01:36

お考えのように極座標系のような必ずしも全単射ではない（ヤコビアンがゼロになる点がある）座標変換を考えると、見かけ上特異点のように振る舞う事があります。

原点に点電荷を置いている場合のように実際に特異点になっている場合を考える事もあるので、一概にこのように対処するというような言い方はしにくい部分ではあるのですが、
原点や両極などが真の特異点ではないのであれば、
基本的には「得られる解が特異的な振るまいを示さない」という条件をあらわに課す必要が出てきます。
より具体的には具体的には原点や両極で何らかの境界条件を課すことになります。

一例として極座標系でスカラー場φについての微分方程式を考えているような場合であれば、
「φ(r=0,θ)がθに依らない定数になる」の辺りが条件になるはずです。
物理では∇^2φを含む微分方程式が登場する事が多いですが、こういう場合
r=0で∇^2φが有限の値であるのならgradφが連続でなければならないので、∂φ/∂r=0の条件も課される事になります。


質問にある数値解析云々の話は無視していますが、メッシュを切った事によって初めて発生する事を問題視しているようには思えなかったのでこんな所で回答になるでしょうか。


>>#2さん
ラプラシアンなどの微分演算子の極座標表示に1/rが登場するという事が念頭にあるのでしょう。

この回答へのお礼

回答ありがとうございました。数値解析を念頭に置いており、缶詰に入っている水の運動を考えてみます。缶詰の真ん中は全く特別な場所ではありません。もっとも平凡な場所（缶詰の側壁か遠く境界の影響が遠い）ともいえると思います。
もう少し具体的にいうと、極座標は0<r<R, 0<θ<2πの区間をdr, dθで区切ったメッシュで考えます（缶詰の半径がR）。この時のメッシュはバームクーヘンとかパイナップルの缶詰を取り出してカットしたようなイメージです。一番内側のメッシュは0<r<dr, ndθ<θ<(n+1)dθ となっていますが、そのメッシュの内側の面(r=0側)でのUr(r方向流速）をどう与えたらいいでしょうか。そこでは流速を定義しないで済むようにして計算するということでしょうか。例えば∂Ur/∂r=0となるように計算するとして、r=0でのUrがいっぱいある（θ方向メッシュ分）ことの違和感には目をつぶるということでしょうか。
r=0がない、芯をくりぬいたパイナップル（2つ円筒に挟まれた領域）のようなところの流体の計算だったら全く問題ないのですが。

お礼日時：2019/07/12 16:04

No.4 回答者： eatern27 回答日時： 2019/07/13 16:43

具体的なアルゴリズムなんていくらでもあるので想定されているアルゴリズムに応じて

適切な原点などの処理をお考えいただくしかありませんが、

#3に書いた事というのは、r=0の所に何らかの境界条件を課すことにすれば、
「r,θが極座標の座標である」という事を忘れて、ただの2変数関数の偏微分方程式だと思って解く事ができる
という事です。（正確にはθ方向の周期境界条件も必要になりますが）

逆に言えば境界条件さえ適切に決めてしまえば、あとは通常のデカルト座標の偏微分方程式の話と何も変わらなくなります。

>一番内側のメッシュは0<r<dr, ndθ<θ<(n+1)dθ となっていますが、そのメッシュの内側の面(r=0側)でのUr(r方向流速）をどう与えたらいいでしょうか。
どう考えるかはいくらでもあるでしょうが、
一番内側のメッシュの範囲をε<r<drだと思って、ε→+0の極限だと思う事にするとか
原点だけ速度場をUxとUyなど基底ベクトルを別途定義する
みたいな感じになるでしょうかね。

>例えば∂Ur/∂r=0となるように計算するとして、r=0でのUrがいっぱいある（θ方向メッシュ分）ことの違和感には目をつぶるということでしょうか。
そこが気になるのであれば、原点を特別扱いしてr=0が一つの点となるように定式化すればよいだけです。
この場合#3に書いたようなr=0での境界条件と言うのは出てこなくなるでしょう。


ただ、本当にここでお聞きになっているような事を考えなきゃいけないような状況であるのなら、そもそも極座標系で考えるのが良かったのかどうかを検討するのが良いように思います。

この回答へのお礼

回答ありがとうございます。r=0での処理の仕方がいろいろあるというご指摘だと思います。でもやはり円筒の中心は一番なにも起こらないところ、何もトリックが必要ないところとも言えると思うのです。いろいろやり方があるということはその選択の仕方でいろんな解がある（似通ったもののはずですが）ということを容認してしまいそうなのですが。
数値計算という立場に立つと円周でのギザギザを恐れずに碁盤の目のような直交デカルト座標でやるという方法もあるようです。ギザギザを少しでも緩和するために円周付近の格子の正方形を4等分するとかさらに4等分するとかです。四分木格子と言うらしいですが。それも方法ですかね。これは円筒座標、極座標の利用をあきらめるということですが。

お礼日時：2019/07/14 19:58

No.5 回答者： masa2211 回答日時： 2019/07/15 10:23

気象シミュレーション用FEMのメッシュの切りかたですが、

現状のもの（気象庁）は
https://www.jma.go.jp/jma/kishou/books/nwptext/4 …
GSM:標準ガウス格子
RGG：適合ガウス格子

次世代型は
http://www.jamstec.go.jp/j/about/press_release/2 …
DFSM：適合ガウス格子を精密化
NICAM：ほぼ正三角形
MSSG：サンソン（極地方抜き）＋横サンソン（極地方）で、ラップ部ほ処理は、よくわからん。

GSM以外に考え方は、
・各メッシュの面積は概ね同等。
を狙っているようです。（うん、まあ、そうしたほうが情報落ちしにくいから計算精度をキープするために式を変な風にいじる必要性が少なくいなる。）
で、質問に対し明確なのは
・座標表記法は関係ない。（個別メッシュの寸法、メッシュ間の接続がわかれば、それ以外の情報は不要となるので。）
・極部分では、南北方向の物理量の出入りはゼロ。
　　DFSMの図を参照。極地方のみ、扇型の先端が尖っているため、南北方向の熱の出入りを式に表現できない。

No.4
>(r=0での処理の仕方が)いろいろやり方があるということはその選択の仕方でいろんな解がある
>ということを容認してしまいそうなのですが。
FEMは、要素の切り方がいろいろ（三角形とか四角形とか、辺の物性分布が一定なのか直線なのか２次式なのか）
あり、その結果いろいろな解があるけど、それを容認しています。
要は、必要な精度（と計算時間）に達していれば何でもアリ。

質問文
＞地球科学では．．．任意の三角形で被覆．．．球なので誰が作っても同じものになるし、
＞球の半径の違いによる値の変換は簡単だと思うのですが。
ここでいう球の半径とは、地表から地下あるいは空中で向かう物理量のこと、ですよね？
それって地球中心までシミュレーションに組み込むことはありえないから、
当初質問の「半径ゼロ部分の物理量が表現できない」と関係ないです。
ですので、物理量が表現できない場合と、単に座標で表現できない場合を一緒くたにしてほしくないです。
両者で別回答となる場合があるので。

No.３
＞缶詰に入っている水の運動を考えてみます。缶詰の真ん中は全く特別な場所ではありません。
缶詰を熱湯で暖めた場合の熱の伝わりかた　のようなことを考えると想定します。
真ん中を無視し、パイナップル形状で何故ダメなのか、必要に応じ、中央の穴を小さくすれば十分な近似になるだろうとしか言えません。
で、缶詰の水の運動とあるので、内部が固体（サバ缶など）でなく、ジュースなど、だよね？
この場合、缶の中で対流がおきて、たとえば、缶の中央が下降流、缶の周囲が上昇流となったら、
缶の中央より、中央と周囲の中央付近、が温度が低いということに。
（どういう対流が起きるかのよるので、この限りでは無いけど。）

で、缶詰の「固体」での熱伝導（凍ったジュース缶のように、小口は無視できる条件で）なら、
別に２次元で解く必要は無く、１次元要素で間に合う（中心側は、勿論、断熱点として考える。）と思うし、
ひょっとして、管路（パイプライン）における水流の流速分布？
うん、管路の中央が最も流速速いから、中央において半径方向の流速が賞実ことは考えなくてヨシ。

No.4
＞数値計算という立場に立つと．．．．．
＞円周でのギザギザを恐れずに直交座標．．．．極座標の利用をあきらめるということですが。
意味解りません。
円周でのギザギザを作るくらいなら、円中央では、半径方向の物理量のやりとり無し（扇型のFEMで普通に組めば必然的ににそうなる）または、中央に微細な穴を開け、パイナップル型にする
ほうが、得られる精度が高いから、そうすると思うが。

この回答へのお礼

懇篤な回答ありがとうございます。

＞極地方のみ、..南北方向の熱の出入りを式に表現できない。

これは計算上OKでしょうか。熱だけでなく、運動量や質量も同じ意味で通過できないのでは？と思うのですが。

＞FEMは、要素の切り方がいろいろ.. あり、その結果いろいろな解がある..

解がやや異なるのは容認できますが、大きく異なるならどっちを信じるのか判断する必要があります。メッシュの切り方で解が異なることを容認するには限度があるとは思うのですが。

＞ここでいう球の半径とは、.. それって地球中心までシミュレーションに組み込むことはありえない..

地震波の場合中心はあり得ますね。大気力学は地表面以上です。地震波計算での地球の中心、大気力学での両極に関する疑問です。
　両極は太陽の熱流入が最低なので特殊な点でしょうか。SF映画で隕石が衝突して熱波や地殻の変動が地球の１点から球面上を同心円状に広がるのがあります。あの場合、両極は特別な点ではないですね。

＞別に２次元で解く必要は無く、１次元要素で間に合う..と思うし、..うん、管路の中央が最も流速速いから、..考えなくてヨシ。

条件から結果を限定して、解を求めるという立場は確かにありますね。しかし、条件が1次元的でも現象が2次元的だったり、2次元の条件で現象が3次元的な場合があります（気象力学の傾圧不安定）。条件が1次元だから数値実験も1次元でよいとは言えないはずですが。

＞円周でのギザギザを恐れずに直交座標..極座標の利用をあきらめるということですが。
＞意味解りません。

最近の動向で直交デカルト座標を用い、境界が複雑なところ（円周も）は座標を2分割（正方形4分割）して局所的に格子を細かくします（四分木格子）。足りないならそれをさらに分割します。方程式を全く単純なまま使え、計算量も少ないです。曲線座標系では項が増えます。しかも、最近はその格子分割を解に適合させ、計算結果に従って動的に変化するという技法(アダプティブメッシュリファインニング(AMR)）も出てきて幅広く使われ始めたようです。
　一方、円筒座標では波の等位相線（波峰線）が屈託なく円の中心を通過するとはどうしても思えないです（SF映画の隕石衝突後の衝撃波の波峰線）。直交デカルト座標では何の問題もなく缶の中心を一本の等位相線が通過できそうです。

お礼日時：2019/07/16 18:52