２次元電位分布の数値計算 の履歴(No.2)

電位決定の基礎方程式†

物質中におけるマクスウェル方程式において

　&math( \begin{cases} \displaystyle\bm\nabla\times\bm E+\frac{\PD\bm B}{\PD t}=\bm o\\ \displaystyle\bm\nabla\cdot\bm B=0\\ \displaystyle\bm\nabla\times\bm H-\frac{\PD\bm D}{\PD t}=\bm i\\ \displaystyle\bm\nabla\cdot\bm D=\rho \end{cases} );

を仮定すれば、

　&math( \begin{cases} \bm\nabla\times\bm E=\bm o\\ \displaystyle -\frac{\PD\bm D}{\PD t}=\bm i\\ \bm\nabla\cdot\bm D=\rho \end{cases} );　　　　(*)

第１式より電位の存在が導かれ、

とともに第３式に代入すれば、

　&math( \bm\nabla\cdot\big(\varepsilon\bm\nabla\phi)=-\rho );

特に、真電荷のない領域では

　&math( \bm\nabla\cdot\big(\varepsilon\bm\nabla\phi)=0 );

これが基礎方程式となる。

電場の z 成分がゼロのとき 、成分に分けて書けば

　&math(\frac{\PD}{\PD x}\left(\varepsilon_x\frac{\PD}{\PD x}\phi\right)+

     \frac{\PD}{\PD y}\left(\varepsilon_y\frac{\PD}{\PD y}\phi\right)=\rho); 

媒質が等方的 ( ) かつ 均質 ( ) であれば、

　&math(\frac{\PD^2}{\PD x^2}\phi+

     \frac{\PD^2}{\PD y^2}\phi=\triangle\phi=\rho/\varepsilon); 

いわゆるポアソン方程式（ ならラプラス方程式）に帰着する。

シート状伝導体における基礎方程式†

(*) 式の第２、第３式から電荷保存則

　

が導かれる。

定常状態では より、

　

ここにシート状伝導体におけるオームの法則

電場と電位の関係式

を代入すると、

　

これがシート状伝導体における電位の基礎方程式となる。

成分に分けて書けば、

　&math(\frac{\PD}{\PD x}\left(\sigma_x\frac{\PD}{\PD x}\phi\right)+

       \frac{\PD}{\PD y}\left(\sigma_x\frac{\PD}{\PD y}\phi\right)=0); 

こちらも、コンダクタンスが等方的かつ一様であればポアソン方程式になる。

中心軸対称な系における円柱座標の基礎方程式†

基礎方程式

　

を円柱座標

　&math( \begin{cases} x=r\cos\theta\\ y=r\sin\theta\\ z=z\\ \end{cases} );

で書き表すと、

　&math( \begin{cases} \vspace{2mm}\displaystyle(\bm\nabla\phi)_r=\frac{\PD}{\PD r}\phi\\ \vspace{2mm}\displaystyle(\bm\nabla\phi)_\theta=\frac{1}{r}\frac{\PD}{\PD \theta}\phi\\ \displaystyle(\bm\nabla\phi)_z=\frac{\PD}{\PD z}\phi\\ \end{cases} );

　&math( \bm\nabla\cdot\bm v=\frac{1}{r}\frac{\PD}{\PD r}(rv_r)+\frac{1}{r}\frac{\PD}{\PD\theta}v_\theta+\frac{\PD}{\PD z}v_z );

より、

　&math( \bm\nabla\cdot(\varepsilon\bm\nabla\phi) =\frac{1}{r}\frac{\PD}{\PD r}\left(r\varepsilon_x\frac{\PD}{\PD r}\right)\phi

1. \frac{\PD}{\PD z}\left(\varepsilon_y\frac{\PD}{\PD z}\phi\right) =-\rho );

この場合には第1項の中に現れる r のため、単純なポアソン方程式にはならず、 方程式は座標 r に依存する

例：一様な抵抗を持つ２次元伝導体への電流注入†

• 計算領域は ｘ, ｙ ともに±１の範囲 ＝ 一辺２の正方領域
• (±0.5, 0.0) を中心に、半径0.05の円形領域から電流を注入した
  • 電位を±１とした

pdesolver2d.h†

LANG:cpp
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <float.h>
#include <algorithm>
#include <time.h>
#include <signal.h>

#ifndef __PDESOLVER2D
#define __PDESOLVER2D

