静岡理工科大学 情報学部 菅沼ホーム 目次 索引

数値計算

  1. 1.連立 1 次方程式と逆行列
    1. 1.1 ガウス−ジョルダンの消去法
  2. 2.非線形方程式
    1. 2.1 二分法
    2. 2.2 Newton 法
    3. 2.3 セカント法
  3. 3.代数方程式と行列の固有値・固有ベクトル
    1. 3.1 ベアストウ法
    2. 3.2 行列の固有値(フレーム法+ベアストウ法)
    3. 3.3 実対称行列の固有値・固有ベクトル(ヤコビ法)
    4. 3.4 最大固有値と固有ベクトル(べき乗法)
  4. 4.数値微分
  5. 5.数値積分
    1. 5.1 台形則
    2. 5.2 シンプソン則
  6. 6.常微分方程式
    1. 6.1 ルンゲ・クッタ法
  7. 7.補間法
    1. 7.1 ラグランジュ補間法
    2. 7.2 スプライン補間法
    3. 7.3 ベジエ曲線
1.連立 1 次方程式と逆行列

1.1 ガウス−ジョルダンの消去法

    C/C++によるプログラム例 → 
    アプレット版及び JavaScript 版では,データを変更して,連立方程式の解,逆行列,行列の乗算,及び,行列式の値を計算し,その結果を画面上で得ることも可能です.

  以下に示すような連立一次方程式の解を求める方法について検討します.

a11x1 + a12x2 + ・・・ a1nxn = b1
a21x1 + a22x2 + ・・・ a2nxn = b2
       ・・・・・
an1x1 + an2x2 + ・・・ annxn = bn

上式は,行列 Ab,及び,x

    

のように定義すれば,以下のように書くことができます.

Ax = b

この方程式の解は,形式的に,A の逆行列を左から掛け,

x = A-1b

のようにして求めることができます.そこで,拡張された行列

を考え,左側の n×n の部分が単位行列になるような操作を行えば,最後の列に解が求まることになります.

  ガウス−ジョルダンの消去法は,この操作をコンピュータ上で行うものです.まず,拡張された行列の 1 行目の各要素を a11 で割ります.次に,a21 の値が 0 になるように,1 行目を何倍かして,2 行目に加えます( 2 行目から引きます).3 行目以降に対しても,同様の操作を繰り返します.この結果,下に示す図の最初の状態になるはずです.

次に,2 行目を a(1)22 で割った後同様の操作を行えば,上図の中央のような状態になります.以下,同様の操作を繰り返すことによって,上図の右端の状態になり,解が求まることになります.

  以上の操作によって,場合によっては,問題が発生します.つまり,a11,a22 等が 0 になる場合です.この場合は,行を入れ替え,aij の値が最大になるものが aii になるようにします.この aii のことを,ピボット要素といいます.ピボット要素が 0 であるとき(ある値より小さいとき)は,逆行列が存在しない(行列 A のランクが n より小さい)ことになります.

  また,拡張された行列において,b の部分を単位行列で置き換えれば,つまり,次のようにすれば,後半の n×n の部分に行列 A の逆行列が求まることになります.

2.非線形方程式

  この節では,非線形方程式 f(x) = 0 の根を求める方法について説明します.

2.1 二分法

    C/C++によるプログラム例 → 
    JavaScript 版では,JavaScript の仕様に適合した形で解を求めたい式を入力することによって,任意の非線形方程式の解を画面上で求めることができます.

  関数 f(x) が区間 [a, b] で連続で,かつ,根が 1 つだけ存在するような場合に利用できる方法です.このとき,以下のような関係が成立する必要があります.

f(a) f(b) < 0

アルゴリズムは以下に示す通りであり,徐々に解の存在区間を狭めることによって解を求める法法です(右図を参照).一般に,収束は遅いですが,上記の条件を満たせば必ず解に到達します.

  1. 区間 [a, b] の中点: c = 0.5 * (a + b) に対し,f(c) の値を計算

  2. f(a) f(c) < 0 なら
    b = c,f(b) = f(c)
    そうでなければ
    a = c,f(a) = f(c)

  3. |b - a| < ε ならば終了する.そうでなければ,1 に戻る.

