ラプラス方程式を差分法で解くプログラムについての質問

2010/01/21 16:42

このQ&Aのポイント

ラプラス方程式を差分法で解くプログラムについての質問です。具体的には、プログラムがどの方法を利用してどのように計算しているのか知りたいです。
ラプラス方程式の差分法による解き方について調べ、それをもとにプログラムを理解しようとしています。そこで、実際にこのプログラムでは、どのような式を用いて計算しているのか、そしてどの方法を利用しているのか教えてください。
ラプラス方程式を差分法で解くプログラムに関して質問があります。具体的には、このプログラムがどのような方法を利用してどのように計算しているのか、そしてどのような式を用いているのか知りたいです。

qwe1232
お礼率17% (7/41)

C・C++・C#
回答数3
ありがとう数1

みんなの回答 （3）
専門家の回答

質問者が選んだベストアンサー

ベストアンサー

麻野なぎ（@AsanoNagi）
ベストアンサー率45% (763/1670)

2010/01/21 17:57 回答No.2

一般に、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) を使っているので、ヤコビ法によるものになります。プログラムの中も、これと、ほとんど同じ式なので、どこでやっているかは見当がつくかなと思います。

その他の回答 (2)

麻野なぎ（@AsanoNagi）
ベストアンサー率45% (763/1670)

2010/01/22 06:57 回答No.3

あと、どう考えても、初期設定が変というか、不自然なんですね。・ラプラス方程式の問題なのに、境界をすべて０にしている。　これだと、計算するまでもなく、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 だけを使っているのも　ちょっと変な気がする（これは、「ちょっと」程度ですけど）このプログラムの出所が気になるところではあります。　