// ２Ｄ正方領域で偏微分方程式 (PDE: Partial Differential Equation) を
// 解くためのクラス
//
// 領域サイズは n ピクセル x n ピクセル
// pixel に NaN を書き込むと、そこは自由端になる
// mask に非ゼロを書き込むと、そこは固定端になる
// 領域の周囲は必ず自由端もしくは固定端にしておかなければならない
//
// PDESolver2D 自体は抽象クラスなので、子クラスにて calc_new_value 
// を実装して、周囲４点の値から新しい値を求める漸化式を与えないと
// 実体化できない
//
class PDESolver2D
{
protected:
  int n;
  double *pixels;             // NaN を代入すると自由端になる
  int *masks;                 // 非ゼロを代入すると固定端になる
  int interrupted;            // 強制終了する
  double speed;               // 加速、減速定数 (speed > 1 で加速)
public:
  PDESolver2D(int n);         // n x n の正方領域を確保する
  virtual ~PDESolver2D();
  int size()      { return n; }
  double pixel(int p, int q)
                  { return pixels[n*q+p]; }
  void pixel(int p, int q, double v)
                  { pixels[n*q+p] = v; }
  int mask(int p, int q)
                  { return masks[n*q+p]; }
  void mask(int p, int q, int v)
                  { masks[n*q+p] = v; }
  void clear()    { memset(pixels, 0, n*n*sizeof(double));
                    memset(masks, 0, n*n*sizeof(int)); }
  void interrupt(){ interrupted = 1; }
  void set_speed(double s) { speed = s; }
  
  // 全ピクセルで更新差分の相対値が delta 未満になるまで繰り返し計算する
  // 更新差分の相対値とは更新差分をピクセル値で割った値のこと
  // 計算が終了すれば非ゼロを返す
  // 途中でキャンセルされればゼロを返す
  // log がゼロでなければ stderr に表示する進行状況を log にも出力する
  int solve(double delta, FILE *log = NULL);

  // 指定されたファイルにデータを書き込む
  void dump(const char *file_name);

  // 現在の計算値を引き継いで２倍の解像度の領域を確保する
  virtual PDESolver2D *create_double(int delete_this = 0);

protected:  // 子クラスで実装すべき virtual 関数

  // 自身の値 (c) と上下左右 (t, b, l, r) と座標値 (p,q) より、新しい値を計算して返す
  virtual double calc_new_value(double c, double t, double r, double b, double l, int p, int q) = 0;

  // サイズ n のインスタンスを作成して返す
  virtual PDESolver2D *create_instance(int n) = 0;

protected:  // 内部で用いるヘルパー関数

  // create_double 内から呼ばれる。２倍領域のピクセル座標を
  // 元画像に与えると、元画像のピクセル値を補完した値を返す
  double interpolate(int p, int q);

  // 自由端処理のためのヘルパー関数
  double replace_nan(double value, double replace)
                  { return isnan(value) ? replace : value; }
  double replace_nan(int p, int q, double replace)
                  { return replace_nan(pixel(p, q), replace); }

  // 進行状況を表示する
  void log_step(int i, double diff, time_t started_at, FILE *log);

  // 奇数ピクセルあるいは偶数ピクセルだけを計算する
  // ピクセル値の更新差分の相対値の最大値を返す
  double half_step(int even_or_odd);

  // 計算を１回分進める
  // 内部では half_step を２回呼ぶ
  // ピクセル値の更新差分の相対値の最大値を返す
  double step();
};
#endif

pdesolver2d.cpp†

LANG:cpp
#include "pdesolver2d.h"

PDESolver2D::PDESolver2D(int n)
{
  speed = 1;
  this->n = n;
  interrupted = 0;
  pixels = new double[n*n];
  masks = new int[n*n];
  clear();
}

PDESolver2D::~PDESolver2D()
{
  delete pixels;
  delete masks;
}

int PDESolver2D::solve(double delta, FILE *log)
{
  time_t started_at = time(NULL);
  
  int i;
  double diff = -1;
  int how_often = 53;
  for(i = 1; !interrupted; i++) {
    diff = step();
    if(log && (i % how_often == 0))
      log_step(i, diff, started_at, log);
    if(diff < delta) break;
  }
  if(log) log_step(i, diff, started_at, log);

  return !interrupted;
}

void PDESolver2D::log_step(int i, double diff, time_t started_at, FILE *log)
{
  if(log!=stderr)
    fprintf(log, "step (%6d, %e, %5.1fs)\n", i, diff, difftime(time(NULL), started_at));
  fprintf(stderr, "step (%6d, %e, %5.1fs)\r", i, diff, difftime(time(NULL), started_at));
}

