偏微分方程式の数値解法

質問

偏微分方程式の込み入った質問です。
2次元(x,y)の空間で2つの関数f(x,y),g(x,y)を考えます。

そこで、それぞれにラプラス方程式を立てました。
fxx+fyy = 0　 (1)
gxx+gyy = 0 (2)
です。これは境界値問題で、差分式からSOR法を使って収束計算によって数値解を求めることができます。f, gはそれぞれ独立という形にはなります。
そこにもう1つ式が出てきました。
fxfy + gxgy = 0 (3)
というものです。f,gをx,yで1回微分してできる式です。
都合3つの式が出てきました。

この数値解を求めるにはどのような方法があるでしょうか。
数値解ですから近似解です。
3つ目の拘束条件の下でのラプラス方程式とみると、ペナルティ関数とかラグランジュの未定係数法とかいろいろあるかもなと思いますが。

3つ目の式は完全に満たすというより、できるだけ満足するようにしたいというものです。

よろしくお願いします。

回答No.3

ddtddtddt

#2です。境界要素法はご存知でしたか(^^)。

＞ラプラス方程式はパラレルの2回微分ですが、クロスの微分（2方向なのでそれも2回微分）項を入れてはどうかということです。

　ラプラス方程式を満たす事：調和関数である事は非常に良い性質なので、それは捨てたくないなぁ～、というのが個人的な印象です(^^;)。収束計算という事で言えば、次のような方法はあると思います（上手く行く保証はすぐには出来ませんが）。

　前回の式で最初に(fj)と(gj)を与えて整理すれば、添付図の最初の3本になりますが、未知数は(fj(n))と(gj(n))の2nに対して条件数は3nとなり、条件過多の系です。そこで左辺の行列の転置を両辺にかけて条件数2nに縮約し、最小2乗型の最適解として(fj(n))と(gj(n))を求めます。

　次にそうやって求めた(fj(n))と(gj(n))を初期値として整理すれば、後半の3本になるので（(fjgj)は、初めの2つから直接計算して代入）、同様に最小2乗型の最適解として(fj)と(gj)を求めて最初に戻る。

　・・・なんとなぁ～く、収束してそれなりに厳密解が求まりそうな気はするんですが(^^;)、望んだ解に制御できるか？という点は心配です。このあたりは、どのようなものが欲しいかで決まるので、自分としてはこれ以上の事は言えません。

　余談ですが、メッシューとしてはアドバンシング・フロントと組み合わせて使用するイメージでした。要素の潰れを防止できるんじゃないかと。もっとも業務では、市販ソフトに頼りますけれど・・・(^^;)。

回答No.2

ddtddtddt

#1です。ちょっと風邪ひいてました(^^;)。

　以下は前回の返答のイメージで書いてますので、誤解かも知れませんが・・・。

　まずあなたのやりたい事は、一定の領域Rの境界上でf，gの境界値（曲線座標系での座標値）を設定し、R内部で∇f・∇g＝0が満たされるように内部関数値を修正していく、という風に読めました。

　しかし上記はたぶん無理です。理由はラプラス方程式の性質から、境界上の関数値を全て与えれば、ラプラス方程式を満たす限り領域内部の関数値も一意に決まってしまうからです。

　勝手に与えたf，gの境界値の解が、R内部で∇f・∇g＝0を満たすとはもちろん限らないので、R内部で∇f・∇g＝0が出来るだけ満たされるように内部関数値を修正すると、ラプラス方程式を満たす事も諦めた事になります。

　よって最初に与える境界条件が重要になります。変動する境界条件を扱う場合、やはり境界要素法の見通しが良いので、以下それで記述させて下さい（添付図）。境界要素法は古典的なグリーン関数法を、数値計算用にアレンジしたものです。「ラプラス方程式の境界要素法」などで検索すれば、式(1)などの導出は、いくらでも出てきます。

　式(1)は、境界方程式と呼ばれます。φがラプラス方程式を満たすべき解関数で、φ*はラプラス方程式の基本解と言われ(3)，(4)です。∂/∂nは、外法線微分値を表します。

