ラプラス方程式

ラプラスの方程式を差分方程式で近似して反復法で解きます。
Fortran95 で導入された forall 代入文の、マスク機能を使って境界値を与えます。これによりラプラスの方程式を解く計算部分と境界条件を分離できます。

二次元ラプラス方程式

熱伝導や静電ポテンシャルなどはラプラス方程式を用いて記述できます。素朴な解き方としてある位置の値を隣接位置での値の平均値で置き換えるという反復法が知られています。
$$
f_{n+1}(x,y) = (f_n(x-h,y) + f_n(x+h,y) + f_n(x,y-h) + f_n(x,y+h)) \ /\ 4h^2
$$
ここでは二次元のラプラス方程式を適当な境界条件の下で反復法で解くことを考えます。

熱伝導の問題として、円形の空洞をもつパイプの中心からずれた位置にもう一本円形のパイプがあるとします。パイプのサイズ、温度、熱伝導率等は適当に規格化されているとします。外側のパイプの温度は100度に、内側のパイプは0度にそれぞれ保たれているとして、そのパイプの間の媒質の熱分布を求めます。媒質の温度の初期値は20度とします。

値が更新される媒質部分はマスクで指定します。マスクで境界を与えることにより、矩形の配列に対する計算ルーチンを用意するだけでよいことになります。下のプログラムでは１行で記述できています。

実行結果　

収束してゆく様子を簡単なキャラクターアニメーションによって表示することにします。コンソールの画面クリアのために、call execute_command_line('cls') を用います。execute_command_line は Fortran2008 で導入された命令で 'cls'は DOS の画面クリアコマンドです。温度の10の桁を書き出しています。100度を超えた部分は * になっています。

収束は遅く収束まで 8000 反復位必要なようです。

プログラムリスト

    program Heat
      implicit none
      integer, parameter :: n = 101
      real    :: mesh(0:n + 1, 0:n + 1)
      logical :: mask(n, n)
      integer :: i, j, iter
  ! define shape : eccentric pipe   
      mask = .false.       
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2) mask(i, j) = .true.
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 71)**2 < 10**2) mask(i, j) = .false.
  ! initial value
      mesh = 100.0
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2) mesh(i, j) = 20.0    
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 71)**2 < 10**2) mesh(i, j) =  0.0
  ! solve Laplace's equation by iterative method
      do iter = 0, 10000
        if (mod(iter, 100) == 0) then        ! character animation
          call execute_command_line('cls')  
          print '(2g/, (52i1))', ' iteration =', iter, int(mesh(::2,::4) / 10.0)
        end if  
        forall (i = 1:n, j = 1:n, mask(i, j)) mesh(i, j) = (mesh(i - 1, j) + mesh(i + 1, j) + mesh(i, j - 1) + mesh(i, j + 1)) / 4
      end do 
    end program Heat

参考：その他の形状

同心円の場合

  ! define shape : concentric circle  
      mask = .false.       
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2) mask(i, j) = .true.
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 10**2) mask(i, j) = .false.
  ! initial value
      mesh = 100.0
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2) mesh(i, j) = 20.0    
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 10**2) mesh(i, j) =  0.0

 iteration =       10000
****************************************************
*****************999999999999999999*****************
*************99999999999999999999999999*************
**********99999998888888888888888889999999**********
********999999888888887777777788888888999999********
******9999988888877777777777777777788888899999******
*****999998888777776666666666666677777888899999*****
****99998888777766665555555555556666777788889999****
***9999888877766655554444444444555566677788889999***
**999988887776655544433322223334445556677788889999**
**999888877766555443221111111122344555667778888999**
*99998887776665543321100000000112334556667778889999*
*99998887776655443221000000000012234455667778889999*
*99998887776655443211000000000011234455667778889999*
*99998887776655443321000000000012334455667778889999*
*99998888776665544332110000001123344556667788889999*
**999988877766655443332221122233344556667778889999**
**999988887776665554443333333344455566677788889999**
***9999888877776665555544444455555666777788889999***
****99999888877776666665555556666667777888899999****
******9999888887777776666666666777777888889999******
*******99999988888877777777777777888888999999*******
*********9999998888888888888888888888999999*********
************9999999988888888888899999999************
***************9999999999999999999999***************
*********************9999999999*********************

中心のずれた円の場合

  ! define shape : eccentric circle   
      mask = .false.       
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2) mask(i, j) = .true.
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 71)**2 < 10**2) mask(i, j) = .false.
  ! initial value
      mesh = 100.0
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2) mesh(i, j) = 20.0    
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 71)**2 < 10**2) mesh(i, j) =  0.0

 iteration =       10000
****************************************************
*****************999999999999999999*****************
*************99999999999999999999999999*************
**********99999999999999999999999999999999**********
********999999999999999999999999999999999999********
******9999999999988888888888888888899999999999******
*****999999999888888888888888888888888999999999*****
****99999998888888888877777777888888888889999999****
***9999999888888877777777777777777788888889999999***
**999999888888877777777766667777777778888888999999**
**999998888887777766666666666666667777788888899999**
*99999988887777766666555555555566666777778888999999*
*99999888877776666555554444445555566667777888899999*
*99999888877766655544444333344444555666777888899999*
*99998888777666555443332222223334455566677788889999*
*99998888777665544332211100111223344556677788889999*
**999988877766554432110000000011234455667778889999**
**999988877766554332100000000001233455667778889999**
***9999888776655433210000000000123345566778889999***
****99998887766554321000000000012345566778889999****
******9998887766554322100000012234556677888999******
*******99988877766554333222233345566777888999*******
*********9998887776665554444555666777888999*********
************9998887777666666667777888999************
***************9999888888888888889999***************
*********************9999999999*********************

内側のパイプが矩形の場合

  ! define shape : circle and square   
      mask = .false.       
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2      ) mask(i, j) = .true.
      forall (i = 1:n, j = 1:n, abs(i - 71) < 15 .and. abs(j - 71) < 15) mask(i, j) = .false.
  ! initial value
      mesh = 100.0
      forall (i = 1:n, j = 1:n, (i - 51)**2 + (j - 51)**2 < 50**2      ) mesh(i, j) = 20.0    
      forall (i = 1:n, j = 1:n, abs(i - 71) < 15 .and. abs(j - 71) < 15) mesh(i, j) =  0.0

 iteration =       10000
****************************************************
*****************999999999999999999*****************
*************99999999999999999999999999*************
**********99999999999999999999999999999999**********
********999999999999988888888888888899999999********
******9999999999988888888888888888888888999999******
*****999999999888888888887777777777788888889999*****
****99999999888888887777777777777777777788888999****
***9999999988888887777777666666666666677777888999***
**999999988888887777776666666555555666666777788999**
**999999988888777776666655555555555555555666778899**
*99999998888877777666655554444444444444455566778899*
*99999988888877776665555444333333223333334455667889*
*99999988888777766655544433222111111111222334567789*
*99999988888777766655444332110000000000000013456789*
*99999988888777666555443321100000000000000002356789*
**999998888877766655544332110000000000000000235679**
**999999888877776655544322110000000000000000245789**
***9999998888777666554433211000000000000000024689***
****99999988887776655543321100000000000000013579****
******9999988887776655543321000000000000000147******
*******99999888877766655433200000000000000038*******
*********9999988887776665543322222222334468*********
************9999988887776665555555566678************
***************9999998888877777778889***************
*********************9999999999*********************