double PDESolver2D::half_step(int even_or_odd)
{
  int p, q, nm1 = n-1;
  double result = 0;
  for(q = 1; q < nm1; q += 1 ){
    for(p = 1 + (even_or_odd  ^ (q % 2)); p < nm1; p += 2 ){
      if(mask(p, q)) continue;    // 固定端
      double value = pixel(p, q);
      if(isnan(value)) continue;  // 自由端
      double new_value = calc_new_value(
                  value,
                  replace_nan(p  , q-1, value),
                  replace_nan(p+1, q  , value),
                  replace_nan(p  , q+1, value),
                  replace_nan(p-1, q  , value),
                  p, q
              );
      pixel(p, q, value + speed * (new_value-value));
      
      double diff = fabs(new_value-value)/(fabs(new_value)+fabs(value));
      if( !isnan(diff) && diff > result ) result = diff;
    }
  }
  return result;
}

double PDESolver2D::step()
{
  return std::max( half_step(0), half_step(1) );
}

PDESolver2D *PDESolver2D::create_double(int delete_this)
{
  int p, q;
  PDESolver2D *result = create_instance(2*(n-1)+1);
  for(p = 0; p < result->size(); p++){
    for(q = 0; q < result->size(); q++){
      result->pixel(p, q, interpolate(p, q));
    }
  }
  if(delete_this) delete this;
  return result;
}

double PDESolver2D::interpolate(int p, int q)
{
  double i = 0;
  double result = 0, pix;
  int p2 = p/2, q2 = q/2;
  if(p % 2){  // need interpolation for p
    if(q % 2){  // need interpolation for q
      if (!isnan(pix = pixel(p2  , q2  )) ) { result += pix; i++; }
      if (!isnan(pix = pixel(p2+1, q2  )) ) { result += pix; i++; }
      if (!isnan(pix = pixel(p2  , q2+1)) ) { result += pix; i++; }
      if (!isnan(pix = pixel(p2+1, q2+1)) ) { result += pix; i++; }
      return i > 0 ? result / i : NAN;
    } else {    // not need interpolation for q
      if (!isnan(pix = pixel(p2  , q2  )) ) { result += pix; i++; }
      if (!isnan(pix = pixel(p2+1, q2  )) ) { result += pix; i++; }
      return i > 0 ? result / i : NAN;
    }
  } else {    // not need interpolation for p
    if(q % 2){  // need interpolation for q
      if (!isnan(pix = pixel(p2  , q2  )) ) { result += pix; i++; }
      if (!isnan(pix = pixel(p2  , q2+1)) ) { result += pix; i++; }
      return i > 0 ? result / i : NAN;
    } else {    // not need interpolation for q
      return pixel(p2, q2);
    }
  }
}

void PDESolver2D::dump(const char *file_name)
{
  FILE *f = fopen(file_name, "w");
  fwrite(pixels, sizeof(double), n*n, f);
  fclose(f);
}

main.cpp†

LANG:cpp
#include "pdesolver2d.h"

class SheetConductorSolver : public PDESolver2D
{
public:
  SheetConductorSolver(int n)
      : PDESolver2D(n)
      { }
  
  PDESolver2D *create_instance(int n) 
      { return new SheetConductorSolver(n); }
  
  double calc_new_value(double c, double t, double r, double b, double l, int p, int q);
};

double SheetConductorSolver::calc_new_value(double c, double t, double r, double b, double l, int p, int q)
{
  return 0.25 * (t + r + b + l);  // Laplace
}

////////////////////////////////////////////////////////

double sqr(double v) { return v*v; }

void set_boundary(PDESolver2D *solver)
{
  for(int p = 0; p < solver->size(); p++){
    for(int q = 0; q < solver->size(); q++){

      // outer boundary
      if( p==0 || q==0 || p==solver->size()-1 || q==solver->size()-1 )
        solver->pixel(p, q, NAN); // free edge

      // current injection
      double x = 2.0 * p / ( solver->size() - 1 ) - 1.0;
      double y = 2.0 * q / ( solver->size() - 1 ) - 1.0;

      // circular electrodes
      if ( sqr(x+0.5) + sqr(y) < sqr(0.05) ) {
        solver->mask(p, q, 1);    // fixed edge
        solver->pixel(p, q, -1.0);
      }
      if ( sqr(x-0.5) + sqr(y) < sqr(0.05) ) {
        solver->mask(p, q, 1);    // fixed edge
        solver->pixel(p, q, +1.0);
      }
    }
  }
}

PDESolver2D *solver = NULL;

// ^C が押されたら現在の計算を中断する
void signal_received(int sig)
{
  if(solver)
    solver->interrupt();
}