　式(1)の全てのパラメータは解析領域Rの境界S上のもので、例えばRを右上の図のように折れ線近似する事で、容易に離散化し連立一次方程式に持ち込めます。右上の図はφ＝fとしたものです。

　折れ線が境界要素に当たり、要素端節点の未知量：fjと∂fj/∂nで要素上の関数値を直線近似するなどして境界積分を離散化します。ki（k）は節点iの内角を表し、滑らかな境界の折れ線近似なら、ki/(2π)＝π/(2π)＝1/2とすれば十分です。基本解のrは節点iと積分点の距離で、θはrの方向，βは積分する境界要素の外法線方向です。

　式(1)を離散化すると、式(5)のようになります。diagはki/(2π)を表す対角行列，Bは∂φ*/∂nに関する積分の離散化結果，Hにはφ*に関する積分の離散化結果が入ります。境界S上にn個の境界要素があるとすると境界節点もn個なので、diag，B，Hはn×nです。一方、境界未知量は境界値(fj)とその法線微分値(fj(n))の2n個です。このままでは不定解になるので、通常はここにn個の境界条件を投入します。

　例えばそれをfのn個の境界値(fj)だとすると、式(5)はn個の法線微分値(fj(n))に関する方程式となり、連立一次方程式解くことで(fj(n))が決まります。(fj)と(fj(n))が決まった後で、式(2)にそれらの結果を使うと、領域R内任意点の関数値を計算できる事になります。式(2)の(ξ，η)は、領域内の任意の内点です。

　以上より、境界上の関数値を決めるとΔf＝0を満たすfは、一意に決定される事になります。

　次にfとgの積fgが、ラプラス方程式を満たすとします。直接微分すると、

　　Δ(fg)＝f Δg＋2∇f・∇g＋g Δf＝0

が得られるので、Δf＝Δg＝0も成り立つとすると、上式より、

　　∇f・∇g＝0

を導けます。よってΔf＝Δg＝Δ(fg)＝0が、望む条件の一つの形です。これを境界要素法で定式化すると、添付図の(5)，(6)，(7)になります。((fjgj)(n))は(fj)，(gj)，(fj(n))，(gj(n))で表せる事に注意して未知数の数と条件数を勘定すると、未知数4n個に対して条件数3nです。

　したがって不定解になりますが、必ずΔf＝Δg＝∇f・∇g＝0を満たす解が存在する事にもなります。しかも関数1個分の自由度が与えられます（4n－3n＝nだから）。ただ(7)は、(fj)と(gj(n))，(gj)と(fj(n))の積が未知数に入ってくるので、このままではかなり難儀です（2次形式を対角化しないと解けそうにない）。そこで関数1個分の自由度を利用して、トライアンドエラーしたらどうだろう？と思いました。厳密解は必ずあるのだから。

　例えばfの境界値(fj)を全て与えて(5)から(fj(n))を計算し、それらを(6)，(7)に代入すれば、未知数2nに対して条件数2nの普通の連立一次方程式です。(fj)の分布を何らかの基準で系統的に変えてみて解の挙動を調べれば、望む形の直交曲線座標系を与えそうな(fj)の分布は予想がつくような気がします。こちらの方が座標系の形を制御するという観点からは、便利そうです。それに(5)，(6)，(7)を不定解の形で直接解いたとしても、解を一意化するために、やっぱり似たような事をやるだろうと思われます。

　以上、自分だったらこうしそうだ、という話でした(^^;)。

　余談ですが、むかし有限要素法の自動メッシュが欲しかった時に、等角写像を用いてマップするという方法がありました。正方形領域に正方形メッシュを切って、それを等角写像で任意領域に写像し、直交分割網を作成するという方法です。等角写像なので、当然共役な調和関数を求めるために連立ラプラス方程式という事になります。やりたかった事は、たぶん今のあなたと同じです。でもその方程式系が余りに複雑だったのであきらめましたが、今回考えてみて、こうすれば案外実用的かな？と思えて来たので、やってみようかな？という気になってます(^^)。

