情報科学 第６回　数値解析(1).

Presentation on theme: "情報科学 第６回　数値解析(1)."— Presentation transcript:

1 情報科学 第６回　数値解析(1)

2 数値計算とは？ 代数演算 数値計算 等式を一定の規則のもとに変形して解析的に求める 反復計算によって”近似的に”解を求める
解析的に求めることが困難な方程式を解く場合に用いる 例）　物理学、化学、天文学などにおける数値積分 　　　や微分方程式の解法など

3 数値解析(1)トピック 数値計算における誤差 数値積分 微分方程式の解法

4 1.数値計算における誤差 実数データ表現 丸め誤差 桁落ち 情報落ち 打切り誤差 浮動小数点表現
0.1の2進数表記・ 単精度表現・倍精度表現 桁落ち 情報落ち 打切り誤差 Taylor展開・級数展開の高次項の影響

5 この表記は、例えば、1.0101…010を２進数として解釈するという意味である
実数データ表現 計算機内の数値データは，有限長のビット列で表現される 大きな数や小さな数を表現するためには， 符号部と指数部，仮数部をもった実数表現　　(浮動小数点表現)が利用される IEEE 754　実数表現に関する規格 符号 仮数 指数 この表記は、例えば、1.0101…010を２進数として解釈するという意味である バイアス（指数が負の値をとれるようにするもの）

6 単精度と倍精度 単精度 32ビットで実数を表現 倍精度 64ビットで実数を表現 IEEE 754 実数表現に関する規格
単精度　32ビットで実数を表現 倍精度　64ビットで実数を表現 指数部と仮数部がすべて０の実数は０（±０）を意味する このほかにも±無限大などの表現が可能 整数の大小比較回路がそのまま使える

