拡散方程式の陰解法 Javaの標準入力(gnuplot 用のデータ形式)

目次 解答例(ないかも)

締切は 2002/07/12. お薦め課題: なし.

課題じゃないけど, きょうの講義の quiz を Mathematica で解いたです. ディスクにセーブして Mathematica で開いてね.

課題12_01 拡散方程式の陰解法

区間 0.0 ≤ x ≤ 1.0, 0.0 ≤ t ≤ 1.0 で拡散方程式 df/dt(x,t)= D · d2f/dx2(x,t) 考えよう. ただし, D=0.2 およびD=0.8.

x=0での境界条件は, ディリクレ境界条件 f(0,t)=tとしよう(注意: はじめて時間に依存します)

x=1での境界条件は, ノイマン境界条件(自由端) df/dx(1,t)=0としよう.

初期条件は f(x,0)=0 としよう.

これを講義で説明する陰解法を用いて解こう. 課題11_01に載っているクラス TridiagMatrix を用いて, 次のようなプログラムにしよう. 注: 課題12_01 は, 課題11_01 のTridiagMatrix.java を完成させなくても( 課題11_01 のプログラム例をそのまま保存しただけで)作成可能です.

/** 陰解法で Dirichlet 境界条件のついた拡散方程式を解くクラス */
public class Sample12_01 extends Sample07_01 {

    // 追加したよ ********************************************************
    /** 陰解法に出てくる1次方程式の係数行列 */
    TridiagMatrix c;	   // こういうふうに, 継承ではメンバーも追加できる.

    /** 初期条件 f(x,tmin) */
     double initial_f(double x){
	 // 埋めてね *****************************************************
     }


    /** ディリクレ境界条件 f(xmin,t) */
    double fmin(double t){ 
	return t;
    }

    /** ディリクレ境界条件 f(xmax,t) */
    // x=xmin 側はノイマンなので使わない.
//       double fmax(double t){ 
//  	 // 埋めてね. *************************************************
//       }


    /** プログラムを実行するとこのメソッドが実行される. */ 
    public static void main(String[] args){
	String myName = "香山リカ";

	Sample12_01 prog= new Sample12_01(); // 変更したよ **************
	RunInfo info= new RunInfo(myName, prog.getClass().getName());
	info.print();

	prog.exec(args);
	return;//nothing
    }

    /** プログラムの本体 */
    public void exec(String[] args){

	initialize_parameters(args);      // コマンドライン引数を読んで
				          // パラメター(dt,dx等)を設定
	print_info();			  // パラメターを印字

	initialize_matrix();	// c を設定.  追加したよ****************

	initialize_fields();	// initial_f から, 最初の current_f を計算
	print(); // current_f を印字

	while( n < tinterval ){
	    evolve();		// 差分公式で next_f を計算し, 
				// current_f <-- next_f, すてる <-- current_f
	    print();		// current_f を印字
	}
	return;//nothing
    }


    void initialize_matrix(){
	final int dim=xinterval+1-1; // 変化したよ

	double[] diag = new double[dim];
	double[] upper = new double[dim];
	double[] lower = new double[dim];

	/* 
	   今の問題にあわせて, d1 などで,
	   連立方程式の係数行列 c を初期化しよう.
	   片側がディリクレ, 片側がノイマン境界条件なことに注意.
	*/
	// 埋めてね ****************************************************

	c=new TridiagMatrix(dim,diag,upper,lower);

    }



    void evolve(){

	/*
	 差分公式(陰解法)で, next_f[0],...next_f[xinterval] を決める.  
	 (i=0, xinterval は上で別に扱った)
	 前回は陽解法, 今回は陰解法なので変更が必要. 

	*/

	// 埋めてね **************************************************

	// 時刻 n を進め, 今計算した next を current とする. 変更なし.
	n++;
	for(int i=0;i <= xinterval;i++){
	    current_f[i]=next_f[i];
	}

	return;//nothing
    }
}

/* emacs setting
Local Variables:
mode: java
compile-command: "javac -deprecation Sample12_01.java"
End:
 */

差分の刻みを, x 方向を20, t 方向を 50 としよう. 拡散係数 D が, 0.2 の場合と 0.8 の場合を比べよう(-d オプションで変更できるのだった)

課題12_02 Javaの標準入力(gnuplot 用のデータ形式)

これまで, Java で作成したプログラムに計算結果をテキストファイ ルとして出力させてきた. 改行のみの行を空行という事にすると, 空行ではさまれた部分(ブロックとよぼう)がひとつの時刻でのデータである. ブロック内の各行は, 特定の位置x=i*dxに対応している. 行内での0番目のフィールドは時刻t 1番目のフィールドは位置x 2番目のフィールドはf(x,t) だった.

これまでは, テキストファイルを, gnuplot や Mathematica に読み込んでグラフ化してきた. テキストファイルとなった計算結果は, その他のプログラムに読み込んで加工することもできる(例. Excel, Word, C で書いたプログラム ... ). ここでは, その例として, 計算結果のテキストファイルを, Java で作成した別のプログラムに読み込んで解析することを行おう.

次のプログラムを(変更せず)コンパイルしよう.

import java.io.*;

/** 
    これまでの数値計算の結果のテキストファイルを
    標準入力から受け取って解析するプログラム
*/
public class Sample12_02 {


    /** 読み込まれる最大のブロック(空行で区切られる部分)数 */
    final int maxblock=200;

    /** 読み込まれる最大の行数(ブロック内での) */
    final int maxrow=21;

    /** 読み込まれる最大の数値(フィールド)の個数(行内) */
    final int maxcol=3;