お礼

回答ありがとうございます。年もあらたまりまりました。風邪をひかれたとのことでどうぞご自愛を。
私は以前、境界要素法を使いましたが、ラプラス方程式を満たす基本解がある、というところが出発して境界の形状等の変化に応じて基本解係数を等を調整していく、というような感じでしたが（数学的厳密性が高い）。またFEMのメッシャーでは解が急激に変化するところで要素分割を多くするということで私の中では曲線座標の直交性に生かせるイメージが付きにくい感じです。
あと、等角写像については調和関数を使うようなので境界要素法に近いかなと思います。

私の考え方はもう少しラフで、解が直交になるように変更を加えるアルゴリズムがあればいいというところがあります。適当な初期条件でラプラス方程式を解いた→直交性はどうだ？→不十分→直交からのずれを使って修正→それを初期条件にして再度ラプラス方程式を解く　という繰り返しです。繰り返すたびにいい方向に行くので適当なところでストップとなります。

と、考えてきて必ずしもラプラス方程式でなくてもいいかもしれないという気になりました。ラプラス方程式はパラレルの2回微分ですが、クロスの微分（2方向なのでそれも2回微分）項を入れてはどうかということです。

※大長考になりますと、このスレッドの時間切れになりそうなのですが。

回答No.1

ddtddtddt

ラプラシアンをΔ＝(∂/∂x)^2＋(∂/∂y)^2で表します。まず、

　　Δf＝0　と　解を定める完全な境界条件．　　(1)

で、fは一意に決まります。g(x，y)についても同様なので、あなたの言うようにf，gはそれぞれ独立に定まり、∇＝(∂/∂x，∂/∂y)として必ずしも、

　　∇f・∇g＝0　　　(2)

を満たすとは限りません（・は内積）。というか、そうならない場合が普通だと思います。fまたはgの境界条件を、(2)を満たすように相手に合わせる必要が生じます。

　ところで(2)を満たすという事は、もしfが決まっていれば、gは(2)から本質的には決められる事になります。(2)の意味は、gの等高線（特性曲線）はfの等高線に直交する、という事なので、gの一般的形は、z＝g(gの等高線の方程式)という形になるからです。gの等高線の方程式＝c（任意定数）に対する関数g(c)の具体的な形を決めるためには、gに関する何らかの境界条件は必要です。

　状況はよくわかりませんが、fとgが(2)を満たすならコーシー・リーマン条件を満たすので、共役な調和関数を定める問題という事になります。

　それから自分はラプラス問題では境界要素法を良く使うのですが、fとgが(2)を満たすなら、添付図のような方法はあると思います。

お礼

ご丁寧な回答ありがとうございます。この問題は実は曲線座標系を作成する場合、直交性を満たすようにしたいということから始まっています。ご解説によりますと、fを決めてgを変動させるということですが、f,gとも変動できます。ｆを借り決めしてgを変動させ、さらにｆをちょこっと変動してさらにgを変動させるというようなことをして直交性の最もいいところを探す、というようなアルゴリズムということでしょうか。
私は以下の様に考えて手詰まりになっています。SOR法による収束計算ですが、少しそれを作り変えて対応できないかと思っています。
1. f, gの境界値のデータ設定（境界値は固定）
2. f, gの領域内部点の初期値を有利なもの（按分とか）で決める
3. f, gを別々に収束させるのではなく、同じループで収束させる
4. 収束ループの中で直交性を評価する。直交だとゼロなので非ゼロ値
5. その値を参照してf, gの値を直交するように移動
6. 収束計算を続ける。ステップ3へ

このようにしたら収束計算をしながら直交するように、直交するように変化していくのではないでしょうか。しかし肝心の５のステップが皆目見当がつかないのですが。
あるいは、
Δf+λ1(∇f・∇g)＝0
Δg+λ2(∇f・∇g)＝0
λ1, λ2は任意の定数として...というような未定定数法を使うとかですが。
どうでしょうか。