int main(int argc, const char *argv[])
{
  int initial_size = 17;
  
  solver = new SheetConductorSolver(initial_size);
  
  // ^C が押されたら signal_received へ飛ぶ
  signal(SIGINT, signal_received);
  
  // 徐々に解像度を上げながら計算する
  char fname[256];
  double delta = 1e-6;
  FILE *log;
  //  0   1   2    3    4    5     6
  // 17, 33, 65, 129, 257, 513, 1025
  for(int i = 0; i < 7; i++) {
    solver->set_speed(1);

    set_boundary(solver);
    sprintf(fname, "solver%dbefore.data", i);
    solver->dump(fname);

    sprintf(fname, "solver%d.log", i);
    log = fopen(fname, "w");
    int success = solver->solve(delta, log);
    fclose(log);

    if(success){
      sprintf(fname, "solver%dafter.data", i);
    } else {
      sprintf(fname, "solver%dinterrupted.data", i);
    }

    solver->dump(fname);

    fprintf(stderr, "\n*********** writing %s \n", fname);

    int delete_original = 1;
    solver = solver->create_double(delete_original);
  }
}

コンパイル＆実行†

LANG:console
$ ls
 main.cpp        pdesolver2d.cpp pdesolver2d.h
$ g++ -O4 main.cpp pdesolver2d.cpp
$ ls
  *a.out         main.cpp        pdesolver2d.cpp pdesolver2d.h
$ ./a.out
 step (   337, 9.730739e-07,   0.0s)
 *********** writing solver0after.data
 step (   948, 9.985805e-07,   0.0s)
 *********** writing solver1after.data
 step (  3011, 9.996066e-07,   1.0s)
 *********** writing solver2after.data
 step (  5274, 9.995926e-07,   1.0s)
 *********** writing solver3after.data
 step (   464, 9.984991e-07,   1.0s)
 *********** writing solver4after.data
 step (   958, 9.997507e-07,   5.0s)
 *********** writing solver5after.data
 step (  1402, 9.997880e-07,  29.0s)
 *********** writing solver6after.data
$

この系は素性がいいので、speed = 1.5 にするとさらに 1.5 倍速くなるけれど、 speed = 2 では発散する。

計算結果ファイル†

solver?after.data

のような名前のファイルが格段の計算結果となる。

LANG:cpp
double pixels[n*n];

の内容を単にダンプしたものなので、sqrt( (ファイルサイズ)/8バイト ) で n を求めて、 n*n の二次元配列として double 値を読み込み、グラフソフトで表示することになる。

Igor Pro なら、

LANG:igor
function show(filename)
string filename;
	GBLoadWave/B/T={4,4}/W=1/O filename
	variable n = sqrt(numpnts(wave0))
	make/O/N=(n,n) image
	wave wave0, image
	image = wave0[p+q*n]
	setscale x,0,2,image
	setscale y,-1,1,image
end

を使って、

LANG:command
*show("z:/electric-field2/solver6after.data")
*newimage image
*ModifyImage image ctab= {*,*,Rainbow,0}

とすればカラー表示できる。

等高線†

[グラフ]-[グラフに追加]-[等高線プロット]

から、

[アピアランス] の

• [ラインカラー] を [すべて同じカラー]-[白]
• [ラベル詳細] で、
  • [背景] を [透明]
  • [ラベルテキスト]-[カラー]-[黒]

とすれば、

LANG:console
*AppendMatrixContour/T image
*ModifyContour image rgbLines=(65535,65535,65535),labelBkg=1,labelRGB=(0,0,0)

と同等の動作になり、

が得られる。

電場強度を見るには†

強度のみであれば をプロットすればいい。

Igor Pro なら、

LANG:console
*duplicate image, image2
*image2 = sqrt((p == 0 ? image[p+1][q]-image[p][q] : image[p][q]-image[p-1][q])^2+(q == 0 ? image[p][q+1]-image[p][q] : image[p][q]-image[p][q-1])^2)
*image2 = image2 == 0 ? NaN : image2

電場の方向を見るには†

function calcfield()
  wave image
  variable n = dimsize(image,0)
  variable d =  20
  variable pp, qq, nn = trunc((n-1)/d)
  variable ma = 0
  make/O/N=(nn*nn) posx, posy
  make/O/N=(nn*nn,2) field
  for(pp=0; pp<nn-1; pp+=1)
    for(qq=0; qq<nn-1; qq+=1)
      posx[pp+qq*nn] = DimOffset(image, 0) + DimDelta(image, 0) * pp * d
      posy[pp+qq*nn] = DimOffset(image, 1) + DimDelta(image, 1) * qq * d
      variable fieldx = -image[pp*d+1][qq*d   ]+image[pp*d][qq*d]
      variable fieldy = image[pp*d   ][qq*d+1]-image[pp*d][qq*d]
      variable fieldr = sqrt(fieldx^2+fieldy^2)
      if (fieldr > ma) 
      	ma = fieldr
      endif
      field[pp+qq*nn][1] = atan2(fieldy, fieldx)
    endfor
  endfor
  field[][0] = 40 // 長さ
