拡散方程式の差分解法

目次 解答例(ないかも)

締切は 2002/06/21. お薦め課題 07_01

課題07_01 拡散方程式の差分解法

差分解法を用いて, 区間 0.0 ≤ x ≤ 1.0, 0.0 ≤ t ≤ 1.0 で拡散方程式 df/dt(x,t)= D • d2f/dx2(x,t) を解こう. ただし, D=0.2

x=0, 1での境界条件は, ディリクレ境界条件(固定端) f(0,t)=0, f(1,t)=0としよう.

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

これを解くのに, 次のようなクラス Sample07_01 を完成させよう.

// Sample07_01
import java.util.*;

/** 拡散方程式を解くクラス */
public class Sample07_01 {

    /** x の区間の下限 */
     double xmin;	

    /** x の区間の上限 */
     double xmax;

    /** t の区間の下限 */
     double tmin;

    /** t の区間の上限 */
     double tmax;

    /** 拡散定数 D */
     double d;

    /** x 方向の分割. 点の個数はこれ + 1 */
     int xinterval;

    /** t 方向の分割. 点の個数はこれ + 1 */
     int tinterval;

    /** 無次元化された拡散定数 D1=D dt/dx/dx */
    double d1;

    /** 刻み幅. xmin, xmax, xinterval から決まる */
     double dx;

    /** 刻み幅. tmin, tmax, tinterval から決まる */
     double dt;

    /** 時刻のカウンタ */
     int n;

    /** 時刻 n での解 */
     double[] current_f;

    /** 時刻 n+1 での解 */
     double[] next_f;

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

    /** ディリクレ境界条件 f(xmin,t) */
     double fmin(double t){ 
	 // 埋めてね *************************************************
     }

    /** ディリクレ境界条件 f(xmax,t) */
     double fmax(double t){ 
	 // 埋めてね *************************************************
     }

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

	Sample07_01 prog= new Sample07_01();
	RunInfo info= new RunInfo(myName, prog.getClass().getName());
	info.print();

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

    /** 
	パラメターの初期化.
	-X xinterval -T tinterval -t tmax -d d
	dx, dt, d1 を計算. 
     */
    public void initialize_parameters(String[] args){
	xmin=0.0; xmax=1.0; xinterval=20;
	tmin=0.0; tmax=1.0; tinterval=200;
	d=0.2;

	int i=0;
	while ( (i < args.length ) && ( (args[i]).startsWith("-"))){
	    switch( args[i].charAt(1) ){
	    case 'X':
		i++;
		xinterval=Integer.parseInt(args[i]);
		break;
	    case 'T':
		i++;
		tinterval=Integer.parseInt(args[i]);
		break;
	    case 't':
		i++;
		tmax=Double.parseDouble(args[i]);
		break;
	    case 'd':
		i++;
		d=Double.parseDouble(args[i]);
		break;
	    default:
		System.err.println("間違ったオプションです.");
		break;
	    }
	    i++;
	}    

	dx=(xmax-xmin)/(double)xinterval;
	dt=(tmax-tmin)/(double)tinterval;
	d1=d*dt/dx/dx;
	return;//nothing
    }

    /**
       関数 f を表す配列確保.
       初期値から, n=0,  での f の値を current_f に代入.
     */
    public void initialize_fields(){
	// 配列を確保.
	next_f = new double [xinterval+1];
	current_f = new double [xinterval+1];
	n=0;
	// current_f[0] .. current_f[xinterval+1] の値を設定

	current_f[0]=fmin(tmin);
	current_f[xinterval]=fmax(tmin);
	for(int i=1;i < xinterval;i++){
	    current_f[i]=initial_f(xmin + (double)i * dx );
	}
	return;//nothing
    }


    /** 
	時刻 n の f  (  current_f ) から,
	時刻 n+1 の f ( next_f ) を計算.
	
	その後, 1つずつずらす
	current_f <-- next_f
	捨てる  <-- current_f
    */
    void evolve(){


	// 境界条件を使って next_f の両端の値を決める.
	// 埋めてね   *************************************************

	// 差分公式を使う. i=0, xinterval は上で別に扱った.
	// 埋めてね   *************************************************

	// 時刻 n を進め, 今計算した next を current とする.
	// 埋めてね   *************************************************

	return;//nothing
    }
    
    /** current_f を gnuplot に読めるように印字 */
    public void print(){
	System.out.println("#" + n );
	for(int i=0; i <= xinterval; i++){
	    System.out.println((tmin+(double)n*dt)  + " " 
			       + (xmin+(double)i*dx) + " "
			       + current_f[i]);
	}
	System.out.println();
	return;//nothing
    }


    /** 計算の条件を出力 */
    public void print_info(){
      System.out.println("# diffusion =" +  d);
      System.out.println("# tinterval =" +  tinterval);
      System.out.println("# xinterval =" +  xinterval);
      System.out.println("# tmax =" +  tmax);
      System.out.println("# dx =" +  dx);
      System.out.println("# dt =" +  dt);
      System.out.println("# d1 =" +  d1);
      return;//nothing
    }

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

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


	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
    }

}










// 以下はいじらなくてよい -----------------------------------------------
// この中で error が出るとしたら, ここより上で } が閉じてない.
class RunInfo {
    String username;
    String progname;

    RunInfo(String u,String p){
	this.username=u;
	this.progname=p;
    }

    public void print(){
	// 学籍番号, 氏名を表示する. 
	System.err.println("# " 
			   + "このプログラムの作成者は "
			   + System.getProperty("user.name") 
			   + " " 
			   + this.username  + "です.");

	// クラスと日時を表示する 
	System.err.println("# "
			   +"課題" + this.progname  + "の"
			   +  (new Date()).toString() // from java.util.Date
			   + "の実行結果です." );
	
    }

}


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

差分の刻みは, x 方向が20, t 方向を 200 でどう? ( コマンドラインオプションで調節してね)

次の looper.gp を利用して, gnuplot でアニメーションしよう.

plot [0:1][-0.1:0.25] "$0" every :::i::i using 2:3, 0 title  "" 
i=i+1
pause 0.7 
if(i< $1) reread
i=0


Copyright © 2002 Saburo Higuchi. All rights reserved.
Saburo HIGUCHI, hig mail address
最後の更新は次の時刻以降 2002/05/19 Sun 16:43