#!/usr/bin/env python # coding: utf-8 #

Table of Contents

#
#
# #
# 誤差(Error) #
#
#
# file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/error #
# https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python #
# cc by Shigeto R. Nishitani 2017-20 #
# # # # # 打ち切り誤差と丸め誤差(Truncation and round off errors) # # # 数値計算のねらいは，できるだけ精確・高速に解を得ることである．誤差 (精度) と収束性 (安定性,速度) # が数値計算のキモとなる．前回に説明した収束判定条件による誤差は打ち切り誤差 (truncation error) # と呼ばれる．ここでは，誤差のもう一つの代表例である，計算機に特有の丸め誤差 (roundoff error) について見ておこう． # # # # ## 整数型と実数型の内部表現 # # 計算機は一般に無限精度の計算をおこなっているわけではない．CPU で足し算をおこなう以上，一般的な計算においては CPUが扱う一番効率のいい数の大きさが存在する．これが,32bit の CPU では 1 ワード,4byte(4x8bits) である．1ワードで表現できる最大整数は,符号に 1bit 必要なので,2^(31)-1 となる．実数は以下のような仮数部と指数を取る浮動小数点数で表わされる． # # # ### caption:浮動小数点数の内部表現 (IEEE754)． # # # |記号| $s \times f \times B^{e-E}$ | # |:----|:----| # |$s$|sign bit(符号ビット:正負の区別を表す) | # |$f$|fraction portion of the number(仮数部) | # |$e$|biased exponent(指数部) | # |$B$|base(基底) で通常は 2　| # |$E$|bias(下駄)と呼ばれる | # |real(単精度)|$s=1, e=8, f=23, E=127$ | # |double precision(倍精度)|$s=1, e=11, f=52, E=1023$| # # $E$は指数が負となる小数点以下の数を表現するためのもの．演算結果は実際の値から浮動小数点数に変換するための操作「丸め (round-off)」が常に行われる．それに伴って現れる誤差を丸め誤差と呼ぶ． # # # # # # 有効桁数(Significant digits) # # # 1 ワードの整数の最大値とその2進数表示． # ``` python # 2**(4*8-1)-1 # => 2147483647 # ``` # # この整数を2進数で表示するように変換するには，bin(n)を用いて， # ``` python # bin(2**(4*8-1)-1) # => 0b1111111111111111111111111111111 # ``` # となり，31個の1が並んでいることが分かる．1 ワードの整数の最大桁は，$n$ の長さを戻すコマンドlen(n)を使って， # ``` python # len(str(2**(4*8-1)-1)) # => 10 # ``` # となり，たかだか10桁程度であることが分かる．一方，64bit の場合の整数の最大桁． # ``` python # len(str(2**(8*8-1)-1)) # => 19 # ``` # である． # # python では多倍長計算するので，通常のプログラミング言語で起こるintの最大数あたりでの奇妙な振る舞いは示さない． # # ``` python # 2147483647+100 # => 2147483747 # ``` # # 単精度の浮動小数点数は,仮数部2進数23bit,2倍長実数で52bitである．この有効桁数は以下の通り． # ``` python # len(str(2**(23))) # => 7 # len(str(2**(52))) # => 16 # ``` # # pythonでは普通に整数を扱う場合には，1ワードによる制約はなく，大きな数を扱える． # # しかし，numpyなどでarrayとして作る場合には最大数に対する制約が現れる． # 下の例では， # - tmpは普通の変数なので，最大数に対する制約はない．しかし， # - xやyはarrayとして宣言した時点で，int64やfloat64に型が決まってしまう． # - 足し算するとOverflowErrorを吐きだす． # - pythonは中ではCに翻訳しているので，"C long"型よりも大きすぎてダメと根をあげている # # のがわかるでしょう． # In[1]: import numpy tmp = numpy.int(2**(8*8-1)-1) print(tmp) # => 9223372036854775807 print(tmp+100) # => 9223372036854775907 x = numpy.array([0, 0]) print(x.dtype) # => int64 y = numpy.append(x, tmp+1) print(y[2]) # => 9.22337203685e+18 print(y[2]+10) print(y.dtype) x[0] = tmp # OK print(x[0]) x[1] = tmp+1 # # 浮動小数点演算による過ち(FloatingPointArithmetic) # # # 「丸め」にともなって誤差が生じる. CやFortran等の通常のプログラミング言語では「丸める」仕様なのでプログラマーが気をつけなければならない．以下のようなC programでは，予期したのとは違う結果となる． # ```c # //プログラムリスト : 実数のケタ落ち # #include # # 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; # } # ``` # ## 分かっているつもりでも，よくやる間違い． # ``` c # //プログラムリスト : 丸め誤差 # #include # # int main(void){ # float x=77777,y=7,y1,z,z1; # y1=1/y; # z=x/y; # z1=x*y1; # printf("%10.2f %10.2f\n",z,z1); # if (z!=z1){ # printf("z is not equal to z1.\n"); # } # printf("Surprising?? \n\n\n\n\n%10.5f %10.5f\n",z,z1); # return 0; # } # ``` # # # これを避けるには，EPSILONという小さな数字を定義しておいて，値の差の絶対値を求めるfabsを使って # # |\hspace{100mm} | # |:----| # | # | # # # とすべき．このときは数学関数であるfabsを使っているので， # ```maple # > gcc -lm test.c # ``` # とmath libraryを明示的に呼ぶのを忘れないように． # # # In[2]: x=77777 # change this digits y=7 y1=1.0/y z=x/y z1=x*y1 print("%30.20f %30.20f" % (z, z1)) if (z != z1): print("no, different") else: print("yes, identical") # In[3]: e = 10**(-4) x=77777 # change this digits y=7 y1=1.0/y z=x/y z1=x*y1 print("%30.20f %30.20f" % (z, z1)) if ( abs(z-z1) < e ): print("yes, identical") else: print("no, different") # # 機械精度(Machine epsilon) # # # 上の例では，浮動小数点数で計算した場合に小さい数の差を区別することができなくなるということを示している．つまり，小さい数を足したときにその計算機がその差を認識できなくなる限界ということ．これは，昔はCPU固有の精度で，今でも機械精度(Machine epsilon)と呼ばれる．今は，言語の実装によって違ってくるが，以下のようにして求めることができる． # # In[4]: e=1.0 w=1.0+e while(w>1.0): print('%-15.10e %-15.10e %-15.10e' % (e,w,w-1.0)) e = e/2.0 w = 1.0+e # In[5]: from decimal import * getcontext().prec=80 e=Decimal(1.0) w=Decimal(1.0)+e while(w>1.0): print('%-15.10e %-15.10e' % (e,w)) e = e/Decimal(2.0) w = Decimal(1.0)+e # # 桁落ち，情報落ち，積み残し(Cancellation) # # 有効数字がそれぞれ5桁で計算した結果を示せ． # # # #### 桁落ち(Cancellation) # # ```maple # 0.723657 # - 0.723649 # ------------ # # ``` # # # # #### 情報落ち(Loss of Information) # # ```maple # 72365.7 # - 1.23659 # ------------ # # ``` # # # # #### 積み残し # # ```maple # 72365.7 # - 0.001 # ------------ # # ``` # # # # # # ## pythonにおける精度制御演算 # # pythonにはroundという標準的な丸める関数がありますが，使用には注意が必要です． # # https://note.nkmk.me/python-round-decimal-quantize/ # In[6]: print('0.4 =>', round(0.4)) print('0.5 =>', round(0.5)) print('0.6 =>', round(0.6)) # In[7]: print('0.5 =>', round(0.5)) print('1.5 =>', round(1.5)) print('2.5 =>', round(2.5)) print('3.5 =>', round(3.5)) # In[8]: print('0.05 =>', round(0.05, 1)) print('0.15 =>', round(0.15, 1)) print('0.25 =>', round(0.25, 1)) print('0.35 =>', round(0.35, 1)) print('0.45 =>', round(0.45, 1)) # # ``` # 注釈 浮動小数点数に対する round() の振る舞いは意外なものかもしれません: 例えば、 round(2.675, 2) は予想通りの 2.68 ではなく 2.67 を与えます。これはバグではありません: これはほとんどの小数が浮動小数点数で正確に表せないことの結果です。 # # 2. 組み込み関数 round() — Python 3.6.3 ドキュメント # ``` # さらに，pythonにはdeimcalという，10進固定及び浮動小数点数演算用のライブラリがあります．この中のquantize関数で直感的な四捨五入が実現できます．こちらの方が良さそう． # In[1]: from decimal import * def pretty_p(result,a,b,operator): print('context.prec:{}'.format(getcontext().prec)) print(' %20.14f' % (a)) print( '%1s%20.14f' % (operator, b)) print('-----------') print( ' %20.14f' % (result)) n = 10 getcontext().prec = n a=Decimal('0.723657').quantize(Decimal(10) ** -n) b=Decimal('0.723649').quantize(Decimal(10) ** -n) pretty_p(a-b,a,b,'-') a=Decimal('0.723657').quantize(Decimal(10) ** -n)*100000 b=Decimal('0.123659').quantize(Decimal(10) ** -n)*10 pretty_p(a-b,a,b,'-') a=Decimal('0.723657').quantize(Decimal(10) ** -n)*100000 b=Decimal('0.1').quantize(Decimal(10) ** -n)/100 pretty_p(a+b,a,b,'+') # In[2]: n = 5 getcontext().prec = n a=Decimal('0.723657').quantize(Decimal(10) ** -n) b=Decimal('0.723649').quantize(Decimal(10) ** -n) pretty_p(a-b,a,b,'-') a=Decimal('0.723657').quantize(Decimal(10) ** -n)*100000 b=Decimal('0.123659').quantize(Decimal(10) ** -n)*10 pretty_p(a-b,a,b,'-') a=Decimal('0.723657').quantize(Decimal(10) ** -n)*100000 b=Decimal('0.1').quantize(Decimal(10) ** -n)/100 pretty_p(a+b,a,b,'+') # # 課題 # # ## 有効数字 # 1. 大きな数どおしのわずかな差は，丸め誤差にとくに影響を受ける. # 23.173-23.094 を有効数字がそれぞれ5桁，4桁，3桁，2桁で計算した結果を示せ. # 同様に，0.81321/(23.173-23.094) を有効数字がそれぞれ5桁，4桁，3桁，2桁で計算した結果を示せ. # 1. 10 進数10桁および3桁の有効桁数をもった計算機になったつもりで，以下の条件で預金を求める計算をおこなえ． # 1. 元本を10000万円とする # 1. 利息0.3%とする # 1. 複利計算で10年でいくらになるか． # ## 2次方程式解の公式の罠 # # 1. 係数を a = 1, b = 10000000, c = 1としたときに, 通常の解の公式を使った解と, 解と係数の関係(下記の記述を参照）を使った解とを出力するプログラムをpythonで作成すると以下の通りとなる．解の有効数字が2種類の計算方法の違いでどう違うか，いくつかの精度で実行させた結果を使って解説せよ． # # 2次方程式 $ax^2+bx+c=0$の # 係数$a, b, c$が特殊な値をもつ場合，通常の解の公式 # # # $$ # x = \frac {-b \pm \sqrt{{b}^{2}-4ac}}{2a} # $$ # にしたがって計算するとケタ落ちによる間違った答えを出す．その特殊な値とは # # # $$ # \sqrt{{b}^{2}-4ac} \approx |b| # $$ # となる場合である． # # ケタ落ちを防ぐには, $b > 0$の場合は, # # # $$ # x_1 = \frac {-b - \sqrt{{b}^{2}-4ac}}{2a} # $$ # として，ケタ落ちを起こさずに求め, この解を使って, 解と係数の関係より # # # $$ # x_2 = \frac {c}{a\, x_1} \,\ldots(1) # $$ # で求める．$b < 0$ の場合は，解の公式のたし算の方 # $$ # x_1 = \frac {-b + \sqrt{{b}^{2}-4ac}}{2a} # $$ # を使って同様に求める． # # In[3]: import numpy as np def solve_by_formula(a,b,c): x0 = (-b + np.sqrt(b**2-4*a*c))/(2*a) x1 = (-b - np.sqrt(b**2-4*a*c))/(2*a) return (x0, x1) def solve_precise(a,b,c): x0 = (-b - np.sqrt(b**2-4*a*c))/(2*a) x1 = c/(a*x0) return (x0, x1) print(solve_by_formula(1,10000000,1)) print(solve_precise(1,10000000,1)) # In[4]: from numpy import sqrt def solve_normal_formula(a,b,c): x0=(-b-sqrt(b**2-4*a*c))/(2*a) x1=(-b+sqrt(b**2-4*a*c))/(2*a) return (x0,x1) def solve_precise_formula(a,b,c): x0=(-b-sqrt(b**2-4*a*c))/(2*a) x1=c/(a*x0) return (x0,x1) print(solve_normal_formula(1,100000000,1)) print(solve_precise_formula(1,100000000,1)) # In[14]: from decimal import * print("prec=40") getcontext().prec = 40 print(solve_normal_formula(Decimal('1'), Decimal('100000000'), Decimal('1'))) print(solve_precise_formula(Decimal('1'), Decimal('100000000'), Decimal('1'))) print("\nprec=20") getcontext().prec = 20 print(solve_normal_formula(Decimal('1'), Decimal('100000000'), Decimal('1'))) print(solve_precise_formula(Decimal('1'), Decimal('100000000'), Decimal('1'))) # ### (1)式の導出 # # (1)式の導出を示しておく．解と係数の関係を導びくと # \begin{eqnarray} # (x-x_1)(x-x_2) &=& x^2 &-&(x_1+x_2)x &+&x_1 x_2　\\ # &=& \frac{a}{a}x^2 &+& \frac{b}{a}x &+& \frac{c}{a} # \end{eqnarray} # となる．(1)式は最後の定数項より， # $$ # \frac{c}{a} = x_1 x_2 # $$ # から導かれる． # In[ ]: