ノイマン境界条件のもとでの波動方程式 拡散方程式の差分解法

目次 解答例(ないかも)

締切は 2002/06/28. お薦め課題 09_01

課題09_01 ノイマン境界条件のもとでの拡散方程式の差分解法

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

x=0, 1での境界条件は, ノイマン境界条件(自由端) df/dx(0,t)=0, df/dx(1,t)=0としよう. ノイマン境界条件の扱い方は授業で説明します.

初期条件は f(x,0)=1-cos(π x) + (1-x)2x2 としよう.

これを解くのに, クラス
Sample07_01 を継承して, 次のようなクラスを完成させよう.
// Sample09_01
import java.util.*;

/** Neumann 境界条件のついた拡散方程式を解くクラス */
public class Sample09_01 extends Sample07_01 {

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


    // 今回はノイマン境界条件を考えるので,
    // fmin, fmax は使わない(のでどう定義してあっても無関係)
    // 特にオーバーライドもしない.
//      /** ディリクレ境界条件 f(xmin,t) */
//       double fmin(double t){ 
//  	 // 埋めてね *************************************************
//       }

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

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

	Sample09_01 prog= new Sample09_01();// ここ書き換えたよ.
	RunInfo info= new RunInfo(myName, prog.getClass().getName());
	info.print();

	prog.exec(args);
	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;

	// fmin, fmax 使わないから*************************************
	// 下のコメント部分書き直さないといけないよね.

        // 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 の両端の値を決める.
	// fmin, fmax は使わずに, 
	// 講義で説明した式を使うようにすればいいはず.
	// 07_01 が出来てる人は, ここだけ変更すればいいはず.
	// 埋めてね   *************************************************

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

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

	return;//nothing
    }
    

}

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


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

課題 08_01 から plotsample.nb をもう一度取ってきて, グラフ化してみよう. 結果はこんな感じかも. 部屋と温度のたとえで言うと何が起きてるの?

課題09_02 ノイマン境界条件のもとでの波動方程式の差分解法

ノイマン境界条件のもとで波動方程式も解こう.

差分解法を用いて, 区間 0.0 ≤ x ≤ 2.0, 0.0 ≤ t ≤ 2.5 で波動方程式 d2f/dt2(x,t)= v2 d2f/dx2(x,t) を解こう. ただし, v=1.0. 注意: x の範囲がはじめて 2.0 までになりました. どこを直せばいいでしょう?

x=0, 1での境界条件は, ノイマン境界条件(自由端) df/dx(0,t)=0, df/dx(1,t)=0としよう. ノイマン境界条件の扱い方は授業で説明します.

初期条件は f(x,0)=

0.2*(sin(2 π x))3(1.0 < x < 2.0),
0 (それ以外)
df/dt(x,0)=
1.2 π (sin(2 π x))2cos(2 π x)(1.0 < x < 2.0),
0 (それ以外)
としよう.

これを解くために, クラス Sample03_01 を継承して, 下のようなクラス Sample09_02 を作ろう. なお, 境界の点i=0,xintervalでは, 拡散方程式の場合と同様, 差分の式を, 何とか_f[-1]何とか_f[+1] と見なして適用, 何とか_f[xinterval+1]何とか_f[xinterval-1] と見なして適用すればよい.

/** ディリクレ境界条件のもとでの波動方程式. 自由端での反射 */
class Sample09_02 extends Sample03_01 {

    /** 初期条件 f(x,tmin) */
    double initial_f(double x){
	// たぶんこれで正しいはず.
	double f;
	if( x > 1.0 && x < 2.00 ){
	    f=Math.sin(2.0 * Math.PI * (x-1.0));
	    f=f*f*f*0.2;
	} else {
	    f=0.0;
	}
	return f;
    }

    /** 初期条件 df/dt(x,tmin) */
    double initial_dfdt(double x){ 
	// たぶんこれで正しいはず.
	double f;
	if( x > 1.0 && x < 2.00 ){
	    f= 3.0 * 2.0 * Math.PI* v 
		* Math.sin(2.0 * Math.PI * ( x-1.0 ))
		* Math.sin(2.0 * Math.PI * ( x-1.0 ))
		* Math.cos(2.0 * Math.PI * ( x-1.0 ))
		*0.2;
	 } else {
	     f= 0.0;
	 }
	 return f;
    }
    
    // 今回はノイマン境界条件を考えるので,
    // fmin, fmax は使わない(のでどう定義してあっても無関係)
    // 特にオーバーライドもしない.
    /** ディリクレ境界条件 f(xmin,t) */
    /*
    double fmin(double t){ 
	return 0;
    }	
    double fmax(double t){ 
	return 0;
    }	 
    */

    /** パラメターの初期化. コマンドラインオプションも読む */
    public void initialize_parameters(String[] args){ 
	// このコメント部分書き直してね ************************************
	//xmin=0.0; xmax=1.0; xinterval=100;

	tmin=0.0; tmax=5.3; tinterval=100;
	v=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 'v':
		i++;
		v=Double.parseDouble(args[i]);
		break;
	    default:
		System.err.println("間違ったオプションです.");
		break;
	    }
	    i++;
	}    
	//  加筆おしまい
	 
	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] の値を設定
	// fmin, fmax 使わないんだから
	// 下のコメント部分(03_01のまま)は書き直しがいるはず.**************
//  	previous_f[0]=fmin(tmin);
//  	previous_f[xinterval]=fmax(tmin);
//  	for(int i=1;i < xinterval;i++){
//  	    previous_f[i]=initial_f( xmin + (double)i *dx );
//  	}


//  	current_f[0]=fmin(tmin + dt);
//  	current_f[xinterval]=fmax(tmin + dt);
//  	for(int i=1;i < xinterval;i++){
//  	    current_f[i]=previous_f[i]
//  		+ initial_dfdt( xmin + (double)i * dx ) * dt
//  		+ 1.0/2.0 * v2*
//  		(previous_f[i+1] - 2.0 * previous_f[i] + previous_f[i-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(){
	// fmin, fmax 使わないんだから
	// 下のコメント部分(03_01のまま)は書き直しがいるはず.
//  	next_f[0]        =fmin(tmin + (double)(n+1)*dt);
//  	next_f[xinterval]=fmax(tmin + (double)(n+1)*dt);

//  	for(int i=1;i < xinterval;i++){
//  	    next_f[i]=
//  		v2 * (current_f[i+1]+current_f[i-1])
//  		+2.0*(1.0-v2)*current_f[i]
//  		-previous_f[i];
//  	}

  	n++;
  	for(int i=0;i <= xinterval;i++){
  	    previous_f[i]=current_f[i];
  	    current_f[i]=next_f[i];
  	}
	return;//nothing
    }


    public static void main(String[] args){
	String myName = "香山リカ";

	Sample09_02 prog = new Sample09_02(); //ここ書き直したよ.
	RunInfo info= new RunInfo(myName, prog.getClass().getName());
	info.print();

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

}

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

差分の刻みは, x方向は100, t方向は150 および 400 で実行しよう.

グラフ化しよう. 同じ初期条件で, 境界条件が異なる(ディリクレ境界条件,固定端)の場合の下の結果と比較しよう. 注: 下のアニメーションは課題 09_02 の答えではありません. 固定端の場合の結果です.


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