ラプラス方程式を差分法で解いてみる

これまた異様な題目ですが(^^;、ラプラス方程式を差分法を使ってコンピュータに解かせてみようと思います。まぁ人生山だらけ、 渡る世間は鬼だらけ、友達とボウリングに行ったらガターだらけ(しかもバツゲーム付き^^;) （最終更新日：というわけで、何を血迷ったかこんな中途半端なページを作ってしまいました。今思えば作るんじゃなかったと思いますが、消すのも面倒だし残してます）。使用言語はＣ言語です。例によって私はこの道のプロではございませんので、至らぬ点はご容赦くださいm(_ _)m。

トップへ戻る

ラプラス方程式というのは、コレです。

ここでｕは x, y, z の関数です。Δは、ラプラシアンといって、次のように定義（≡）される演算子です。

式(1)を成分で表記すると、こうなります。

よってさらに詳しく書くと、こうなります。

この式は、いろいろな物理量の平衡状態を表す方程式で、電界や磁界の計算、温度分布の計算などに使われるそうです。この式の右辺がさらにｘ、ｙ、ｚの関数となっている場合、つまり

の形をしている微分方程式をポアソン方程式といいます。ラプラス方程式とポアソン方程式は、だ円型微分方程式の典型例なんだそうです。(←勉強しながら書いている^^;)

式(3)を２次元で書くとこうなります(以後(x,y,z)は省略します)。

まず、これをコンピュータに計算させてみます。これを行うには、偏微分記号 ∂/∂x、∂/∂y、∂/∂z の意味を知っている必要があります。たとえば、∂/∂x は、ｘ方向の微分、つまりものすごく細かい区間でのｘ方向の傾きを表しています。∂^2/∂x^2 は、 ∂/∂x したものをさらに ∂/∂x します。数学的には(と私のような者が言ったら数学者の方に怒られるかもしれませんが)、たとえば∂u/∂xは、次のような意味だったと思います。

ｈは微小な区間を表していて、ｈ→0はその微小な区間を限りなく小さくすることを表しています (しかしゼロではない)。今回は２階の偏微分を扱うので、式(7)をもう一度 ∂/∂xする必要があります（2010.1.11: 符号訂正）。

コンピュータに計算させる場合は、h をある程度小さな間隔で近似します。この間隔をそれぞれ Δｘ、Δｙ(注：これはデルタです。同じ記号を用いていますが、これはラプラシアンではありません。ｈと同様、微小な幅を意味します。)とすると、式(6)は、式(8)を用いて次のように近似できます。

しかし、コンピュータに計算させる場合は精度上の問題から式(8)を、

のように変更し、式(9)を

で近似します。式(9)は前進(前方)差分を用いた式、式(11)は中心差分を用いた式と言います。ちなみに後退(後方)差分を用いた式は次のようになります。

トップへ戻る

とりあえず中心差分の式(11)を解くプログラムを考えます。ラプラス方程式は、流体力学や電磁気学、伝熱学などで現れる式 (このうち私が現れることを確認しているのは流体力学と伝熱学) らしいですが、ここでは伝熱学を考えてみます。つまり、ラプラス方程式で扱う物理量を熱とします。そこで、便宜上式(11)を温度Ｔを用いて、

のように書き換えることにします。

いま、図に示すような、熱伝導率が一定の矩形形状をした平板の熱伝導を考えます。板の左と下をそれぞれT=1.0、T=0.5で加熱するものとし、板の右と上に一定温度T=0.0 を与えるものとします(境界条件)。また、板に熱を加える前は板上の全ての領域でT=0.0であるとします(初期条件)。

このように板に熱を加えると、最終的にどのような温度分布に落ち着くのかを (つまり熱平衡状態の温度分布を)、上のラプラス方程式を解くことによって知ることができそうです。

上の図のように、平板の右端が NX、上端が NY となるように等しく分割すると、 Δｘ=1/(NX-1)、Δｙ=1/(NY-1) になります。こうして作成した微小な要素を格子(Grid) と言います。ここでは、温度Tを線と線が交わる点で定義します。そして、T(i,j) という表記方法を用いることにします。ちなみに、T(i,j) とは上の図の赤の点における温度の値です。この表記を用いて、式(11)を書き直すと、

となります。コンピュータに計算させるためには、式(14)を T(i,j) について解き、

を得ます。さらに、Δｘ、Δｙは計算中は定数なので、式(15)をさらに、

と書き直すとスッキリして良いと思います。

式(16)は、４×４くらいの分割数で、T(i,j)を考えてみると (実際にi,jに値を順に入れてゆくと)、連立方程式になることが分かると思います。そこで、この連立方程式を反復解法で解いてやることにします。反復解法には、 ヤコビ法、ガウス・ザイデル法、ＳＯＲ法などがありますが、ここではガウス・ザイデル法を使います。

やや話が飛躍してしまい申し訳ないのですが、中心差分を使ったプログラムです。
[プログラムのソースコードをダウンロードする]
[2007/12/8]
プログラム中の const double dx = 1.0 / NX; を const double dx = 1.0 / (NX-1); に修正しました(dy も同様)。
まちがえてました(^_^;
都合により laplace.lzh → laplace.zip に変更、あと実行ファイルをアーカイブから削除しました。
コンパイルは gcc laplace.c -o laplace で laplace (UNIX系) または laplace.exe (Windows) ができます。
実行は ./laplace > result.txt (UNIX系) または laplace.exe > result.txt (Windows) です。

↓とりあえずプログラム見たい人用

     1|/******************************************************************************
     2|** [laplace.c]
     3|**   ラプラス方程式を差分法により解きます。
     4|**
     5|**   使い方： laplace.exe > result.txt
     6|******************************************************************************/
     7|#include <stdio.h>
     8|
     9|/*
    10|** 計算領域の分割数 NX, NY を決めます。理論上は分割数を細かくすれば
    11|** するほど計算誤差は減るはずですが、コンピュータが扱う浮動小数点数
    12|** は有限値なので、あまり細かくしすぎると、桁落ちによって逆に
    13|** 計算精度が落ちると思います。いろいろな DX, DY で結果を比較して、
    14|** 計算精度が充分に良いところで、かつ計算速度も妥当な分割数を選ぶと
    15|** 良いと思います。ここでは仮に 41 x 41 とします。
    16|** NX と NY は同じ数にした方が、微分方程式の特性に合うらしいです。
    17|*/
    18|
    19|#define NX 41
    20|#define NY 41
    21|
    22|
    23|/*
    24|** for (y のループ) { for (x のループ) { (x 方向の計算); } }
    25|** とする場合は、Ｃ言語の配列のメモリ配置を考慮して、T(x, y)
    26|** ではなく、T(y, x) とします。
    27|*/
    28|double T[NY][NX];
    29|
    30|
    31|/*-----------------------------------------------------------------------------
    32|** [SetBoundaryCondition]
    33|** 境界条件を設定します。
    34|** ２次元配列の扱い方はいろいろあります。書籍や web page 等を
    35|** あたってみてください。
    36|**---------------------------------------------------------------------------*/
    37|void SetBoundaryCondition(double T[NY][NX])
    38|{
    39|    int x;
    40|    int y;
    41|
    42|    /* 領域の底の境界値 0.5 */
    43|    for (x = 0; x < NX; x++)
    44|    {
    45|        T[0][x] = 0.5;
    46|    }
    47|
    48|    /* 領域の上の境界値は 0.0 */
    49|    for (x = 0; x < NX; x++)
    50|    {
    51|        T[NY-1][x] = 0.0;
    52|    }
    53|
    54|    /* 左の壁の境界値は 1.0 */
    55|    for (y = 0; y < NY; y++)
    56|    {
    57|        T[y][0] = 1.0;
    58|    }
    59|
    60|    /* 右の壁の境界値は 0.0 */
    61|    for (y = 0; y < NY; y++)
    62|    {
    63|        T[y][NX-1] = 0.0;
    64|    }
    65|}
    66|
    67|
    68|int main(void)
    69|{
    70|    const double dx = 1.0 / (NX-1);
    71|    const double dy = 1.0 / (NY-1);
    72|
    73|    /* 式(16)から、C1 と C2 の値を先に計算しておきます。 */
    74|    const double C1 = 0.5 * dx*dx / (dx*dx + dy*dy);
    75|    const double C2 = 0.5 * dy*dy / (dx*dx + dy*dy);
    76|
    77|    int i;
    78|    int j;
    79|
    80|    double errorMax;
    81|    double previousValue;
    82|
    83|    /*
    84|    ** 境界条件を設定します。この境界条件は
    85|    ** ディリクレ条件なので、はじめに１度だけ設定
    86|    ** すればＯＫです。
    87|    */
    88|    SetBoundaryCondition(T);
    89|
    90|    /*
    91|    ** 前回の計算値と今回の計算値との最大誤差が
    92|    ** 0.00001 を下回るまで繰り返します。
    93|    */
    94|    do
    95|    {
    96|        errorMax = 0.0;
    97|
    98|        /* 計算領域は境界の１つ内側です。 */
    99|        for (i = 1; i < NY-1; i++)
   100|        {
   101|            for (j = 1; j < NX-1; j++)
   102|            {
   103|                previousValue = T[i][j];
   104|
   105|                /* ガウス・ザイデル法 */
   106|                T[i][j] = C1*(T[i+1][j] + T[i-1][j]) + C2*(T[i][j+1] + T[i][j-1]);
   107|
   108|                if (errorMax < fabs(T[i][j] - previousValue))
   109|                {
   110|                    errorMax = fabs(T[i][j] - previousValue);
   111|                }
   112|            }
   113|        }
   114|    } while (errorMax > 0.00001);
   115|
   116|    /* 全領域の計算結果を出力します。 */
   117|    for (i = 0; i < NY; i++)
   118|    {
   119|        for (j = 0; j < NX; j++)
   120|        {
   121|            printf("%d\t%d\t%15.15f\n", j, i, T[i][j]);
   122|        }
   123|    }
   124|
   125|    return 0;
   126|}
   127|
   128|

なお、ＳＯＲ法を使った計算部分は、

    previousValue = T[i][j];

    /* ガウス・ザイデル法 */
    T[i][j] = C1*(T[i+1][j] + T[i-1][j]) + C2*(T[i][j+1] + T[i][j-1]);

の部分を

    previousValue = T[i][j];

    /* ＳＯＲ法 */
    T[i][j] += SOR_COEF*(C1*(T[i+1][j] + T[i-1][j]) + C2*(T[i][j+1] + T[i][j-1]) - previousValue);

とすればできます。ここで、SOR_COEFの実用的な値の範囲は 1 < SOR_COEF < 2 です。つまり、最新の値と前の値との差(最終値への変化幅)を少し大きくすることで、最終値に到達するのを早めてやるのがＳＯＲ法です。SOR_COEF = 1.17 とか、 1.2 といった値にすることが多いようですが、この値の根拠は詳しくは知りません。

トップへ戻る

このページの数式は次のソフトウェアで作成しています：