2.2 Newton 法

    C/C++によるプログラム例 → 
    JavaScript 版では,JavaScript の仕様に適合した形で解を求めたい式及びその微分を入力することによって,任意の非線形方程式の解を画面上で求めることができます.

  Newton 法によるアルゴリズムは以下の通りです.右図に示すように,まず,初期値 x0 を与えます.次に,f(x0) を計算し,(x0, f(x0)) における接線が x 軸と交わる点を x1 とします.再び,x1 から同様の手続きを繰り返し,収束するまで続けます.一般的に記述すれば,以下のようになります.適切な初期値を設定すれば,非常に速く解を求めることが可能です(勿論,関数の形にもよりますが).

xn+1 = xn - f(xn) / f'(xn)

収束判定方法としては,以下のような方法があります.

|xn+1 - xn| < ε
|f(xn+1) - f(xn)| < ε
|f(xn)| < ε

    多次元(x がベクトル)の場合に対する C/C++ によるプログラム例 → 
    JavaScript 版では,JavaScript の仕様に適合した形で解を求めたい式及びその微分を入力することによって,任意の非線形方程式の解を画面上で求めることができます.例えば,次数が 3 である場合は,解を求めたい式 f1(x1, x2, x3),f2(x1, x2, x3),f3(x1, x2, x3) と,その微分,
       
    を 1 行 1 列目,1 行 2 列目,・・・,3 行 3 列目の順で入力して下さい.また,各変数は配列の形で記述して下さい.

2.3 セカント法

    C/C++によるプログラム例 → 
    JavaScript 版では,JavaScript の仕様に適合した形で解を求めたい式を入力することによって,任意の非線形方程式の解を画面上で求めることができます.

  二分法と同様,関数 f(x) が区間 [a, b] で連続で,かつ,根が 1 つだけ存在するような場合に利用できる方法です.ただし,二分法とは異なり,中点の関数値ではなく,Newton 法における計算方法の微分を差分で置き換えた方法を利用します.アルゴリズムは以下に示す通りです.

  1. c = b - f(b) * (b - a) / (f(b) - f(a))

  2. f(a) f(c) < 0 なら
    b = c,f(b) = f(c)
    そうでなければ
    a = c,f(a) = f(c)

  3. |b - a| < ε ならば終了する.そうでなければ,1 に戻る.

3.代数方程式と行列の固有値・固有ベクトル

3.1 ベアストウ法

    C/C++によるプログラム例 → 
    アプレット版及び JavaScript 版では,任意の実係数 n 次代数方程式の根を画面上で求めることができます.

  実係数 n 次代数方程式,

f(x) = a0xn + a1xn-1 + ・・・ + an-1x + an = 0

のすべての根を求める方法としてベアストウ法があります.ベアストウ法は,対象とする式を,2 次式,

x2 + px + q

と (n - 2) 次式の積の形に変換することを繰り返して解く方法です.p と q の値を決めるために,その初期値 p0 と q0 を与え,f(x) を x2 + p0x + q0 で割った余りが 0 になるように,p0 と q0 を徐々に修正していきます.具体的なアルゴリズムは,以下に示す通りです.

  1. p0, q0 を与え,k = 1 とおく

  2. bi = ai - pk-1bi-1 - qk-1bi-2,   i = 0, 1, ・・・, n,  b-1 = b-2 = 0

  3. ci = bi - pk-1ci-1 - qk-1ci-2,   i = 0, 1, ・・・, n,  c-1 = c-2 = 0

  4. D = cn-22 - cn-3(cn-1 - bn-1)
    Δp = (bn-1cn-2 - bncn-3) / D
    Δq = (bncn-2 - bn-1(cn-1 - bn-1)) / D

  5. もし,|Δp| < ε, かつ, |Δq| < ε であれば,

    1. 2 次方程式 x2 + pk-1x + qk-1 = 0 を解く

    2. b0, ・・・, bn-2 を,新たに,a0, ・・・, an-2 とおき,かつ,n = n - 2 とし,1 へ戻る.ただし,n が 1 または 2 の場合は,その方程式を解く.

    そうでなければ,

    1. pk = pk-1 + Δp
      qk = qk-1 + Δq

    2. k = k + 1 とし,2 へ戻る.

