tridiag.c

tridiag.c

/*　三重対角行列を係数行列に持つ連立一次方程式の数値解法(掃き出し法)　*/
/*　係数行列A=(a_{i,j})はa_{i,j}=2(if i==j),-1(if |i-j|==1),0(otherwise)である。*/
/*　これは１次元微分方程式の離散化ではよく見かける形式である　*/

#include<stdio.h>

main(){
	int i,j;
	long N=10,NRHS=1,IPIV[N],LDA=N,LDB=N,INFO;
	double DL[N-1],D[N],DU[N-1],b[N],h=1./N,hs=h*h;
	/* 変数の説明 N:変数の数 NRHS:右辺の列数 IPIV[N]:Pivot選択の情報
	LDA:左辺の行数 LDB:右辺の行数 INFO:ルーチンが成功したかどうかを返す
	D:Aの対角成分
	DL:Aの対角部分の一つ下の部分 (DL={a_{2,1},a_{3,2},...,a_{N,N-1}})
	DU:Aの対角部分の一つ上の部分 (DU={a_{1,2},a_{2,3},...,a_{N-1,N}})
	*/
	
	for(i=0;i<N;i++)
		D[i]=2;
	
	for(i=0;i<N-1;i++){
		DL[i]=-1;
		DU[i]=-1;
	}
	
	for(i=0;i<N;i++)
		b[i]=hs;
	
	/*
		実際に三重対角行列方程式を解くサブルーチンは
		dgtsv_(*long,*long,*double,*double,*double,*double,*long,*long)
		である。引数の順に注意
	*/
	
	dgtsv_(&N,&NRHS,DL,D,DU,b,&LDB,&INFO);
	
	/*
	これはf''(x)=1,f(0)=0,f(1)=0という常微分方程式の境界値問題をx[i]=(i+1)*h,(0 \le i \le N-1 ,h=1/(N+1))で離散化したときの
	差分法による数値解法である。方程式の厳密解はf(x)=(e^{1-x}+e^{x})/(1+e)である。
	数値解と厳密解を比較してみよう
	*/
	
	for(i=0;i<N;i++)
		printf("%lf\n",b[i]);
	
}