http://www.geocities.jp/laprog321/

質問

http://www.geocities.jp/laprog321/

この上記URLのプログラムについて質問です。

このプログラムはラプラス方程式を差分法で解くプログラムと書いていたので、
自分はまずプログラムを理解しようとするために、
インターネットなどで、ラプラス方程式の差分法による解き方についていろいろ調べて、式を考えたりしました。

そして、このプログラムがいう差分法とは中央差分方程式を用いて解いているのかなと思ったのですが、
何卒知識が少ないもので、実際に合ってるかどうかが分かりません。
また、ラプラス方程式をコンピュータに解かせる場合、解き方の反復法として、
いろいろな種類（ガウスザイデル法、ヤコビ法、SOR法など）の方法があることがわかりました。

ここで質問させていただきます。

プログラムでは、実際にどの方法を利用してどのように計算しているのか教えてほしいです。

例えば、ここの部分（行数など）で「～式」、「～法」を用いてどのように計算してるかなど教えてほしいです。

プログラムの流れについては、違う質問でご教授いただいたので、上記した質問のように
プログラム本文を使ってではなく、実際にこのプログラムでは、ラプラス方程式をどのような式を用いて、
また、どのような方法で解いているかを、矛盾しているようですが、プログラムと対比して教えていただきたいです。

同じような質問を繰り返し失礼いたします。
文章力のない質問で申し訳ありませんがよろしくお願いします。

回答No.2

麻野 なぎ（@AsanoNagi） · ベストアンサー

一般に、Ux（U の x による微分）を、『中心差分近似』すると、
Ux = (u(x+h / 2, y)-u(x-h / 2, y)) / h になり、これを用いて、
Uxx（U の x による２階偏微分）を『中心差分近似』すると、
Uxx = (Ux(x+h / 2, y)-Ux(x-h / 2, y)) / h
= (u(x+h / 2+h / 2 , y) / h-u(x+h / 2-h / 2, y) / h) / h （前半分）
- (u(x-h / 2+h / 2 , y) / h-u(x-h / 2-h / 2, y) / h) / h （後半分）

（前半分） = (u(x+h, y)-u(x, y)) / h^2
（後半分） = (u(x, y)-u(x-h, y)) / h^2
つまり、
Uxx = (u(x+h, y)-u(x, y)) / h^2-(u(x, y)-u(x-h, y)) / h^2
= (u(x+h, y)-2u(x, y)+u(x-h, y)) / h^2 となります。
同様に、
Uyy = (u(x, y+h)-2u(x, y)+u(x, y-h)) / h^2
これを、ラプラス方程式 Uxx+Uyy = 0 に代入して整理すると、

u(x+h, y)-2u(x, y)+u(x-h, y)+(u(x, y+h)-2u(x, y)+u(x, y-h) = 0
u(x+h, y)+u(x-h, y)+u(x, y+h)+u(x, y-h)-4u(x, y) = 0
u(x, y) = (u(x+h, y)+u(x-h, y)+u(x, y+h)+u(x, y-h))/4

h = 1 として、（さらに、x, y を i, j にして）配列で書き直せば、

u[i][j]=(u[i-1][j]+u[i+1][j]+u[i][j-1]+u[i][j+1])/4
各点におけるこの関係を満たすように、u を決定すればいいのですが、この連立方程式を解くのは手間なので、
u(new)[i][j]=(u(old)[i-1][j]+u(old)[i+1][j]+u(old)[i][j-1]+u(old)[i][j+1])/4
の計算を反復させることで近似値を得ます。

この反復計算は、すべての u(new) 算出に、u(old) を使っているので、ヤコビ法によるものになります。
プログラムの中も、これと、ほとんど同じ式なので、どこでやっているかは見当がつくかなと思います。

回答No.3

麻野 なぎ（@AsanoNagi）

あと、どう考えても、初期設定が変というか、不自然なんですね。

・ラプラス方程式の問題なのに、境界をすべて０にしている。
　これだと、計算するまでもなく、U(x, y) = 0 に収束するはず。

・それを別にしても、sin((j-1)/(SIZE*M_PI)) という式は変。
　普通、M_PI （円周率ですよね？）は、（割るのじゃなくて）掛けるはず
　もっとも、単純に、sin((j-1)/SIZE*M_P) にすると、int / int
　のほうが結合力が強いので、これはこれで、期待しない結果になる。sin(M_P*(j-1)/SIZE) あたりかな？

・そもそも、初期値を設定するのに、（i を使わずに） j だけを使っているのも
　ちょっと変な気がする（これは、「ちょっと」程度ですけど）

このプログラムの出所が気になるところではあります。

回答No.1

Tacosan

「ヤコビ法」とか「ガウス=ザイデル法」, 「SOR法」などがそれぞれどのようなものであるかは理解できていますか?
そして, そもそもあなたは何がしたいのですか?