    /*
      maxblock=3
      maxrow=2
      maxcol=5   で読み込めるぎりぎりのデータは次のようなもの

      1 2 3 2 1
      3 3 8 9 9

      9 9 3 8 8
      8 8 7 8 8

      8 8 8 7 8   ← この7 は, 3番目のブロックの 0 行目の 4番目のフィールド
      8 0 9 8 8
     */


    /** メインメソッド. 実行するとこれが呼ばれる */
    public static void main(String[] args) {
	String myName = "香山リカ";

	Sample12_02 prog= new Sample12_02(); // 変更したよ **************
	RunInfo info= new RunInfo(myName, prog.getClass().getName());
	info.print();

	prog.exec(args);
	return;//nothing
    }

    /** プログラムの実体 */
    void exec(String[] args){

	double m[][][]= new double [maxblock][maxrow][maxcol];

	/* initialize array m */
	for(int b=0; b < maxblock ; b++){
	    for(int r=0; r < maxrow; r++){
		for(int c=0; c < maxcol; c++){
		    m[b][r][c]=0.0;
		}
	    }
	}

	readData(m);

	// ***解析開始ポイント*** ****************************************

	analyzeData(m);
	
    }

    /** 標準入力かデータを受け取り,  配列 m に格納する */
    void readData(double [][][] m){
	InputStreamReader f=new InputStreamReader(System.in);
	StreamTokenizer s= new StreamTokenizer(f);

	s.parseNumbers();	// 数値も理解してね
	s.eolIsSignificant(true); // 改行もトークンとして扱う
	s.ordinaryChar((int)'#'); // # は空白文字として扱わない
				// 以下で, comment の開始という意味を与える

	int b=0;		// いま何番目のブロックをみているか
	int r=-1;		// いま何番目の行をみているか
	int c=0;		// いま何番目のフィールドをみているか

	int type;		// 今みている部分の種別. 数値, 単語など

	boolean inComment=false; // 状態変数. # と 改行の間をみている
	int consecutiveNewLines=1; // 現在, 改行が何個連続しているか

	while( true ){
	    try {
		if( (type=s.nextToken())==StreamTokenizer.TT_EOF ){
		    break;
		}
	    } catch ( IOException e ){
		e.printStackTrace();
		break;
	    }

	    switch ( type ){
	    case StreamTokenizer.TT_EOF:
		break; // this cannot happen
	    case StreamTokenizer.TT_EOL:
		if ( !inComment ){
		    consecutiveNewLines++;
		}
		inComment=false;
		break;
	    case  StreamTokenizer.TT_NUMBER:
		if( inComment ){
		    break;
		}

		if( consecutiveNewLines >=2 ){
		    b++;
		    c=0;
		    r=0;
		} else if ( consecutiveNewLines==1 ){
		    r++;
		    c=0;
		}
    
		if( b < maxblock && r>=0 && r < maxrow &&  c < maxcol ){
		    m[b][r][c]=s.nval ;
		}

		consecutiveNewLines=0;
		c++;
		
		break;
	    case '#':
		inComment=true;
		break;
	    default:
		break;
	    }
		
	}
    }

    /** 引数 m の配列に納められたデータを解析して標準出力に送る */
    void analyzeData(double [][][] m){
    /* 
       現状では順に print するだけ.
       この中全部書き換えてね. ********************************************
    */
	
//  	for(int b=0; b < maxblock; b++){
//  	    for(int r=0; r < maxrow; r++){
//  		for(int c=0; c < maxcol; c++){
//  		    System.out.println(b + " " + r + " " + c + " " + m[b][r][c] );
//  		}
//  	    }
//  	}


	return ;
    }
}

/* emacs setting
Local Variables:
mode: java
compile-command: "javac -deprecation Sample12_02.java"
End:
*/

そして次のようにして実行しよう.
s1609h017% java Sample07_01 > datafile.1
s1609h017% cat datafile.1  | java Sample12_02 > datafile.2
または
s1609h017% java Sample07_01 | java Sample12_02  > datafile.2

readData()メソッドが完了した時点( ***解析開始ポイント*** という行) で, 3次元配列 m に, すべてのデータが格納されている. その後のanalyzeData() メソッドでは, 単に,

ブロック番号 行番号 フィールド番号 数値
を出力している.

このプログラムの, analyzeData()メソッドを書き換えて, 次のような出力を行うようにしよう.

(時刻 t)  (時刻 t での最小値 minxf(x,t)) (時刻 t での最大値 maxxf(x,t))  (時刻 t でのf の平均値) 
0.5 0.8 1.5 1.2                                       
1.0 0.9 1.4 1.2
...
f の平均値は, ∫ f(x,t) dx/(xmax-xmin)で, {0.5*f(0,t) + 1.0 * f(dx,t) + 1.0*f(2dx,t) + ... + 1.0 *f((xinterval-1)dx,t) + 0.5 * f(xinterval*dx,t)}/xinterval で求めよう.

完成したら, さらに Sample07_01 (ディリクレ境界条件)と Sample09_01 (ノイマン境界条件) のそれぞれ出力したテキストファイルについて, 解析を行おう. さらに,

s1609h017% gnuplot
gnuplot> plot "datafile.2" using 1:2, "" using 1:3 , "" using 1:4
によって可視化しよう. 意味を考えよう.

さらに興味と暇のある人は, データを配列 m[][][] に保存するよなメモリ浪 費をせずに同様の出力を得る方法を工夫しよう.


Copyright © 2002 Saburo Higuchi. All rights reserved.
Saburo HIGUCHI, hig mail address
最後の更新は次の時刻以降 2002/06/28 Fri 10:04