Skip to content

Instantly share code, notes, and snippets.

@gogotanaka
Created January 18, 2016 15:55
Show Gist options
  • Save gogotanaka/84fb00e7838f5ce28441 to your computer and use it in GitHub Desktop.
Save gogotanaka/84fb00e7838f5ce28441 to your computer and use it in GitHub Desktop.
tridiag.c
/* 三重対角行列を係数行列に持つ連立一次方程式の数値解法(掃き出し法) */
/* 係数行列A=(a_{i,j})はa_{i,j}=2(if i==j),-1(if |i-j|==1),0(otherwise)である。*/
/* これは1次元微分方程式の離散化ではよく見かける形式である */
#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]);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment