波動方程式の差分解法

目次 解答例(ないかも)

見やすいように課題番号を書く場所を変更し, Sample のコメントを 加筆しました. 内容に変化ありません.2002/04/26 Fri 19:16

コメント のような部分は 実際には現れないここだけの説明. 1234 のような部分は, 人や場合によって変化したり, 自分で決めたりする部分. given のような部分は, 最初からあるので新たに入力する必要のない部分.

課題03-01-a ディリクレ境界条件のもとでの波動方程式の差分解法(陽解法)

差分解法を用いて, 区間 0.0 < x < 1.0, 0.0 < t < 3.0 で波動方程式 d2f/dt2(x,t)= v2 d2f/dx2(x,t) を解こう. ただし, v=2.0.

境界条件はディリクレ境界条件(固定端) f(0,t)=0, f(1.0,t)=0, とする.

初期条件は, f(x,0)=sin(πx), df/dt(x,0)=0. とする.

差分の刻みは, x 方向が40, t 方向を300 とする. つまり, 0=x0 < x1 < x2 < .. < x40=1.0. 0=t0 < t1 < t2 < .. < t300=3.0.

差分解を次のような書式で印字し, ファイルに保存しよう(例えば, datafile.1)

t   x   f(x,t)t,x の順に注意
t0 x0 0.12283      [最後の値 は f(x0,t0) の値]
t0 x1 0.12233
...(くり返し)
t0 xn 0.12277
(改行のみの行)
t1 x0 0.12299
t1 x1 0.122  
...
t1 xn 0.1221 
(改行のみの行)
t2 x0 0.12893
...
...
...
tm xn 0.2938 
		    
プログラムは, 下の例を Sample03_01.java として保存して, 埋めてね と書いてあるところを埋めれば完成するはず.
// Sample03_01
import java.util.*;

/** 波動方程式を解くクラス */
public class Sample03_01 {

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

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

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

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

    /** 速度 */
     double v;

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

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

    /** (v dt/dx)*(v dt/dx). 速度の2乗を無次元化したもの */
    double v2;			

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

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


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

    /** 時刻 n-1 での解 */
     double[] previous_f;

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

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

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

    /** 初期条件 df/dt(x,tmin) */
     double initial_dfdt(double x){ 
	 // 埋めてね
     }

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

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



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

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

	prog.initialize_parameters();
	prog.initialize_fields();
	prog.exec();
	return;//nothing
    }

    /** 
	パラメターの初期化.
	dx, dt を計算. 
     */
    public void initialize_parameters(){
	xmin=0.0; xmax=1.0; xinterval=40;
	tmin=0.0; tmax=3.0; tinterval=300;
	v=2.0;

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

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

	n=0;
	// previous_f[0] .. previous_f[xinterval+1] の値を設定
	// current_f[0] .. current_f[xinterval+1] の値を設定
	// 埋めてね
	return;//nothing
    }


    /** 
	時刻 n-1, n の f  ( previous_f, current_f ) から,
	時刻 n+1 の f ( next_f ) を計算.
	
	その後, 1つずつずらす
	current_f <-- next_f
	previous_f <-- current_f
	捨てる  <--- previous_f
    */
    void evolve(){
	// 埋めてね
        n++;	      
	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
    }

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

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

    /** ループで実行 */
    public void exec(){
	print_info();

	print_previous();

	n++;

	print();

	while( 
	      //埋めてね
	      ){
	    evolve();
	    print();
	}
	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
			   + "の実行結果です." );
	
    }

}



  • 上で, datafile.1に保存したデータを, gnuplot を使ってグラフやアニメーションにしてみよう.

    まず初期条件が正しいか? x 対 f のグラフを描いてみよう.

    s1609h111% gnuplot
    gnuplot> plot "datafile.1" every :::0::0 using 2:3 
    # datafile.1 の 0番目のブロックを, 2列目を x, 3列目を y と見なしてグラフ作成
    
    

    コマ送り. 上向きカーソルキーを使うと簡単.

    s1609h111% gnuplot
    gnuplot> i=0
    gnuplot> plot "datafile.1" every :::i::i using 2:3  ; i=i+1
    gnuplot> plot "datafile.1" every :::i::i using 2:3  ; i=i+1
    gnuplot> plot "datafile.1" every :::i::i using 2:3  ; i=i+1
    ....
    

    全部の時刻をまとめて書いてみよう.

    s1609h111% gnuplot
    gnuplot> plot "datafile.1" using 2:3 
    

    次にアニメーションする. 次の内容を, ファイル looper.gp として保存.

    plot [0:1][-2:2] "$0" every :::i::i using 2:3  ,-1,1
    i=i+1
    pause 0.4 
    if(i< $1) reread
    i=0
    
    次に実行.
    s1609h111% gnuplot
    gnuplot> i=0
    gnuplot> call "looper.gp" "datafile.1" 300
    gnuplot> 
    # looper.gp で, datafile.1 を,  300コマのアニメーション
    

    課題03-01-b差分解法の安定性

    t 方向の差分の刻みを, 300 から 240, 239, 200 に変更して, 同じことをやってみよう. プログラムのファイル名は Sample03_01.java のままで中身を変更してよい.

    課題03-02-a クラスの継承

    クラス Sample03_01 の一部分だけを変更, 追加して 別のプログラムにするには, 継承(inheritance)という仕組みを使うとよい.

    下のプログラム例 Sample03_02.java の指定の場所に加筆して, 差分解法を用いて, 区間 0.0 < x < 1.0, 0.0 < t < 2.0 で波動方程式 d2f/dt2(x,t)= v2 d2f/dx2(x,t) を, 初期条件だけ変えて解こう. ただし, v=2.0.

    境界条件はディリクレ境界条件 f(0,t)=0, f(1.0,t)=0 のままとする.

    初期条件は, 0.5< x < 0.75 f(x,0)=(1-cos(8π(x-0.5)))/2, df/dt(x,0)=8πsin(8π(x-0.5)), それ以外で, f(x,0)=0, df/dt(x,0)=0 とする.

    差分の刻みは, x 方向に 40, t 方向に 160 および 200 を試そう.

    差分解を上と同じ書式で印字し, ファイルに保存しよう(例えば, datafile.2)

    // Sample03_02.java
    
    
    /** Sample03_01 を継承して, 別の条件で波動方程式を解く */
    public class Sample03_02 extends Sample03_01 {
    // extends が継承のキーワード
    
        /** 初期条件 f(x,tmin) */
        double initial_f(double x){ 
    	// 埋めてね
        }
    
        /** 初期条件 df/dt(x,tmin) */
         double initial_dfdt(double x){ 
    	// 埋めてね
        }
    
        /** ディリクレ境界条件 f(xmin,t) */
        double fmin(double t){ 
    	// 埋めてね
        }
    
        /** ディリクレ境界条件 f(xmax,t) */
        double fmax(double t){ 
    	// 埋めてね
        }
    
        public void initialize_parameters(){
    	// 変更して実験してね
    	xmin=0.0; xmax=1.0; xinterval=40;
    	tmin=0.0; tmax=1.0; tinterval=300;
    	v=2.0;
    
    	dx=(xmax-xmin)/(double)xinterval;
    	dt=(tmax-tmin)/(double)tinterval;
    	v2=(v*dt/dx)*(v*dt/dx);
    
        }
    
        public static void main(String[] args){
    	String myName = "香山リカ";
    	Sample03_02 prog= new Sample03_02();
    	RunInfo info= new RunInfo(myName, prog.getClass().getName());
    	info.print();
    
    	prog.initialize_parameters();
    	prog.initialize_fields();
    	prog.exec();
        }
    
    
        // 他のメソッド, フィールドは, Sample03_01 のものが使われる.
    }
    
    
    
    
    		    

    課題03-02-b

    同様にアニメーションしてみよう. どんな様子が見えるでしょう?
    Copyright © 2002 Saburo Higuchi. All rights reserved.
    Saburo HIGUCHI, hig mail address
    最後の更新は次の日時よりあとのはず…: 2002/04/12 Fri 22:25