end

として、

LANG:console
*calcfield()
*appendtograph/t/l posy vs posx
*ModifyGraph mode=3,arrowMarker(posy)={field,0.5,0,0.5,1},msize=0.1
*ModifyGraph rgb(posy)=(0,0,0)

などとすれば、下記のようなグラフを得る。

平面での漸化式†

停留条件は、

　&math( \frac{\PD U}{\PD \phi_{p,q}} &=\frac{\varepsilon_{p-1/2,q-1/2}}{2}\big(2\phi_{p,q}-\phi_{p-1,q}-\phi_{p,q-1}\big)\\ &+\frac{\varepsilon_{p+1/2,q-1/2}}{2}\big(2\phi_{p,q}-\phi_{p+1,q}-\phi_{p,q-1}\big)\\ &+\frac{\varepsilon_{p-1/2,q+1/2}}{2}\big(2\phi_{p,q}-\phi_{p-1,q}-\phi_{p,q+1}\big)\\ &+\frac{\varepsilon_{p+1/2,q+1/2}}{2}\big(2\phi_{p,q}-\phi_{p+1,q}-\phi_{p,q+1}\big)\\ &=0 );

これを整理した漸化式は、

　&math( &\big[2(\varepsilon_{p-1/2,q-1/2}+\varepsilon_{p+1/2,q-1/2}+\varepsilon_{p-1/2,q+1/2}+\varepsilon_{p+1/2,q+1/2})\big]\phi_{p,q}\\ &=\varepsilon_{p-1/2,q-1/2}\big(\phi_{p-1,q}+\phi_{p,q-1}\big)+\varepsilon_{p+1/2,q-1/2}\big(\phi_{p+1,q}+\phi_{p,q-1}\big)\\ &\ +\varepsilon_{p-1/2,q+1/2}\big(\phi_{p-1,q}+\phi_{p,q+1}\big)+\varepsilon_{p+1/2,q+1/2}\big(\phi_{p+1,q}+\phi_{p,q+1}\big) );

となる。

これも上下左右の格子点の値しか使っていないため、加えて epsilon(x,y) を与えてやることで 上記のソルバーを使って解ける。

円柱座標での漸化式†

円柱座標の場合、対応する停留条件は、

　&math( \frac{\PD U[\phi]}{\PD \phi_{p,q}}&=\frac{\pi \Delta}{2}\Big\{ p\varepsilon_{p+1/2,q+1/2}\big(-\phi_{p+1,q}+\phi_{p,q}-\phi_{p,q+1}+\phi_{p,q}\big)+\\ &\hspace{11mm}\ (p-1)\varepsilon_{p-1/2,q+1/2}\big(\phi_{p,q}-\phi_{p-1,q}-\phi_{p,q+1}+\phi_{p,q}\big)+\\ &\hspace{11mm}\ p\varepsilon_{p+1/2,q-1/2}\big(-\phi_{p+1,q}+\phi_{p,q}+\phi_{p,q}-\phi_{p,q-1}\big)+\\ &\hspace{11mm}\ (p-1)\varepsilon_{p-1/2,q-1/2}\big(\phi_{p,q}-\phi_{p-1,q}+\phi_{p,q}-\phi_{p,q-1}\big)\Big\}\\ &=0 );

整理すると、

　&math( &2\big[p\varepsilon_{p+1/2,q+1/2}+(p-1)\varepsilon_{p-1/2,q+1/2}+p\varepsilon_{p+1/2,q-1/2}+(p-1)\varepsilon_{p-1/2,q-1/2}\big]\phi_{p,q}\\ &=p\varepsilon_{p+1/2,q+1/2}\big(\phi_{p+1,q}+\phi_{p,q+1}\big)+ (p-1)\varepsilon_{p-1/2,q+1/2}\big(\phi_{p-1,q}+\phi_{p,q+1}\big)+\\ &\hspace{5mm}p\varepsilon_{p+1/2,q-1/2}\big(\phi_{p+1,q}+\phi_{p,q-1}\big)+ (p-1)\varepsilon_{p-1/2,q-1/2}\big(\phi_{p-1,q}+\phi_{p,q-1}\big)\\ );

となる。