誤差

誤差は数値計算のツボ

数値計算のねらいは，できるだけ精確・高速に解を得ることである．誤差と収束性，安定性が数値計算のツボとなる．ここでは，丸め誤差(roundoff error)について見ておこう．打ち切り誤差(truncation error)，収束性と安定性については次章でみる．

整数型と実数型の内部表現

32bitのCPUでは1ワードは4byte(4x8bits)．符号に1bit必要なので，2^(31)-1が1ワードで表現できる最大整数．実数も浮動小数点で表わす場合は1ワードを基本として，以下のような式で表わされる．

Bはbase(基底)で通常は2，Eはbias(下駄)と呼ばれる．64bitのCPUでは1ワードは8bytesになり，realも64bitになるが，現在ではあまりこうは言わないようだ．上記はIEEE規格で現在の主流．

機械精度(Machine epsilon)

小さい数を足したときにマシンがそれを認識できなくなる限界．

#include <stdio.h>

int main(void){
  float e,w;
  e=1; w=1+e;
  while(w>1){
    printf("%-15.10g %-15.10g %-15.10g\n",e,w,w-1);
    e/=2;
    w=1+e;
  }
}

1               2               1              
0.5             1.5             0.5            
0.25            1.25            0.25           
0.125           1.125           0.125          
0.0625          1.0625          0.0625         
0.03125         1.03125         0.03125        
0.015625        1.015625        0.015625       
0.0078125       1.0078125       0.0078125      
0.00390625      1.00390625      0.00390625     
0.001953125     1.001953125     0.001953125    
0.0009765625    1.000976562     0.0009765625   
0.00048828125   1.000488281     0.00048828125  
0.000244140625  1.000244141     0.000244140625 
0.0001220703125 1.00012207      0.0001220703125
6.103515625e-05 1.000061035     6.103515625e-05
3.051757812e-05 1.000030518     3.051757812e-05
1.525878906e-05 1.000015259     1.525878906e-05
7.629394531e-06 1.000007629     7.629394531e-06
3.814697266e-06 1.000003815     3.814697266e-06
1.907348633e-06 1.000001907     1.907348633e-06
9.536743164e-07 1.000000954     9.536743164e-07
4.768371582e-07 1.000000477     4.768371582e-07
2.384185791e-07 1.000000238     2.384185791e-07
1.192092896e-07 1.000000119     1.192092896e-07

桁落ち，情報落ち，積み残し

 0.723657
-0.723649
---------
 0.000008

 72365.7
+    1.23859126
---------
 72366.9

 72365.7
+    0.001
---------
 72365.7

実数，桁落ちのプログラム例

#include <stdio.h>

int main(void){
  float a,b,c;
  double x,y,z;

  a=1.23456789;
  printf("  a= %17.10f\n",a);

  b=100.0;
  c=a+b;
  printf("%20.10f %20.10f %20.10f\n",a,b,c);

  x=(float)1.23456789;
  y=(double)100;
  z=x+y;
  printf("%20.12e %20.12e %20.12e\n",x,y,z);

  x=(double)1.23456789;
  y=(double)100;
  z=x+y;
  printf("%20.12e %20.12e %20.12e\n",x,y,z);

  return 0;
}

#include <stdio.h>
#include <math.h>

int main(void){
  double a,b,c,xp,xm,xpp,xpm;
  a=1;
  b=-10000000;
  c=1;
  xp=(-b+sqrt(pow(b,2)-4.0*a*c))/(2.0*a);
  xm=(-b-sqrt(pow(b,2)-4.0*a*c))/(2.0*a);
  xpp=2*c/(-b-sqrt(pow(b,2)-4.0*a*c));
  xpm=2*c/(-b+sqrt(pow(b,2)-4.0*a*c));
  printf("%25.15f %25.15f \n%25.15f %25.15f\n",xp,xm,xpp,xpm);
	   
  return 0;
}

Fortranの文法

変数型

• complexはreal+ I real，複素数は旧Cには標準ではないが，通常の数値計算ライブラリではheaderなどで独自に定義されていた．C99では標準となる．
• 変数名の先頭は英数，2字目以降は英数字，6字以内．
• プログラムの冒頭で宣言してから使うのが原則．ただし宣言を省略しても，名前がi,j,k,l,m,nで始まる変数は整数として扱われる．
• implicit real(a-h,o-z)という宣言は，前述のi-nで始まる変数が整数，その他で始まる変数が単精度実数を意味することを暗黙に(implicitly)宣言している．double precisionも同様．

注釈，継続行，参照ラベル

Fortranの各行のあたまには通常6文字分の空白をあける．そこになにかが表記されているときは，注釈，継続行，および参照ラベルを意味する．注釈は第1桁目に文字*あるいは文字cがあると注釈行となる．6ケタ目に0以外の文字がある場合，この行は継続行とみなされる．参照ラベルは，すぐ後に出てくる入出力の書式指定や，制御文に使われる．

入出力

scanfとprintfに対応．

write(*,*)あるいはread(*,*)だと書式指定なしでデフォルトへの出力あるいは入力となる．

      read(5,100) a,b,c
      write(6,101) a,b,c
  100 format(3f10.5)
  101 format("Input:",2f10.5," Answer:",f10.5)

算術演算子と関数

• 整数型の変数，定数の演算では結果の小数点以下は切り捨てられる．
• 混合演算：実数型の数と整数型の数との演算も文法違反ではないが，以下のrealあるいはint関数で明示的におこなったほうがよい．

制御文

繰り返しなどの範囲は参照ラベルで指定．

      if (条件文) then
         実行文
      else 
         実行文
      end if

      do 100 i=1,10,2
        実行文
  100 continue 
1ずつ上がるときは最後の2を省略．Cでは
for (i=1; i<=10; i+=2){
    実行文
}

      do 200 while (条件文)
        実行文
  200 end do

関係演算子・論理演算子