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

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

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

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

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

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

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

丸め誤差の例 誤差 # 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", "01111111011", "10011001100…110011010"] 0捨1入 ["0", "01111111101", "001100110011…00110100"] 誤差 ["0", "01111111101", "0011001100…1100110011"] 丸め誤差は # （検証）irbで以下を試してみよ 0.1*3-0.3 2**-54 0が52個

..0110011010 ..1100110100 ------------------- ..0011001110 ..00110100

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

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

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

桁落ち誤差の例 桁落ち誤差 打切り誤差 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" } 　　　を小さくしていくと、 ある程度までは打切り誤差が 減少していくが、より微小になると 　　　　　と　　　　　の差も微小となり 有効桁数が減少する その結果、桁落ち誤差が大きくなる ※打切り誤差については 後ほど述べる

情報落ち誤差の例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 = 1.0000001 n = 100 d = 0.01 sum = 0.99999934 n = 1000 d = 0.0010 sum = 0.9999907 n = 10000 d = 1.0E-4 sum = 1.0000535 n = 100000 d = 1.0E-5 sum = 1.0009902 n = 1000000 d = 1.0E-6 sum = 1.0090389 n = 10000000 d = 1.0E-7 sum = 1.0647675 n = 100000000 d = 1.0E-8 sum = 0.25 n = 1000000000 d = 1.0E-9 sum = 0.03125 実行結果 丸め誤差 情報落ち誤差

情報落ち誤差の例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,…,1000000 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=1 d=1.0 sum=1.0 n=10 d=0.1 sum=1.0 n=100 d=0.01 sum=1.0 n=1000 d=0.001 sum=1.0 n=10000 d=0.0001 sum=0.999999999999906 n=100000 d=1.0e-05 sum=0.999999999998084 n=1000000 d=1.0e-06 sum=1.00000000000792 丸め誤差

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

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

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

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

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

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

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

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

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

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

簡単な例 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

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

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

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

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

簡単な例 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

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

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

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

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

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

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

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

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 } ※簡単のためここでは　　　　としている

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

Runge-Kutta法の例（続き） 実行結果 であったから は質点の軌跡を表す をgnuplotやexcelなどを用いて Ruby 実行結果 t= 0.1 u[0]= 0.995004158961472 u[1]= 0.0998333321663404 u[2]= -0.0998334893797951 u[3]= 0.995004158786358 t= 0.2 u[0]= 0.980066566148056 u[1]= 0.19866916135847 u[2]= -0.19866948092869 u[3]= 0.980066565047173 t= 0.3 u[0]= 0.955336472286007 u[1]= 0.295519952547433 u[2]= -0.295520442490704 u[3]= 0.955336467918296 t= 0.4 u[0]= 0.921060971123942 u[1]= 0.389418004245675 u[2]= -0.389418675043335 u[3]= 0.921060958695548 t= 0.5 u[0]= 0.877582530631772 u[1]= 0.479425117337853 u[2]= -0.479425981331435 u[3]= 0.877582502133123 であったから は質点の軌跡を表す 　　　　　をgnuplotやexcelなどを用いて プロットしてみると図のような軌跡が 描かれる 解析解 と同様の結果が得られていることがわかる 初期条件を変えて軌跡がどのように 変化するか試してみよ

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])

[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: ステップ数

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

backup

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