3重対角行列クラス 拡散方程式の陰解法

目次 解答例(ないかも)

11_01(お薦め)の締切は 2002/07/12. 11_02(お薦め)の締切は 2002/07/12(最終回).

課題11_01 3重対角行列クラス

次の, 3重対角行列の計算をするためのクラスを考える. TridiagMatrix.java として保存しよう.

/** double の3重対角行列を表すクラス. 疎行列として扱われる */
class TridiagMatrix {
    /** 行列の次元*/
    private int dim;

    /** 対角に近い成分を表す配列. 長さ dim */
    private double d[],u[],l[];
    /*-
      d[0] u[0]    0    0    0 ...                0
      l[1] d[1] u[1]    0    0 ...                0
      0    l[2] d[2] u[2]    0 ..                 0
      0       0 l[3] d[3] u[3]                    0
                                                  0
                                  d[dim-2] u[dim-2]
      0   ....                    l[dim-1] d[dim-1]

	   l[0], u[dim-1] は 0.0 に保たれる.
     */


    /** n x n で, 対角が d, 上が u, 下が l の3重対角行列のコンストラクタ */
    TridiagMatrix(int n, double d[], double u[], double l[]){

	this.dim=n;

	if( dim < 2 || dim != d.length ||  dim != u.length || dim != l.length ){
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}

	this.d = new double[dim];
	this.u = new double[dim];
	this.l = new double[dim];



	for(int i=0; i < n ; i++){
	    this.d[i]=d[i];
	    this.u[i]=u[i];
	    this.l[i]=l[i];
	}
	this.l[0]=0.0; this.u[this.dim-1]=0;

    }

    /** 
    対角が d, 上が u, 下が l の3重対角行列のコンストラクタ. 
    次元は配列のサイズから取る
    */
    TridiagMatrix(double d[], double u[], double l[]){
	this( d.length, d,u,l );
    }

