速度が位置に依存する波動 減衰のある波動

目次 解答例(ないかも)

締切は 2002/06/07. 重点なし.

課題06_01 速度が位置に依存する波動

密度 ρ(x)が場所によって変化するような弦(片方が太くて片方が細いような弦) を伝わる波動を考えよう. 速度v と密度とは, v=( T/ρ(x))1/2 で関係している. ただし, T は定数で, 張力に対応する.

差分解法を用いて, 区間 0.0 ≤ x ≤ 1.0, 0.0 ≤ t ≤ 2.0 で波動方程式 d2f/dt2(x,t)= v(x)2 × d2f/dx2(x,t) を解こう.

ただし, (v(x))2=
0.25 * 0.2 ( x < 0.3),
0.25 * (0.6-0.4*cos((x-0.3)/0.3 * π )) (0.3 ≤ x ≤ 0.6)
0.25 * 1.0 (x > 0.6)

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

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

0( x ≤ 0.8 )
2sin3(10 π (x-0.8))(0.8 < x ≤ 1.0)
, df/dt(x,0)=
0( x ≤ 0.8 )
30 π sin2(10 π (x-0.8)) cos(10 π (x-0.8))(0.8 < x ≤ 1.0)
としよう.

これを解くのに, 課題 05_02 のクラス を継承して, 次のようなクラス Sample06_01 をつくろう.

/** 弦の線密度が変化する場合 */
class Sample06_01 extends Sample05_02 {

    /** 初期条件 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){ 
	// 埋めてね *************************************************
    }	 

    /** 場所 x によって変化する v2 . 新設*/
    double v2(double x){
	//* 埋めてね *****************************************************
    }

    /**
       関数 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] の値を設定
	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++){
// 以下のコメント部分を, 
// いまの問題にあわせて v2(x) を使って書き直してね. ***********************
	    /*
	    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(){
	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++){
// 以下のコメント部分を, 
// いまの問題にあわせて v2(x) を使って書き直してね. ***********************
	    /*
	    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 = "香山リカ";

	Sample06_01 prog = new Sample06_01(); //ここ書き直したよ.
	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 Sample06_01.java"
End:
*/

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

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

g(x)=( x < 0.3 ? 0.2 : ( x < 0.6 ? 0.6- 0.4*cos((x-0.3)/0.3 *3.14) : 1.0))
plot [0:1][-4:4] "$0" every :::i::i using 2:3  , 0.1/g(x) t "string" with line 2 , -0.1/g(x) t "" with line 2
i=i+1
pause 0.7 
if(i< $1) reread
i=0

課題06_02 減衰のある波動

差分解法を用いて, 区間 0.0 ≤ x ≤ 1.0, 0.0 ≤ t ≤ 2.0 で減衰項(抵抗)のある波動方程式 d2f/dt2(x,t)= v2× d2f/dx2(x,t) - c df/dt(x,t) を解こう. ただし, v=1.0.

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

初期条件は , 06_01と同じで, f(x,0)=

0(x ≤ 0.8 )
2sin3(10 π (x-0.8))(0.8 < x ≤ 1.0)
, df/dt(x,0)=
0( x ≤ 0.8 )
30 π sin2(10 π (x-0.8)) cos(10 π (x-0.8))(0.8 < x ≤ 1.0)
としよう.

これを解くのに, 課題 03_01 のクラスを(継承でなく)コピーして修正した, 次のようなクラス Sample06_02 をつくろう (継承でもよかったんだけど). コマンドラインオプション "-c" で, 抵抗 c の値を設定できるようにしよう.

import java.util.*;

/** 抵抗をいれた場合 */
class Sample06_02 {
    /** x の区間の下限 */
    double xmin;	

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

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

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

    /** 速度 */
    double v;

    /** 減衰項の係数 (抵抗) */
    double c;			// 追加したよ ***********************

    /** 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 = "香山リカ";

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

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

    public void initialize_parameters(String[] args){ 
	xmin=0.0; xmax=1.0; xinterval=100; // 適切な値とは限らないよ
	tmin=0.0; tmax=2.0; tinterval=100; // 適切な値とは限らないよ
	v=1.0;
	c=0.0;			// 加筆したよ

	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;
// この間に, 上の例をまねて*******************************************
// -c オプションを増設しよう. ****************************************
	    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
    }


    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] の値を設定
	// 埋めてね
	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++){
 // 変更しよう. 以下のコメント部分を, c が入るように書き換えよう ******
	    /*
	      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(){
	// 埋めてね
	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++){
 // 変更しよう. 以下の comment 部分を, c が入るように書き換えよう ******
	    /*
	      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
    }

    /** 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("# tmax =" +  tmax); // 加筆したよ
      System.out.println("# dx =" +  dx);
      System.out.println("# dt =" +  dt);
      System.out.println("# v2 =" +  v2);
      System.out.println("# c =" +  c); // 加筆したよ
      return;//nothing
    }

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

	print_previous();

	n++;

	print();

	while( tmin+(double)n*dt <= tmax ){
	    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
			   + "の実行結果です." );
	
    }

}


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

c=0.0, c=2.0, c=-1.0について, アニメーションを行おう. 差分の刻みは, x 方向が100, t 方向は, 200 くらいでどう?( コマンドラインオプションで調節してね)

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

plot [0:1][-4:4] "$0" every :::i::i using 2:3  , -2 w l 2, 2 w l 2 
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