今週のプログラムその4
今週は行列計算の応用としてガウスの消去法のプログラムと、
関数副プログラムの簡単な使い方として非線形方程式の解を求める二分法の
プログラム。
ただし、消去法のプログラムはピボット選択をしていないので、このままでは
問題が生じることもあります。あくまでもプログラミングの練習ということで。
ポイント
- 実数(浮動小数点)演算の精度は計算機の語長(数字ひとつが何ビットで表されるか)で決まる。
一語32ビットの機械の場合、実数普通10進数で7桁程度の精度がある。
しかし、計算によってはこれでは精度が不足する場合があるので、数字ひとつを2語で
表す(32ビットの機械なら64ビットで)ように指定する倍精度実数を使うことができる。
たとえば、数aと配列b(要素数10)を倍精度数として宣言するには
double a,b[10]; (Cの場合)
real*8 a,b(10) (FORTRANの場合)
などとする。
昔と違って今の計算機は倍精度実数の計算だからといって特に計算時間がかかるわけではないので、科学計算の場合は常に倍精度を指定すべきだろう。
- 数字を0.000000000123などと書いていては大変なので、たとえば1.23*10-10は1.23e-10という書きかたができる
- 二分法のプログラムでは、解を得たい関数について変数の値を変えながら何度も計算する。
このような場合、関数を副プログラムとしてプログラム本体からわけることにより、
関数に名前をつけ、あたかもsin(x)やexp(x)などのあらかじめ用意された
関数と同様に扱うことができる。
たとえば、一変数関数
x3 + 3*x2 + 5*x + 1という関数にfuncという名前をつける
にはFORTRANでは
function func(x)
real*8 func 関数が倍精度の値を返すことを宣言
real*8 x 変数が倍精度数であることを宣言
func= x**3 + 3*x**2 + 5*x + 1 この関数が返す値
end
とする。ここで変数の名前xは、倍精度の変数を関数の中でxと呼ぶという宣言で
あって、実際に関数を使う際に主プログラム中で同じxという名前を使う必要は
ない。また、主プログラム中でxという変数を使っていても、それは関数の計算に
影響を与えない。その意味で、関数副プログラム中に現われる変数名は局所的である。
同様にCでは
double func(double x) 倍精度変数に対して倍精度数を返す関数を宣言
{
return pow(x,3) + 3*x*x + 5*x + 1; return 以下が主プログラムに返す値
}
とすればよい(powはべき乗を計算する関数)。
ただし、Cでは、使用する関数の名前と種類をmain()の前にあらかじめ宣言しておく
必要がある(関数プロトタイプ宣言)。今の例では
double func(double);
がプロトタイプ宣言。
- Cには大きくわけてK&RとANSIの二種類があり、古いCコンパイラでは何も指定
しないとK&Rとみなされる場合がある。
例題はすべてANSI準拠なので(といっても、前回までのプログラムはどちらでも
同じ)、コンパイラによってはANSI Cであることを宣言する必要がある。
例えば
cc -Aa -o prog prog.c -lm
など(ここで -Aa がANSI Cを宣言するオプション)。
コンパイル時に警告メッセージが出ないなら気にする必要はない。
- Cプログラムでは
printf("solution is between %17.12f and %17.12f\n",x1,x2)
のように、%fではなく%12.7fと書いたところがある。これは出力のしかたを指定するもので、今の例では実数ひとつを17桁分の幅の中に書き、小数点以下12桁を表示するという意味。17桁には小数点やマイナス記号も含まれる。
この手の書式指定は自習してもらうつもりだったが、%fではあまりにも出力される桁数が少ない場合があるので、例外的に使用した。FORTRANでも同じことはできるが、こちらは(*,*)でも充分な桁数を出力する場合が多いようなので、あえて解説しない。
発展
- Gaussの消去法にピボット選択を組み込んでみる(ピボット選択については、そのへんの解説書を見てください)
- 二分法のプログラムでは、解への収束を相対誤差で判定している。しかし、この方法は解が0に近い場合にはうまくいかない。そこで解が0に近い場合を扱えるように工夫してみる。
- いろいろな関数を二分法で解いてみる。
まずはC版
program 10
..................
#include<stdio.h>
#define NSIZ 3 プログラム中のNSIZという文字をすべて3に置き換える
/* 行列演算の応用
Gauss の消去法(ピボット選択なし)
a x = b の解 x を求める */
main()
{
double a[NSIZ][NSIZ],b[NSIZ],x[NSIZ]; 倍精度数の配列を宣言
double awork[NSIZ][NSIZ],bwork[NSIZ],check[NSIZ];
double dd,xx; 倍精度の変数を宣言
int i,j,k;
/* ファイルから行列aに値を読み込む */
for( i=0; i<NSIZ; ++i){
for( j=0; j<NSIZ; ++j){
scanf("%lf",&a[i][j]); 倍精度数に読み込むときは"%f"ではなく"%lf"とする(lはlongの意味)
}
}
/* ファイルからベクトルbに値を読み込む */
for( i=0; i<NSIZ; ++i){
scanf("%lf",&b[i]); ここも倍精度なので"%lf"
}
/* 作業用の行列aworkとベクトルbworkにaとbをコピー */
for( i=0; i<NSIZ; ++i){
for( j=0; j<NSIZ; ++j){
awork[i][j]=a[i][j];
}
}
for( i=0; i<NSIZ; ++i){
bwork[i]=b[i];
}
/* 前進消去 */
for( i=1; i<NSIZ; ++i){ iは1からNSIZ-1まで。行列は0行目から始まっていることに注意
for( j=i; j<NSIZ; ++j){ jはiからNSIZ-1まで
dd=awork[j][i-1]/awork[i-1][i-1]; このへんは前進消去の意味が分かっていれば理解できるはず
awork[j][i-1]=0;
bwork[j] -= dd*bwork[i-1];
for( k=i; k<NSIZ; ++k){ kはiからNSIZ-1まで
awork[j][k] -= dd*awork[i-1][k];
}
}
}
/* 消去結果の出力(右上三角行列になったことを確認) */
printf("result of forward elimination\n");
printf("matrix awork =\n");
for( i=0; i<NSIZ; ++i){
for( j=0; j<NSIZ; ++j){
printf(" %f",awork[i][j]);
}
printf("\n");
}
printf("vector bwork =\n");
for( i=0; i<NSIZ; ++i){
printf(" %f",bwork[i]);
}
printf("\n");
/* 後退代入 */
for( i=NSIZ-1; i>=0; --i){ iの値をNSIZ-1から0まで、1ずつ減らしてゆく
xx=0;
for( j=i+1; j<NSIZ; ++j){
xx += x[j]*awork[i][j]; このへんは後退代入の意味がわかっていれば理解できるはず
}
x[i] = (bwork[i]-xx)/awork[i][i];
}
/* 代入結果の出力 */
printf("result of backward substitution\n");
printf("vector x =\n");
for( i=0; i<NSIZ; ++i){
printf(" %f",x[i]);
}
printf("\n");
/* 結果のチェック( a x を計算) */
for( i=0; i<NSIZ; ++i){
check[i] = 0;
for( j=0; j<NSIZ; ++j){
check[i] += a[i][j]*x[j]; 行列とベクトルの掛け算
}
}
/* チェック結果の出力(a x - b = 0であることを確認) */
printf("result of consistency check\n");
printf("a x - b =\n");
for( i=0; i<NSIZ; ++i){
printf(" %f",check[i]-b[i]); check[i]-b[i]を出力。このようにprintfの中で計算してもよい
}
printf("\n");
}
..................
program 11
..................
#include<stdio.h>
#include<math.h> 数学関数を使うおまじない
#define NMAX 30 プログラム中のNMAXをすべて30に置き換え
/* 関数サブプログラムの利用
2分法による非線形方程式の解 */
double func(double); 倍精度変数を入れると倍精度数を返す関数funcを使うという宣言(関数プロトタイプ宣言)
main() ここから主プログラム
{
double x1,x2,xtemp; 倍精度変数x1,x2,xtempの宣言
int irepeat; 整数変数irepeatの宣言
/* 解探索の範囲を指定 */
scanf("%lf%lf",&x1,&x2); ふたつの倍精度数を読み込む
/* 2分法を最大NMAX回くりかえす */
for(irepeat=1; irepeat<=NMAX; ++irepeat){
/* x1とx2の差が与えた精度以内なら終了 */
if( fabs((x1-x2)/x1) <= 1e-10){ x1とx2の相対差が1e-10以下なら収束とみなす
printf("solution is between %17.12f and %17.12f\n",x1,x2);
return; プログラムの実行を終了
}
/* func(x1)とfunc(x2)の符号が逆なら範囲を二分する */
else if(func(x1)*func(x2) < 0){ func(x1)とfunc(x2)が逆符号なら
xtemp = (x1+x2)/2; x1とx2の中点の値をxtempに
if(func(x1)*func(xtemp) < 0){ func(x1)とfunc(xtemp)が逆符号なら
x2=xtemp; x1はそのままで、xtempを新しいx2にする
}
else{ func(x1)とfunc(xtemp)が同符号なら
x1=xtemp; x2はそのままで、xtempを新しいx1にする
}
printf("%d step x1= %17.12f x2= %17.12f\n",irepeat,x1,x2); x1,x2の値を出力
}
/* func(x1)とfunc(x2)が同符号なら、この範囲に解がないと判断して終了 */
else{
printf("no solution is found between %17.12f and %17.12f\n",x1,x2);
return; プログラムの実行を終了
}
} ここまでがforで指定した繰り返し範囲
/* NMAX回のくりかえしで解に収束しなかったときは終了 */
printf("not converged after %d steps\n",NMAX);
printf("current values of x are %17.12f and %17.12f\n",x1,x2);
} ここまでが主プログラムの範囲("main(){"からこの"}"まで
/* 解を得たい関数 */
double func(double x) ここからが関数副プログラム。関数の名前はfunc、ひとつの倍精度変数xに対して倍精度の値を返す
{
return pow(x,3) + 3*x*x + 5*x + 1; return 以下が主プログラムに返す値
} 関数副プログラムの終わり
FORTRAN版
program 10
..................
C 行列演算の応用
C Gauss の消去法(ピボット選択なし)
C a x = b の解 x を求める
program main
parameter(nn=3) プログラム中のnnの値をすべて3に
real*8 a(nn,nn),b(nn),x(nn) 倍精度配列の宣言
real*8 awork(nn,nn),bwork(nn),check(nn)
real*8 dd,xx 倍精度変数の宣言
C ファイルから行列aに値を読み込む
do i=1,nn
read(*,*)(a(i,j),j=1,nn)
enddo
C ファイルからベクトルbに行列を読み込む
read(*,*)(b(i),i=1,nn)
C 作業用の行列aworkとベクトルbworkにaとbをコピー
do i=1,nn
do j=1,nn
awork(i,j)=a(i,j)
enddo
enddo
do i=1,nn
bwork(i)=b(i)
enddo
C 前進消去
do i=2,nn 第二行目からnn行目まで
do j=i,nn 第i行目からnn行目まで
dd=awork(j,i-1)/awork(i-1,i-1) このへんは前進消去の意味が分かっていれば理解できるはず
awork(j,i-1)=0
bwork(j)=bwork(j)-dd*bwork(i-1)
do k=i,nn 第i行目からnn行目まで
awork(j,k)=awork(j,k)-dd*awork(i-1,k)
enddo
enddo
enddo
C 消去結果の出力(右上三角行列になったことを確認)
write(*,*)"result of forward elimination"
write(*,*)"matrix awork ="
do i=1,nn
write(*,*)(awork(i,j),j=1,nn)
enddo
write(*,*)"vector bwork ="
write(*,*)(bwork(i),i=1,nn)
C 後退代入
do i=nn,1,-1 iの値をnnから1まで、1ずつ減らしてゆく
xx=0
do j=i+1,nn
xx=xx+x(j)*awork(i,j) このへんは後退代入の意味が分かっていれば理解できるはず
enddo
x(i)=(bwork(i)-xx)/awork(i,i)
enddo
C 代入結果の出力
write(*,*)"result of backward substitution"
write(*,*)"vector x ="
write(*,*)(x(i),i=1,nn)
C 結果のチェック( a x を計算)
do i=1,nn
check(i)=0
do j=1,nn
check(i)=check(i)+a(i,j)*x(j) 行列とベクトルの掛け算
enddo
enddo
C チェック結果の出力(ax-b=0であることを確認)
write(*,*)"result of consistency check"
write(*,*)"a x - b ="
write(*,*)(check(i)-b(i),i=1,nn) check(i)-b(i)を出力。このようにwriteの中で計算してもよい
end
..................
program 11
..................
C 関数サブプログラムの利用
C 2分法による非線形方程式の解
program main
parameter(nmax=30) プログラム中のnmaxの値をすべて30に
real*8 x1,x2,xtemp 倍精度変数の宣言
real*8 func 倍精度関数funcを使う宣言
C 解探索の範囲を指定
read(*,*)x1,x2
C 2分法を最大nmax回くりかえす
do irepeat = 1,nmax
C x1とx2の差が与えた精度以内なら終了
if(abs((x1-x2)/x1).le. 1e-10)then x1とx2の相対差が1e-10以下なら収束とみなす
write(*,*)"solution is between ",x1," and ",x2
stop プログラムの実行を終了
C func(x1)とfunc(x2)の符号が逆なら範囲を二分する
else if(func(x1)*func(x2).lt.0)then func(x1)とfunc(x2)が逆符号なら
xtemp = (x1+x2)/2 x1とx2の中点の値をxtempに
if(func(x1)*func(xtemp).lt.0)then func(x1)とfunc(xtemp)が逆符号なら
x2=xtemp x1はそのままで、xtempを新しいx2にする
else func(x1)とfunc(xtemp)が同符号なら
x1=xtemp x2はそのままで、xtempを新しいx1にする
endif
write(*,*)irepeat," step x1=",x1," x2=",x2
C func(x1)とfunc(x2)が同符号なら、この範囲に解がないと判断して終了
else
write(*,*)"no solution is found between ",x1," and ",x2
stop プログラムの実行を終了
endif
enddo ここまでがdoで指定した繰り返し範囲
C nmax回のくりかえしで解に収束しなかったときは終了
write(*,*)"not converged after ",nmax," steps"
write(*,*)"current values of x are ",x1," and ",x2
end 主プログラムの終わり
C 解を得たい関数
function func(x) ここから関数副プログラム。関数名はfunc、ひとつの変数xに対して値を返すことを宣言
real*8 func 関数funcは倍精度の値を返す
real*8 x 変数xは倍精度であることを宣言
func= x**3 + 3*x**2 + 5*x + 1 xを使ってfuncの値を計算。これが主プログラムに返される値
end 関数副プログラムの終わり