課題5プログラム例

課題5−１

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 100         // 格子数
#define VX0 5.0       // 左境界の電位
#define VX1 10        // 右境界の電位
#define VY 15.0      // 上境界の電位　（下はゼロ）
#define MIN_CHANGE  0.001  // 終了誤差条件

void outResult(double v[N+2][N+2]);
double v[N+2][N+2], v_tmp[N+2][N+2];   // 0とN+1が境界

int main() {
  int i,j, iteration=0;
  double change= 0.0,  change_tmp;
  // 初期値のセット
  for(i=1; i<=N; i++) {
    v[0][i] = VX0;
    v[N+1][i] = VX1;
    v[i][0] = VY;
    v[i][N+1] = VY;

    for(j=1; j<=N; j++) {
      v[i][j] = (VX0 + VX1 +VY*2 )*0.25;  // 内部格子上の初期値
    }
  }

  // 緩和法の繰り返し
  // v[][]を使って次の値v_tmp[][]を計算するとともに、差の一番大きかったものを記憶
  // 差の最大値がある値を下回ったら終了
  do {
    change = 0.0;
    for(i=1; i<= N; i++) {
      for(j=1; j<= N; j++) {
	v_tmp[i][j] = (v[i-1][j] + v[i+1][j] + v[i][j-1] + v[i][j+1])*0.25;
	change_tmp = fabs(v_tmp[i][j] - v[i][j])/v[i][j];
	//　printf("tmp=%f %f %f,", change_tmp, v[i][j], v_tmp[i][j]); //検証用
	if(change_tmp > change) change = change_tmp;
      }
    }
    for(i=1; i<= N; i++) {
      for(j=1; j<= N; j++) {
	v[i][j] = v_tmp[i][j];
      }
    }
    printf("%f interation=%i\n",change,iteration); // 変化の最大値のチェック出力
    iteration++;
  } while(change > MIN_CHANGE);

  outResult(v, N);
}

void outResult(double v[N+2][N+2]) {
  int i, j;
  // 結果出力
  FILE *FP;
  FP = fopen("e_field.csv","w");
  for(i=1; i<= n; i++) {
    for(j=1; j< n; j++) {
      fprintf(FP, "%f,", v[i][j]);
    }
    fprintf(FP,"%f\n", v[i][n]);
  }
  fclose(FP);
}

(注): main内で二次元配列を使っているとき、関数にその配列を渡す場合、配列の引数に要素数を書かざるを得ないのは格好がわるい。何個の場合でも対応可能にするには、main文のなかでもポインタを使い、ポインタを渡すようにする必要がある。

何人かのレポートにもあったが、

void outResult(double v[N+2][N+2]);

結果を直接gnuplotに渡して描画するプログラム例

上記のプログラムで出力関数outResultを次のように変えると直接gnuplotを起動し、図示できる。

あるプログラムから別のプログラムへデータを流し込むことをパイプ(pipe)と呼ぶ。下では、fopenの代わりにpopenを使っている。

void outResult(double v[][N+2], int n) {
  int i, j;
  // 結果出力
  FILE *FP;
  FP = popen("gnuplot -persist","w");
  fprintf(FP, "splot '-' matrix with pm3d\n");
  // splotにmatrixオプションをつけると、matrix型でz値をならべたデータに対応できる。
  for(i=1; i<= n; i++) {
    for(j=1; j<= n; j++) {
      fprintf(FP, "%f ", v[i][j]);
    }
    fprintf(FP,"\n");

  }
  pclose(FP);
}

残念ながら、教室のLinuxに載っているgnuplotは、パイプで渡した場合、三次元プロットをマウスで操作できない。 よってマウスでグリグリ動かしたい場合は、授業で説明したように、gnuplotを端末内で立ち上げ、ファイルを読み込んで描画する方法を用いてほしい。

なお、gnuplot-x11版を使うと、上記の問題は生じない。（手元のコンピュータではそちらを使っていたので、気が付かなかった。）

課題5−２

ヒント

内側の境界を繰り返し計算の範囲から除くのではなく、内側の境界は無視して計算し、新しい値を計算した後で、強制的に内側の境界上の値を境界の値にセットする方が楽です。（このような効率的なアルゴリズムを考えるのが大事です。）

つまり、緩和法の繰り返し計算のループの中で、

1. 外側境界内の全体のサイトを緩和法によって値を更新
2. 内側境界上の値を初期の値にセットし直す
3. 更新前後の変化の最大値を求める

という順序で実行するのが簡便でしょう。

課題のページヘ

/*
  長方形同軸円筒の電位分布
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NX 100         // 格子数（外側）
#define NY 100
#define MX  30         // 内側境界の大きさ
#define MY  30
#define V_IN 5.0       // 内側境界の電位
#define V_OUT 10.0      // 外側の電位

#define MIN_CHANGE  0.00001  // 終了誤差条件

void outResult(double v[][NY+2], int nx, int ny);
double v[NX+2][NY+2], v_tmp[NX+2][NY+2];   // 0とN+1が境界

int main() {
  int i,j, iteration=0;
  /* 内側境界の位置 */
  int inRight=(NX-MX)/2, inLeft=(NX+MX)/2, inBottom=(NY-MY)/2, inTop=(NY+MY)/2;
  double change= 0.0,  change_tmp;


  for(i=1; i<=NX; i++) v[i][0] = v[i][NY+1] = V_OUT;
  for(i=1; i<=NY; i++) v[0][i] = v[NX+1][i] = V_OUT;

  for(i=1; i<=NX; i++) {
    for(j=1; j<=NY; j++) {
      v[i][j] = (V_IN+V_OUT)/2.0;  // 内部格子上の初期値
    }
  }
  do {
    change = 0.0;
    // 内側境界は毎回設定（「内側境界は元に戻す」で行っているので必要ない）
    for(i=inRight; i<=inLeft; i++) v[i][inBottom] = v[i][inTop] = V_IN;
    for(i=inBottom; i<=inTop; i++) v[inRight][i] = v[inLeft][i] = V_IN;

    // 繰り返し計算
    for(i=1; i<= NX; i++) {
      for(j=1; j<= NY; j++) {
	v_tmp[i][j] = (v[i-1][j] + v[i+1][j] + v[i][j-1] + v[i][j+1])*0.25;
	//	printf("tmp=%f %f %f,", change_tmp, v[i][j], v_tmp[i][j]);
      }
    }

    // 内側境界は元に戻す
    for(i=inRight; i<=inLeft; i++) v_tmp[i][inBottom] = v_tmp[i][inTop] = V_IN;
    for(i=inBottom; i<=inTop; i++) v_tmp[inRight][i] = v_tmp[inLeft][i] = V_IN;

    // 終了条件のチェック
    for(i=1; i<= NX; i++) {
      for(j=1; j<= NY; j++) {
	change_tmp = fabs(v_tmp[i][j] - v[i][j])/v[i][j];
	if(change_tmp > change) change = change_tmp;
      }
    }

    for(i=1; i<= NX; i++) {
      for(j=1; j<= NY; j++) {
	v[i][j] = v_tmp[i][j];
      }
    }
    printf("%f interation=%i\n",change,iteration);
    iteration++;
  } while(change > MIN_CHANGE);

  outResult(v, NX, NY);
}

void outResult(double v[][NY+2], int nx, int ny) {
  int i, j;
  // 結果出力
  FILE *FP;
  FP = fopen("e_field_coax.csv","w");
  for(i=1; i<= nx; i++) {
    for(j=1; j< ny; j++) {
      fprintf(FP, "%f,", v[i][j]);
    }
    fprintf(FP,"%f\n",v[i][ny]);
  }
  fclose(FP);
}