    /** 
    次元がn の3重対角行列のコンストラクタ. 
    成分は初期化されない
    */
    TridiagMatrix(int n){
	if( n>1 ){
	    this.dim=n;
	    this.d = new double[n];
	    this.u = new double[n];
	    this.l = new double[n];
	    this.l[0]=0.0; this.u[this.dim-1]=0;
	} else {
	    /* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}
    }

   
    /** 次元を返す */
    public int getDim(){
	return dim;
    }

    /** (i,j) 成分を返す. i,j=0,1,2,...,dim-1 */
    public double getElem(int i, int j){
	double ret=0.0;
	if( 0 <= i && 0 <= j && i< dim && j < dim ){
	    if( i==j ){
		ret= d[i];
	    } else if ( i-j==1){
		ret=l[i];
	    } else if ( i-j==-1){
		ret=u[i];
	    } else {
		ret=0.0;
	    }

	} else {
	    /* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	    ret=0.0;
	}

	return ret;
    }

    /** x と y の和を生成して返す */
    public static TridiagMatrix add(TridiagMatrix x, TridiagMatrix y){
	if ( x.dim != y.dim ){
	    /* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}

	TridiagMatrix z= new TridiagMatrix( x.dim );
	
	for( int  i= 0;  i < z.dim ; i++){
	    z.d[i]=x.d[i]+y.d[i];
	    z.u[i]=x.u[i]+y.u[i];
	    z.l[i]=x.l[i]+y.l[i];
	}
	z.l[0]=0.0; z.u[z.dim-1]=0.0;
	
	return z;
    }

    /** 自分自身に引数 y を加えた3重対角行列を生成して返す */
    public TridiagMatrix add(TridiagMatrix y){
	if ( this.dim != y.dim ){
	    /* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}

	TridiagMatrix z= new TridiagMatrix( y.dim);
	// 埋めてね  ************************************************* 
	return z;
    }

    /** 自分自身と, 引数の配列 v の表すベクトルをかけて得たベクトルを返す*/
    public double[] multiply(double[] v){
	if( this.dim !=v.length ){
	    /* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}

	double [] vret = new double [dim];
	// 埋めてね  ************************************************* 
	return vret;
    }

    /** 
	自分自身と, 引数の行列 y をかけて得た行列のうち,
	対角と隣の斜め線を取り出して, 3重対角行列と見なして返す.
    */
    public TridiagMatrix multiplyProject(TridiagMatrix y){

	if ( this.dim != y.dim ){
	    /* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}

	TridiagMatrix z= new TridiagMatrix( this.getDim());
	// 埋めてね *************************************************
	return z;
    }

    /**
       自分自身 M を
       3重対角行列 L, U に
       LU 分解.
       すなわち, M=LU.
       {L,U} の形の配列を返す. 自分自身を変えない 
       L の対角成分は 1 に固定
    */
    public TridiagMatrix[]  decomposeLU(){

	TridiagMatrix[]  ret = new TridiagMatrix[2];

	for(int i=0; i < 2 ; i++){
	    ret[i]= new TridiagMatrix(this.dim);
	}

	TridiagMatrix lmat=ret[0];
	TridiagMatrix umat=ret[1];


	/* あらかじめ決まってるところは代入してしまう*/
	for(int i=0; i < this.dim ; i++){
	    lmat.d[i]=1.0;
	    lmat.u[i]=0.0;
	    umat.u[i]=this.u[i];
	    umat.l[i]=0.0;
	}
	
	/* 最初だけちょっと特別だが */
	umat.d[0]=this.d[0];
	/* あとは順に解いていく */
	for(int i=1; i < this.dim ; i++){
	    lmat.l[i]=this.l[i]/umat.d[i-1];
	    umat.d[i]=this.d[i]-lmat.l[i]*umat.u[i-1];
	}

	/* 最初から決まっている部分 */
	umat.l[0]=0.0;
	umat.u[dim-1]=0.0;
	lmat.l[0]=0.0;
	lmat.u[dim-1]=0.0;

	return ret;
    }


    /** 1次方程式  this * x  == rhs  を解いて, 配列 x を返す. */
    public double[] applyInverse(double[] rhs){
	if( this.dim!=rhs.length ){
	    /* エラー*/
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}

	/* まず LU 分解する */

	TridiagMatrix[] lumats=decomposeLU();

	TridiagMatrix lmat= lumats[0];
	TridiagMatrix umat= lumats[1];

	double[] v1 = lmat.applyInverseLower(rhs);
	
	double[] v2 = umat.applyInverseUpper(v1);

	return v2;
    }

    /** 
    this が下三角かつ3重対角のときに,
   1次方程式  this * x  == rhs  を解いて, 配列 x を返す. 
    */
    public double[] applyInverseLower(double[] rhs){
	if( this.dim != rhs.length ){
	    /* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}

	for(int i=1; i < dim ; i++){
	    if( u[i]!=0.0 ){
		/* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	    }
	}

	double [] v = new double[this.dim];
	
	v[0]=rhs[0];
	for(int i=1; i < this.dim; i++){
	    v[i]=rhs[i]-this.l[i]*v[i-1];
	}
	
	return v;
    }

    /** 
    this が上三角かつ3重対角のときに,
    1次方程式  this * x  == rhs  を解いて, 配列 x を返す. 
    */
    public double[] applyInverseUpper(double[] rhs){
	if( this.dim != rhs.length ){
	    /* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	}

	for(int i=0; i < dim-1 ; i++){
	    if( l[i]!=0.0 ){
		/* エラー */
	    System.err.println("Tridiag Matrix のエラー");
	    /* 本当は, "例外を投げる" べき */
	    }
	}

	double [] v = new double[dim];

	v[dim-1] = rhs[dim-1]/this.d[dim-1];
	for(int i=dim-2; i>=0 ; i--){
	    v[i]=(rhs[i]-this.u[i]*v[i+1])/this.d[i];
	}

	return v;
    }
    
    /** 文字列に変換. 表型の表示 */   
    public String toString(){
	String buf="";

	for(int i=0; i < dim ; i++){
	    for(int j=0; j< dim ; j++){
		buf += getElem(i,j) ;
		buf += " ";
	    }

	    buf += "\n";
	}

	return buf;
    }


}

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

メソッド add, multiply, multiplyProjectは不完全である. 次のテストプログラムが正しく動作するように, 完成させよう. 注: 課題11_02 は, TridiagMatrix.java を完成させなくても(下に載ってるのをそのまま保存しただけで)実行可能です.
/**
   TridiagMatrix クラスのテストのためのプログラム
*/
class Sample11_01 {

    /** メインメソッド */
    public static void main(String[] args){


	double u1[] = {1.0,1.0,1.0,1.0,0.0};
	double d1[] = {4.0,4.0,4.0,4.0,4.0};
	double l1[] = {0.0,2.0,2.0,2.0,2.0};

	/* これらのデータから, 5x5 の行列 m1 を生成し, 印字 
	   この行列は 高橋:数値計算法 p173[2](2) に出てくる.
	 */
	int dim=u1.length;
	TridiagMatrix m1 = new TridiagMatrix(dim,d1,u1,l1);
	System.out.println("m1=");
	System.out.println(m1);




	double u2[] = {1.0,0.0,1.0,0.0,0.0};
	double d2[] = {1.0,2.0,3.0,4.0,5.0};
	double l2[] = {0.0,2.0,3.0,2.0,3.0};
	/* これらのデータから, 5x5 の行列 m2 を生成して印字 */
	// 埋めてね ************************************************

	/* m1+m2 を計算して印字 */
	System.out.println("m1+m2=");
	System.out.println(m1.add(m2));


	TridiagMatrix [] lu;

	/* m1 を LU 分解して, L と U を上の配列にいれ, L と U を印字 */
	lu= m1.decomposeLU();
	System.out.println("L=");
	System.out.println(lu[0]);
	System.out.println("U=");
	System.out.println(lu[1]);

	/* L と U の積を計算して印字 */
	System.out.println("LU=");
	System.out.println(lu[0].multiplyProject(lu[1]));

	// このベクトルを考える
	double v[] = {2.0,-1.0,-1.0,0.0,3.0};
	// 答えを入れるベクトル
	double x[] = new double[v.length];
	// m1 * x ==v を解いて印字
	x=m1.applyInverse(v);
	System.out.println("x=");
	for(int i=0; i < x.length ; i++){
	    System.out.print(x[i] + " ");
	}
	System.out.println();

	/* 検算. v2=m1 * x  して印字する  */	
	double v2[] = {0,0,0,0,0};
	v2=m1.multiply(x);
	System.out.println("m1 * v2=");	
	for(int i=0; i < v2.length ; i++){
	    System.out.print(v2[i] + " ");
	}
	System.out.println();
	    
	return;

    }

}



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

たぶん, 次のような結果が正しいはず.
m1=
4.0 1.0 0.0 0.0 0.0 
2.0 4.0 1.0 0.0 0.0 
0.0 2.0 4.0 1.0 0.0 
0.0 0.0 2.0 4.0 1.0 
0.0 0.0 0.0 2.0 4.0 

m2=
1.0 1.0 0.0 0.0 0.0 
2.0 2.0 0.0 0.0 0.0 
0.0 3.0 3.0 1.0 0.0 
0.0 0.0 2.0 4.0 0.0 
0.0 0.0 0.0 3.0 5.0 

m1+m2=
5.0 2.0 0.0 0.0 0.0 
4.0 6.0 1.0 0.0 0.0 
0.0 5.0 7.0 2.0 0.0 
0.0 0.0 4.0 8.0 1.0 
0.0 0.0 0.0 5.0 9.0 

L=
1.0 0.0 0.0 0.0 0.0 
0.5 1.0 0.0 0.0 0.0 
0.0 0.5714285714285714 1.0 0.0 0.0 
0.0 0.0 0.5833333333333333 1.0 0.0 
0.0 0.0 0.0 0.5853658536585366 1.0 

U=
4.0 1.0 0.0 0.0 0.0 
0.0 3.5 1.0 0.0 0.0 
0.0 0.0 3.428571428571429 1.0 0.0 
0.0 0.0 0.0 3.416666666666667 1.0 
0.0 0.0 0.0 0.0 3.4146341463414633 

LU=
4.0 1.0 0.0 0.0 0.0 
2.0 4.0 1.0 0.0 0.0 
0.0 2.0 4.0 1.0 0.0 
0.0 0.0 2.0 4.0 1.0 
0.0 0.0 0.0 2.0 4.0 

x=
0.6517857142857143 -0.6071428571428571 0.12499999999999997 -0.2857142857142857 0.8928571428571429 
m1 * v2=
2.0 -0.9999999999999998 -1.0 1.1102230246251565E-16 3.0 

課題11_02 拡散方程式の陰解法

課題09_02と同じ偏微分方程式(ノイマン境界条件)を考える. すなわち,区間 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 としよう.

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

/** 陰解法で拡散方程式を解くためのクラス */
public class Sample11_02 extends Sample07_01 {

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

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




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

	Sample11_02 prog= new Sample11_02(); // 変更したよ **************
	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
    }


    /**
       関数 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
    }

    void initialize_matrix(){
	final int dim=xinterval+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 は上で別に扱った)
	 前回は陽解法, 今回は陰解法なので変更が必要. 

	*/

	// 埋めてね **************************************************
	// と言ってもたった1行書けばいいだけ.


	// 時刻 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 Sample11_02.java"
End:
 */

課題09_01 では, 差分の刻みを, x 方向が20, t 方向を 250 としたが, 今回は, t 方向を 50 にしても安定に計算されることを確かめよう. 課題 09_01 をやった人は, 課題 09_01 で t 方向の刻み 50 にするとどうなるか試そう.


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