3.2 行列の固有値(フレーム法+ベアストウ法)

    C/C++によるプログラム例 → 
    アプレット版及び JavaScript 版では,任意の行列の固有値を画面上で求めることができます.

  行列 A の固有ベクトルを x,固有値を λ とすると,

Ax = λx

という関係が成立します.この方程式が零ベクトル以外の解をもつための必要十分条件は,

|A - λI| = 0

となることであり,この方程式を A特性方程式固有方程式)といいます.そして,この方程式の n 個の根 λ1,λ2,・・・,λn は,A の固有値になっています.

  従って,固有値を求めることは,この特性方程式を解くことになります.特性方程式は,

|A - λI| = λn + b1λn-1 + ・・・ + bn-1λ + bn = 0

と書き表せます.b1 から bn の値が決まれば,ベアストウ法によって上式の根を求めることができます.b1 から bn をフレーム法で求め,ベアストウ法で特性方程式の根を求めるといった方法で固有値を求めるのが以下に示すアルゴリズムです.

  1. b1 = -trA とする(「tr」はトレースを表し,対角要素の和である)

  2. Hk-1 = AHk-2 + bk-1I
    bk = -tr(AHk-1) / k   k = 2, ・・・ n, ただし,H0 = I

  3. ベアストウ法で固有値を求める.

3.3 実対称行列の固有値・固有ベクトル(ヤコビ法)

    C/C++によるプログラム例 → 
    アプレット版及び JavaScript 版では,任意の実対称行列の固有値・固有ベクトルを画面上で求めることができます.各固有値に対応する固有ベクトルは,各列に表示されます.

  n 次の実対称行列(要素がすべて実数である対称行列であり,固有値はすべて実数となる) A の固有ベクトルを xi,固有値を λi とすると,

Axi = λixi   i = 1, 2, ・・・, n

という関係が成立します.今,

R = [x1, x2, ・・・, xn]

とし,各 xi が正規直行系をなしていれば,

とすると,

R-1AR = RTAR = Λ

という関係が成立します.従って,直交行列 R を求めれば,固有値と固有ベクトルが求まることになります.

  直交行列 R を求めるために,まず,行列 A の非対角要素の内,絶対値が最大になる要素を 0 にするような直交行列 Rk を順次作成していき,最終的に R を求めます.つまり,

A1 = R1TAR1
A2 = R2TA1R2 = R2TR1TAR1R2
     ・・・
Ak = (R1R2・・・Rk)TA(R1R2・・・Rk)

と計算していき,もし,Ak が対角行列になっていれば,求める R は,R1R2・・・Rk となります.

  具体的アルゴリズムは,以下のようになります.

  1. k = 1, A0 = AX0 = I とおく.

  2. 行列 Ak-1 の非対角要素の内,絶対値が最大になる要素 apq を見つける.

  3. もし,|| < ε ならば,
    1. 終了する(Ak-1 の対角要素が固有値で,かつ,Xk-1 の各列が対応する固有ベクトル)
    そうでなければ,
    1. 要素 apq を 0 にするような直交変換
      Ak = RkTAk-1Rk
      を行う.このとき,Ak の各要素の値は,以下のようになる.
      s = -apq
      t = 0.5(app - aqq)
      v = |t| / (s2 + t2)0.5
      sinθ = (0.5(1 - v))0.5・sign(st)
      cosθ = (1 - sin2θ)0.5
      apj(k) = apj(k-1)cosθ - aqj(k-1)sinθ
      aqj(k) = aqj(k-1)cosθ + apj(k-1)sinθ
      aip(k) = aip(k-1)cosθ - aiq(k-1)sinθ
      aiq(k) = aiq(k-1)cosθ + aip(k-1)sinθ
      aij(k) = aij(k-1)   以上において, i,j ≠ p,q
      app(k) = app(k-1)cos2θ + aqq(k-1)sin2θ - 2apq(k-1)sinθcosθ
      aqq(k) = app(k-1)sin2θ + aqq(k-1)cos2θ + 2apq(k-1)sinθcosθ
      apq(k) = 0
      aqp(k) = 0
    2. Xk を,
      Xk = Xk-1Rk
      として計算する.このとき,Xk の各要素の値は,以下のようになる.
      xjp(k) = xjp(k-1)cosθ - xjp(k-1)sinθ
      xjq(k) = xjq(k-1)cosθ + xjq(k-1)sinθ   j = 1, 2, ・・・, n
    3. 2 へ戻る.

3.4 最大固有値と固有ベクトル(べき乗法)

    C/C++によるプログラム例 → 
    アプレット版及び JavaScript 版では,任意の行列の最大固有値と固有ベクトルを画面上で求めることができます.適当にパラメータを変更することによって,すべての固有値及び固有ベクトルを求めることも可能です.

  べき乗法は,任意の初期ベクトル v0 から始め,

vk = Avk-1 / ‖Avk-1‖   k = 0, 1, 2, ・・・

という手続きを繰り返すと,‖Avk-1‖,及び,vk が,行列 A の最大固有値の絶対値 |λ|,及び,その対応する固有ベクトル v に収束することを利用した方法です(初期ベクトルの与え方によっては,必ずしも最大固有値に収束しない).最大固有値から,順に,複数の固有値,固有ベクトルを求めることが可能です.具体的アルゴリズムは,以下のようになります.

  なお,この方法は,固有値が 0 の場合や,重複固有値がある場合に対しては適用できません.また,A-1 の固有値は,A の固有値の逆数になりますから,最小固有値を求める方法としても利用できます(もちろん,逆行列が存在する場合に限ります).

  1. v0 を与え,B = I とする.

  2. l0 = ‖v0‖,k = 0 とする.

  3. vk = Bvk / ‖Bvk‖  (丸め誤差を取り除くためであり,必ずしも毎回行わなくても良い)

  4. lk+1 = ‖Avk‖,  vk+1 = Avk / lk+1

  5. もし,|(lk+1 - lk) / lk| < ε ならば,

    1. vk+1 の要素の内,0 でないものを見つける.

    2. 対応する vk の要素の符号が,vk+1 のものと異なっていれば,lk+1 = -lk+1 とする.

    3. 固有値及び固有ベクトルとして,lk+1 及び vk+1 を記憶する.

    4. B = (A - lk+1I)B,  v0 = Bv0 として,2 へ戻る(次の固有値,固有ベクトルを求める手続き).

    そうでなければ,

    1. k = k + 1 として,3 へ戻る.

4.数値微分

  数値微分は,微分の定義に従い,以下のような形で計算します.

5.数値積分

5.1 台形則

    C/C++によるプログラム例 → 
    JavaScript 版では,JavaScript の仕様に適合した形で積分したい式を入力することによって,任意の関数の積分を画面上で求めることができます.

  台形則では,積分区間 [a, b] を n 等分し,右図に示すように,台形の面積の和で近似して計算します.つまり,以下の式を計算することになります.

     h = (b - a) / n,  xi = a + ih

5.2 シンプソン則

    C/C++によるプログラム例 → 
    JavaScript 版では,JavaScript の仕様に適合した形で積分したい式を入力することによって,任意の関数の積分を画面上で求めることができます.

  シンプソン則による積分方法は以下のようになります.積分区間 [a, b] を 2n 等分し,

h = (b - a) / 2n,  xi = a + ih

としたとき,たとえば,3 点,

(x0, f(x0)), (x1, f(x1)), (x2, f(x2))

用いて,f(x) を 3 点を通る 2 次式 g(x) で近似すると,ラグランジュの補間公式によって,以下のようになります.

  この g(x) の積分によって,区間 [x0, x2] における f(x) の積分を近似すると以下のようになります.

以上より,区間 [a, b] の積分は,以下の式によって近似することができます.

6.常微分方程式

6.1 ルンゲ・クッタ法

    C/C++によるプログラム例 → 
    JavaScript 版では,JavaScript の仕様に適合した形で微分方程式を入力することによって,任意の微分方程式の解を画面上求めることができます.具体的には,
       
    の右辺を,変数には配列を使用して入力することになります.

  ルンゲ・クッタ法について説明する前に,オイラー法によって微分方程式を解く基本的な考え方について説明します.微分方程式,

が与えられたとします.微分の定義によって,以下のような近似式が得られます.

  オイラー法では,この式に従って微分方程式を解いていきます.まず,初期値 t0,x(t0) (= x0) を与えます.すると,与えられた微分方程式から,t0 における接線の傾き f(t0) を計算できます(右図参照).これらのデータから,t0+h における x(t0+h) の推定値を計算します.

  以下,t0+h を 新たに t0 と考え,この手続きを繰り返していきます.図からも明らかなように,h が大きいと,真の x(t0+h) とその推定値の間には大きな差がありますが,h を十分小さくとることによって,この差を小さくすることができます.

  オイラー法では,t0 における傾き f(t0) だけを利用して,x(t0+h) を推定していますが,ルンゲ・クッタ法では,以下に示すように,複数点における傾きを利用しています.

g1 = hf(t0, x0)
g2 = hf(t0 + h / 2, x0 + g1 / 2)
g3 = hf(t0 + h / 2, x0 + g2 / 2)
g4 = hf(t0 + h, x0 + g3)
x(t0+h) = x0 + (g1 + 2g2 + 2g3 + g4) / 6

7.補間法

7.1 ラグランジュ補間法

    C/C++によるプログラム例 → 
    アプレット版及び JavaScript 版では,任意のデータに対して,n 次補間多項式による計算結果を画面上で求めることができます.

  ラグランジュ補間法は,(n+1) 個の点 (xi, yi) (i = 0, 1, ・・・, n) が与えられたとき,与えられた点を通る n 次多項式,

f(x) = a0xn + a1xn-1 + ・・・ + an-1x + an

によって,補間する方法です.x が与えられたとき,f(x) は以下の式によって計算されます.

7.2 スプライン補間法

    C/C++によるプログラム例 → 
    アプレット版及び JavaScript 版では,任意のデータに対して,スプライン補間法による計算結果を画面上で求めることができます.

  スプライン補間法は,与えられた点を通るなるべくなめらかな曲線によって,補間しようとするものです.具体的には,以下のような多項式によって補間します.

  具体的なアルゴリズムは以下のようになります.

  1. ステップ 1

    hi = xi+1 - xi   i = 0,・・・,n-1
    bi = 2(hi + hi-1)   i = 1,・・・,n-1
    di = 3{(yi+1 - yi) / hi - (yi - yi-1) / hi-1}   i = 1,・・・,n-1

  2. ステップ 2

    g1 = h1 / b1
    gi = hi / (bi - hi-1gi-1)   i = 2,・・・,n-2
    u1 = d1 / b1
    ui = (di - hi-1ui-1) / (bi - hi-1gi-1)   i = 2,・・・,n-1

  3. ステップ 3

    rn-1 = un-1
    ri = ui - giri+1   i = n-2,・・・,1
    r0 = 0

  4. ステップ 4  xi < x < xi+1 のとき

    pi = yi = f(xi)
    qi = (yi+1 - yi) / hi - hi(ri+1 + 2ri) / 3
    si = (ri+1 - ri) / 3hi
    y = pi + qi(x - xi) + ri(x - xi)2 + si(x - xi)3

7.3 ベジエ曲線

    C/C++によるプログラム例 → 
    アプレット版及び JavaScript 版では,任意のデータに対して,ベジエ曲線上の座標を画面上に出力することができます.

  ここで述べるベジエ曲線は,先の節で述べたような補間法とは異なります.つまり,与えられた点を通る曲線を描くのではなく,与えられた点(B0,B1,・・・,Bn)で構成される定義多角形(ベジエ多角形)から,次式によって決まる曲線を描きます.

P(t) = ΣBiJn , i(t)   0 ≦ t ≦ 1

ここで,P(t) は曲線上の点の座標であり,t = 0 の時は B0 に,また,t = 1 の時は Bn に一致します.また,Jn , i(t) は,ベジエ基底関数と呼ばれ,以下の式で計算されます.

  たとえば,n = 3 とし,ベジエ多角形を B0 = (1 1),B1 = (2 3),B2 = (4 3),B3 = (3 1) とすると,ベジエ曲線は以下のようになります(プログラム例は,この曲線を描くためのものです).

  なお,ベジエ曲線には,以下のような性質があります.

  ベジエ曲線を複数個つなげて,下の図に示すような,複雑な曲線を描くこともできます.その際,接続点で,1 階の導関数が一致するためには,図に示すように,最初の定義多角形の最後の辺の方向と大きさが,2 番目の定義多角形の最初の辺の方向と大きさに一致していなければなりません(B3 - B2 = C1 - C0).

静岡理工科大学 情報学部 菅沼ホーム 目次 索引