9 (応用例)ポアソン方程式 緩和法

9.1 概要

今度はポアソン方程式を緩和法(PointJacobi法)で解きます．

9.2 アルゴリズムの解説

計算の繰り返し数をnとし，n回目の値とn+1回目の値が1つに収束して行く時，

$$\phi _j^{n+1}=\frac{1}{2}\left(\phi _{j+1}^n+\phi _{j-1}^n-\Delta x^2\rho _j^n\right)$$

はポアソン方程式の差分式となります．
これを，境界条件 $\phi _1(0)=0$， $\phi _N(1)=0$を考慮して行列表現すると，次の様になります．

$$\left( \begin{array}{c} \phi _2^{n+1}\ \phi _3^{n+1}\ \vd... ...ho _j\ \vdots\ \rho _{N-2}\ \rho _{N-1} \end{array} \right) \right)$$

これを初期分布 $\bm{\phi ^0}$を適当に決め， $\parallel \bm{\phi ^{n+1}}-\bm{\phi ^n}\parallel \leq \epsilon$ となるまで繰り返し計算を行えば，分布$\phi$が求められます．
より詳しい解説は他を御覧下さい．

9.3 例題

section8と同じ問題を解く．


コードを以下に示します．

//====================================================================[include]
#include ``cpplapack.h''

//=======================================================================[main]
/*! main */
int main(int argc, char** argv)
{
  //// declare objects ////
  const int N(201);
  const double dx(1/(double)(N-1));
  const double pai(3.141592);
  const double k(2*pai);
  const double eps(1.0e-6);
  CPPL::dgbmatrix A(N-2,N-2,1,1);
  CPPL::dcovector rho(N-2);
  CPPL::dcovector phi(N-2),phi_old(N-2);

  //// make dgbmatrix A ////
  A.zero();
  for(int i=0; i<A.n-1; i++){
    A(i,i+1) = 1;
    A(i+1,i) = 1;
  }

  //// make dcovector rho ////
  for(int i=0; i<rho.l; i++){
    rho(i) = -k*k*sin((i+1)*k*dx);
  }
  rho *= dx*dx;

  //// solve ////
  do{
    phi_old = phi;
    phi = A*phi_old-rho;
    phi *= 0.5;
  }while(nrm2(phi-phi_old)>eps);

  //// print ////
  std::cout << ``x       phi(x)      解析解'' << std::endl;
  std::cout << ``0 0 0'' << std::endl;
  for(int i=0; i<rho.l; i++){
    std::cout << (i+1)*dx << `` `` << phi(i) << `` `` << sin((i+1)*k*dx) 
              << std::endl;
  }
  std::cout << ``1 0 0'' << std::endl;

}

実行結果は以下の様になります．

x       phi(x)      解析解
0 0 0
0.005 0.031407 0.0314108
0.01 0.062783 0.0627905
0.015 0.094097 0.0941083
0.02 0.125318 0.125333
0.025 0.156416 0.156434
0.03 0.187359 0.187381
0.035 0.218117 0.218143
0.04 0.24866 0.24869
0.045 0.278958 0.278991
0.05 0.30898 0.309017
0.055 0.338697 0.338738
0.06 0.36808 0.368124
0.065 0.3971 0.397148
0.07 0.425728 0.425779
0.075 0.453936 0.45399
0.08 0.481696 0.481754
0.085 0.50898 0.509041
0.09 0.535762 0.535827
0.095 0.562016 0.562083
0.1 0.587715 0.587785
0.105 0.612833 0.612907
0.11 0.637347 0.637424
0.115 0.661232 0.661312
0.12 0.684465 0.684547

$\vdots$