7 丸め誤差 10進数を2進数で表現する場合、有限桁では近似的にしか表せないことに起因する誤差 例）
0.1 (10) = … = … (2) x 20 = … (2) x 2-4 0.1(10)は有限桁では近似的にしか表現できない (10進数の小数で1/3や1/7を表す場合に相当） 1/3 = … 1/7 = …

8 丸め誤差の例 誤差 # irbで以下を試してみる (0.1*3==0.3) # 上記と同じことを実現するのに以下の書き方もある
Ruby # irbで以下を試してみる (0.1*3==0.3) # 上記と同じことを実現するのに以下の書き方もある (0.1*3-0.3).zero? これは0.1の3倍が0.3に等しいかどうかを真偽で評価するための表現で 数学的には真(true)という結果になるはずであるが、実際は偽(false)となる 符号部 指数部 仮数部 0捨1入 倍精度表現では ["0", " ", " … "] 0捨1入 ["0", " ", " … "] 誤差 ["0", " ", " … "] 丸め誤差は # （検証）irbで以下を試してみよ 0.1*3-0.3 2**-54 0が52個

9

10 単精度表現 仮数部23ビット 指数部8ビット 仮数の表現範囲は有効桁数は24ビット これは10進数で7桁程度の精度しかないことを意味する
1.00…0 (2) =1 から 1.11…1(2) = ≈2-1.2x10-7 ≈ 2 指数の表現範囲は 1.2 x から1.7 x 1038 (2) -127 =2-126 ≈ 1.2 x 10-38 (2) -127 =2127 ≈ 1.7 x 1038 (2) と、 (2) は 特別な意味を持つので 除外してある 最大値の限界は　 ≈ 2 x 2127=2128 ≈ 3.4 x 1038 となる

11 倍精度表現 仮数部52ビット 指数部11ビット 仮数の表現範囲は有効桁数は53ビット これは10進数で16桁程度の精度
1.00…0(2) =1 から 1.11…1(2) = ≈ 2-2.2x10-16 ≈ 2 指数の表現範囲は 2.2 x から9.0 x 10307 2000…01 (2) = ≈ 2.2 x 2111…10 (2) =21023 ≈ 9.0 x 10307 2000…000 (2) と 2111…111 (2) は 特別な意味を持つので 除外してある 最大値の限界は　2 x 21023=21024 ≈ 1.7 x となる

12 桁落ち 計算機の数値計算では誤差が拡大する可能性がある 桁落ち誤差 情報落ち誤差 桁落ち誤差と情報落ち誤差がある 近い数値間の減算による誤差
ほぼ同じような数値の差をとると有効桁数が減少する 情報落ち誤差 大きさの異なる数値間の加減算による誤差 大きさの異なる数値の加減算では、小さな数値は大きな数値の有効桁範囲外になり無視されてしまう

13 桁落ち誤差の例 桁落ち誤差 打切り誤差 def f(x) return Math.log(x); # f(x)=ln(x) end
Ruby h を十分小さく取れば微係数は、 　　　　　　　　　で近似できる 桁落ち誤差 h を変えてx=1における微係数の誤差を調べる 打切り誤差 def f(x) return Math.log(x); # f(x)=ln(x) end x=1.0; 15.times { |i| h =0.1**i #h=0.1, 0.01, …, 1e-15 df=(f(x+h)-f(x))/h err=(df-1.0).abs; print "h=",h,"\t f'(",x,")=",df," err=",err,"\n" } 　　　を小さくしていくと、 ある程度までは打切り誤差が 減少していくが、より微小になると 　　　　　と　　　　　の差も微小となり 有効桁数が減少する その結果、桁落ち誤差が大きくなる ※打切り誤差については 後ほど述べる

14 情報落ち誤差の例1 1をn分割しn回足し算して1になるかどうか検証する err2.java 実行結果 丸め誤差 情報落ち誤差 Java
public class err2 { public static void main(String[] args) { int n = 1; for (int i = 0; i < 10; i++) { float sum = 0.0f; float d = 1.0f / n; for (int j = 0; j < n; j++) sum += d; System.out.println("n = " + n + "\t d = " + d + "\t sum = "+sum); n *= 10; }}} err2.java n = 1 d = 1.0 sum = 1.0 n = 10 d = 0.1 sum = n = 100 d = sum = n = d = sum = n = d = 1.0E-4 sum = n = d = 1.0E-5 sum = n = d = 1.0E-6 sum = n = d = 1.0E-7 sum = n = d = 1.0E-8 sum = 0.25 n = d = 1.0E-9 sum = 実行結果 丸め誤差 情報落ち誤差

15 情報落ち誤差の例2 1をn分割しn回足し算して1になるかどうか検証する err2.rb 7.times{ |m|
Ruby 1をn分割しn回足し算して1になるかどうか検証する 7.times{ |m| n=10**m # n=1,10,100,…, sum=0.0 d=1.0/n n.times{ sum=sum+d } 　# 中はsum+=dでもいい print "n=",n,"\t d=",d,"\t sum=",sum,"\n" } err2.rb Rubyでは実数の表現に倍精度浮動小数点数 を用いているため、情報落ち誤差が生じにくい Javaなどで単精度浮動小数点数を用いると、 情報落ち誤差が生じる場合がある 実行結果 n= d=1.0 sum=1.0 n= d=0.1 sum=1.0 n= d=0.01 sum=1.0 n= d= sum=1.0 n= d= sum= n= d=1.0e sum= n= d=1.0e sum= 丸め誤差

16 打切り誤差 ある関数の値を無限級数を用いて数値計算する場合、有限の項数で打切って近似することにより生じる誤差 （たとえ，実数が無限精度で表現できても生じる） 例）指数関数の原点の周りのテイラー展開 無限回の計算を行うことはできないので、有限回（n回）で 打切って近似を行う スライド12では 微係数を求める場合の 打切り誤差の例を紹介した 誤差の主要成分は

17 2. 数値積分 台形公式 シンプソン公式

18 台形公式 積分を区分線形(piecewise linear)で近似する方法 全積分区間をn等分し各部分区間
　の積分を台形の面積として計算し総和をとる

19 台形公式（続き） ある関数 の から までの 積分の近似値を求めるには、 各部分区間 における台形の総和で近似する

20 シンプソン公式 積分を1次式ではなく2次式で近似する方法 台形公式(1次式近似) シンプソン公式(2次式近似)

21 ラグランジュ補間 ある関数　　　を3点 で補間する2次関数

22 シンプソン公式（続き） 部分区間 とし 区間 の部分の積分近似値は となるため、 から までの 積分近似値の総和は

23 3．常微分方程式の解法 Euler Method Runge-Kutta (Gill) Method

24 常微分方程式の数値解法 高階(n階)の常微分方程式は一般に n組の連立1階常微分方程式に帰着される 1階の常微分方程式の解法が基本となる
(1階)常微分方程式の数値解法にはTaylor展開を利用する場合が多い 常微分方程式 を初期条件 で解く場合 初期条件近傍の解 はTaylor展開で求まる

25 Euler法 常微分方程式 を初期条件 で解く場合 を起点としてその近傍 を以下のように逐次求める (Eulerの公式) 任意の値
はこの漸化式を数値計算することにより 数列として求めることができる 問題点：精度が良くない・誤差が蓄積する

26 簡単な例 def euler(n) y = 1.0 h = 1.0/n for i in 0..n-1
Ruby 簡単な例 def euler(n) y = 1.0 h = 1.0/n for i in 0..n-1 y = y + (1.0/(2*y))*h end y

27 Taylor展開・数値解法の問題 アイデア：Euler法（１次）よりも高次の微分係数を用いる方法が考えられる
問題点：計算機では高次の導関数を求める数式処理が一般に容易ではない 例） ではどうするか？

28 Taylor級数（下の式）と一致するように
Runge-Kutta法 別の方法でTaylor級数と同値なものを計算する方法を考える アイデア：高次導関数を求めるのではなく、高次導関数の値を差分で近似する ※高次導関数が 現れていないこと に注意
上の式で の値が、 の2次の項までは、 Taylor級数（下の式）と一致するように を決めたい

29 Runge-Kutta法（続き） に注意すると そこで と の両式が に関係なく等しくなるためには 従って (Taylor展開 2次まで）
の項は に含まれる そこで と (Taylor展開 2次まで） の両式が に関係なく等しくなるためには 従って

30 Runge-Kutta法（続き） 従って、以下の方法での数値解法は、 の2次の項まではTaylor級数と一致する
同様の手法用いることにより、 高次のRunge-Kutta法を導くこともできる

31 簡単な例 def runge_kutta(n) y = 1.0 h = 1.0/n for i in 0..n-1 k1 = ...
Ruby 簡単な例 def runge_kutta(n) y = 1.0 h = 1.0/n for i in 0..n-1 k1 = ... k2 = ... y = y + (k1+k2)/2.0 end y

32 収束に関するコメント Euler Runge-Kutta hの項まで近似している 各ステップの誤差はh2に比例 ステップ数は1/hに比例

33 Runge-Kutta法（４次） Runge-Kutta法は原理的に高次のものを用いることができる
4次の場合は、 の項を計算するが、　このうち　　　　　　と　　にかかる係数に任意性がある（Taylor級数と一致するという要件では一意に決まらない） 4次の Runge-Kutta法 （ここでは結果だけ示すが、 自分で構築して見よ）

34 Runge-Kutta-Gill法（４次）
何故このような係数を用いるのか？ 計算の方法で各ステップで３つの変数の値を記憶するだけで良い

35 Gill法の特徴 通常Runge-Kuttaでは、 を計算するのに5つの変数 を記憶する必要がある
記憶する変数 step1 step2 step3 step4

36 同様に高階の微分方程式も連立1次微分方程式に帰着できる
連立微分方程式 物理現象の多くは2階の微分方程式で記述できる 一般に 2階微分方程式 は 変数変換をすることにより、連立１次微分方程式 （各変数に関する１次微分方程式の集合）に帰着できる 例）調和振動子の運動方程式 において とおくと となる 同様に高階の微分方程式も連立1次微分方程式に帰着できる

37 連立微分方程式の数値解法 連立微分方程式の場合も同様にRunge-Kutta法を適用できる

38 Runge-Kutta法の例 平面上で質点Mの周りを運動する質点mの軌跡 質点Mの座標系で質点mの座標を とする このとき運動方程式は、
Ruby 平面上で質点Mの周りを運動する質点mの軌跡 質点Mの座標系で質点mの座標を　　　　　とする このとき運動方程式は、 となる ここで変数変換を行い、 ４連１次微分方程式が得られるので、これを4次のRunge-Kutta法で解いてみよう

39 Runge-Kutta法の例 （続き） メインのプログラム 配列u と時間変数ｔとから各微係数を求めて、新しい配列new_uを返す関数
Ruby 　メインのプログラム h=0.1　　#時間変数ｔの差分 t=0　　　 # 時間変数の初期値 u=[1,0,0,1]　# 配列uで4変数u0,u1,u2,u3を定義（初期条件を代入） 100.times{　#計算のステップを１００回繰り返す u=rk(t,u,h,&func)　#Runge-Kutta法の計算を1ステップ分行う t+=h　#時間変数をhだけ増加させる print "t= ", t0, " u[0]= ",u[0]," u[1]= ",u[1]," u[2]= ",u[2]," u[3]= ",u[3]," \n" } 配列u と時間変数ｔとから各微係数を求めて、新しい配列new_uを返す関数 Runge-Kutta法の関数へ引数として渡される func=lambda {|t,u| new_u=[0,0,0,0] r=(u[0]**2+u[1]**2)**1.5 new_u[0]=u[2]; new_u[1]=u[3]; new_u[2]=-u[0]/r; new_u[3]=-u[1]/r; return new_u } ※簡単のためここでは　　　　としている

40 Runge-Kutta法の例（続き） （補足）ここではプログラムを簡潔に書くため、種々のテクニックを使用している
Ruby （補足）ここではプログラムを簡潔に書くため、種々のテクニックを使用している 興味のある者はRubyの参考書を見ながらプログラムを読んでみよう Runge-Kutta法の関数(スライド32の計算を4つの変数に拡張して行っている） def rk(t,u,h,&f) k1=muladd([[f.call(t,u),h]]) k2=muladd([[f.call(t+h/2.0,muladd([[u,1.0],[k1,1.0/2.0]])),h]]) k3=muladd([[f.call(t+h/2.0,muladd([[u,1],[k2,1.0/2.0]])),h]]) k4=muladd([[f.call(t+h,muladd([[u,1.0],[k3,1.0]])),h]]) new_u=muladd([[u,1],[k1,1.0/6.0],[k2,1.0/3.0],[k3,1.0/3.0],[k4,1.0/6.0]]) return new_u end uは4つの変数を含む配列であるが、実数の1変数と見なすと　 1変数のRunge-Kutta法の計算に等しい 上記Runge-Kutta法の関数で使う配列を線形結合する関数 配列aと実数kのペア(a,k)のリストpairsを引数とする 各要素(a,k)に対し、配列aの要素をk倍したものを 加算し、結果の配列を返す （ベクトルの線形結合と思えば良い） def muladd(pairs) r=Array.new(pairs.first.first) r.fill(0) pairs.each{|a,k| r.size.times{|i| r[i]+=k*a[i]}} return r end

41 Runge-Kutta法の例（続き） 実行結果 であったから は質点の軌跡を表す をgnuplotやexcelなどを用いて
Ruby 実行結果 t= 0.1 u[0]= u[1]= u[2]= u[3]= t= 0.2 u[0]= u[1]= u[2]= u[3]= t= 0.3 u[0]= u[1]= u[2]= u[3]= t= 0.4 u[0]= u[1]= u[2]= u[3]= t= 0.5 u[0]= u[1]= u[2]= u[3]= であったから は質点の軌跡を表す 　　　　　をgnuplotやexcelなどを用いて プロットしてみると図のような軌跡が 描かれる 解析解 と同様の結果が得られていることがわかる 初期条件を変えて軌跡がどのように 変化するか試してみよ

42 lin([[v11,v12],[v21,v22],[v31,v32]],[k1,k2,k3]) 
Ruby ベクトルの線形結合を計算する。 lin([[v11,v12],[v21,v22],[v31,v32]],[k1,k2,k3])  [v11*k1+v21*k2+v31*k3, v12*k1+v22*k2+v32*k3] def lin(vs,xs) r = Array.new(vs[0].size) for i in 0..vs[0].size-1 sum = 0 for j in 0..vs.size-1 sum += vs[j][i]*xs[j] end r[i] = sum r def vx(v,x) lin([v],[x])

43 [u[2], u[3], -u[0]/r, -u[1]/r] } def runge_kutta(f,h,t,u,n)
Ruby func = lambda{|t,u| r = (u[0]**2+u[1]**2)**1.5 [u[2], u[3], -u[0]/r, -u[1]/r] } def runge_kutta(f,h,t,u,n) for i in 0..n-1 k1 = vx(f.call(t,u),h) k2 = vx(f.call(t+h/2,lin([u,k1],[1,1.0/2])),h) k3 = vx(f.call(t+h/2,lin([u,k2],[1,1.0/2])),h) k4 = vx(f.call(t+h,lin([u,k3],[1,1])),h) u = lin([u,k1,k2,k3,k4],[1,1.0/6,1.0/3,1.0/3,1.0/6]) t = t+h print u[0], "\t", u[1], "\n" end u runge_kutta(func,0.1,0,[1,0,0,1],100) f: 関数 h: 時間変数の差分 t: 時間変数の初期値 u: 初期値 n: ステップ数

44 数値解析(1) まとめ 数値計算には種々の誤差が生じる 数値積分の数値解法 常微分方程式の数値解法 実数表現 丸め誤差 打切り誤差 桁落ち
情報落ち 数値積分の数値解法 台形公式 シンプソン公式 常微分方程式の数値解法 Euler法 Runge-Kutta (Gill)法

45 backup

46 打切り誤差の例 Ruby 指数計算（自然対数の底）