授業資料/03 のバックアップの現在との差分(No.1)

バックアップ一覧
差分を表示
ソースを表示
バックアップを表示
授業資料/03 へ行く。
- 1 (2015-01-19 (月) 07:49:57)
追加された行はこの色です。
削除された行はこの色です。
aa
* 連立一次方程式の求解 [#tdae5197]

** LU 分解 [#f759ffe4]

掃き出し法を線形代数できちんと捉え直し、かつ、なるべく便利な形に直すと、LU 分解というものになる．

*** サンプルプログラム [#f0e1cf05]
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# LU decomposition test program.
# Solve a problem ax=b where a is a matrix and b is a vector.

include Math

# size of the problem
n = 4

# matrix
a = [
  [3.0, 3.0, -5.0, -6.0],
  [1.0, 2.0, -3.0, -1.0],
  [2.0, 3.0, -5.0, -3.0],
  [-1.0, 0, 0, 1.0 ]
]

# vector
b = [1.0, 2.0, 3.0, 4.0]

# LU decomposition without pivot operation. 
# Note that we don't need the vector b and computation cost order is n^3 because it is triple loop.
for k in (1..(n-1))
  for i in (k+1..n)
    a[i-1][k-1] /= a[k-1][k-1]
    for j in (k+1..n)
      a[i-1][j-1] -= a[i-1][k-1]*a[k-1][j-1]
    end
  end
end

# We obtain y=L^(-1)b by a forward substitution.
# Its loop is double and it means its computation cost order is n^2.
y = [0.0, 0.0, 0.0, 0.0]
for i in (1..n)
  y[i-1] = b[i-1]
  if i < n
    for j in ((i+1)..n)
      b[j-1] -= a[j-1][i-1]*b[i-1]
    end
  end
end

# Solve the solution x=U^(-1)y by a backward substitution.
# Its computation cost order is also n^2.
x = [0.0, 0.0, 0.0, 0.0]
n.downto(1) { |i|
  x[i-1] = y[i-1]/a[i-1][i-1]
  if i > 2
    (i-1).downto(1){ |j|
    y[j-1] -= (a[j-1][i-1]/a[i-1][i-1])*y[i-1]
    }
  end
}

# result
p x

# Sample result
# [-4.999999999999999, 8.881784197001252e-16, -1.9999999999999993, -1.0]
}}

** 最急降下法 [#l1699d2e]

行列 A が正定値対称の場合、Ax = b の真の解 x に対して S(a) := (a-x)^t A (a-x) とすると、
この S(a) は非負であり、a = x の時のみゼロになる．
これを利用して、S(a) が小さくなるようにベクトル a を反復改良していって、S(a) がほぼゼロになればその a は x の近似値とできる．

この考え方に沿って、S(a) を「最も早く小さくする方向」に a を改善していく手法である．&br;
&ref(/materials/warning.png); S(a) をもっとも早く小さくする方向とは S(a) の勾配の逆方向で、計算してみると b - Aa である．

ちなみに、この方向を「もう少し賢く」選んだものが後述の CG法であるとみなせる．

*** サンプルプログラム [#v5866189]
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# Steepest descent method to solve Ax = b where A: matrix, b: vector.
include Math

##### contants #####
# Allowed error which is evaluated by norm.
Eps = 1.0e-08

# Matrix
A = [
  [7.0, 1.0, 2.0, 3.0],
  [1.0, 10.0, 3.0, 4.0],
  [2.0, 3.0, 13.0, 6.0],
  [3.0, 4.0, 6.0, 16.0 ]
]

# Vector on the right hand side
B = [1.0, 2.0, 3.0, 4.0]

# size of the problem
N = A.size

##### Functions ######
### Function: inner product
def inner_product(x,y)
  z = 0.0
  for i in (1..N)
    z += x[i-1]*y[i-1]
  end
  return z
end

### Function: scalar multiply
def vc_coeff(c,b)
  z = Array.new(N).collect{ 0.0 }
  for i in (1..N)
    z[i-1] = c * b[i-1]
  end
  return z
end

### Function: sum of vectors
def vc_sum(a,b)
  c = Array.new(N).collect{ 0.0 }
  for i in (1..N)
    c[i-1] = a[i-1]+b[i-1]
  end
  return c
end

### Function: matrix * vector
def mtrx_vc_product(a,b)
  c = Array.new(N).collect{ 0.0 }
  for i in (1..N)
    for j in (1..N)
      c[i-1] += a[i-1][j-1]*b[j-1]
    end
  end

  return c
end

### Function: A * vector
def a_vc(b)
  return mtrx_vc_product(A,b)
end

### Function: residual = b - Ax
def residual(x)
  return vc_sum(B, vc_coeff(-1.0, a_vc(x)))
end

### Function: The norm to use in the Jacobi method.
def norm(x)
  v = 0.0
  for i in (1..N)
    v += x[i-1]**2
  end
  return sqrt(v)
end

### Function: print results
def print_result(x)
  print "x = ["
  for i in (1..N-1)
    printf("%.6f",x[i-1])
    print ", "
  end
  printf("%.6f",x[N-1])
  print "], "
  print "residual= ",norm(residual(x)),"\n"
end

### Function: Steepest descent method iteration
def steepest_descent_iteration(x)
  r = residual(x)
  w = inner_product(r,r) / inner_product(r, a_vc(r))

  return vc_sum(x, vc_coeff(w,r))
end

##### main program #####
# initial values
x = Array.new(N).collect{ 1.0 }
print_result(x)

# loop to obtain a good x.
loop_times = 0
loop do
  if (norm(residual(x)) < Eps) then
      break
    else
      x = steepest_descent_iteration(x)
      print_result(x)
      loop_times += 1
    end
end

# result
print "\n\nYou obtain a numerical solution after ",loop_times," loops.\n\n"
print_result(x)
print "\n"
}}

** CG 法 [#fc9cfaef]

S(a) を最も早く小さくする方向に改善する最急降下法は欲張りすぎでかえって収束が遅いのが常と言って良いので、数学的にきちんと考えなおしたのが CG 法である．
具体的には、一回前に a を改善してきた「方向」に対し、次の改善方向を A-直交するようにしたものが CG 法である．&br;
&ref(/materials/warning.png); 最急降下法は、先の改善方向に対し次の改善方向が必ず直交している．&br;
&ref(/materials/warning.png); 二つのベクトル x,y が A-直交であるとは、x^t A y = 0 であること．

さて、このように改善方向を A-直交になるようにすると、「原理的には」A の次元数を n とするとベクトル a を n 回改善するかそれまでにかならず a = x となるという大変良い性質がもたらされる．

*** サンプルプログラム [#oe1111a1]
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# CG method to solve Ax = b where A: matrix, b: vector.
include Math

##### contants #####
# Allowed error which is evaluated by norm.
Eps = 1.0e-08

# Matrix
A = [
  [7.0, 1.0, 2.0, 3.0],
  [1.0, 10.0, 3.0, 4.0],
  [2.0, 3.0, 13.0, 6.0],
  [3.0, 4.0, 6.0, 16.0 ]
]

# Vector on the right hand side
B = [1.0, 2.0, 3.0, 4.0]

# size of the problem
N = A.size

##### Functions ######
### Function: inner product
def inner_product(x,y)
  z = 0.0
  for i in (1..N)
    z += x[i-1]*y[i-1]
  end
  return z
end

### Function: scalar multiply
def vc_coeff(c,b)
  z = Array.new(N).collect{ 0.0 }
  for i in (1..N)
    z[i-1] = c * b[i-1]
  end
  return z
end

### Function: sum of vectors
def vc_sum(a,b)
  c = Array.new(N).collect{ 0.0 }
  for i in (1..N)
    c[i-1] = a[i-1]+b[i-1]
  end
  return c
end

### Function: matrix * vector
def mtrx_vc_product(a,b)
  c = Array.new(N).collect{ 0.0 }
  for i in (1..N)
    for j in (1..N)
      c[i-1] += a[i-1][j-1]*b[j-1]
    end
  end

  return c
end

### Function: A * vector
def a_vc(b)
  return mtrx_vc_product(A,b)
end

### Function: residual = b - Ax
def residual(x)
  return vc_sum(B, vc_coeff(-1.0, a_vc(x)))
end

### Function: The norm to use in the Jacobi method.
def norm(x)
  v = 0.0
  for i in (1..N)
    v += x[i-1]**2
  end
  return sqrt(v)
end

### Function: print results
def print_result(x)
  print "x = ["
  for i in (1..N-1)
    printf("%.6f",x[i-1])
    print ", "
  end
  printf("%.6f",x[N-1])
  print "], "
  print "residual= ",norm(residual(x)),"\n"
end

### Function: CG iteration
def cg_iteration(x,r,p,r_R,q)
  q_Q = inner_product(p,q)
  w = r_R/q_Q
  x_new = vc_sum(x, vc_coeff(w,p))
  r_new = vc_sum(r, vc_coeff((-1.0*w),q))

  r_R_new = inner_product(r_new, r_new)
  beta = r_R_new / r_R
  p_new = vc_sum(r_new, vc_coeff(beta,p))
  q_new = a_vc(p_new)
  
  return [x_new, r_new, p_new, r_R_new, q_new]
end

##### main program #####
# initial values
x = Array.new(N).collect{ 1.0 }
r = residual(x)
p = vc_coeff(1.0, r)
r_R = inner_product(r,r)
q = a_vc(p)

print_result(x)

# loop to obtain a good x.
loop_times = 0
loop do
  if (norm(residual(x)) < Eps) then
      break
    else
      x,r,p,r_R,q = cg_iteration(x,r,p,r_R,q)
      print_result(x)
      loop_times += 1
    end
end

# result
print "\n\nYou obtain a numerical solution after ",loop_times," loops.\n\n"
print_result(x)
print "\n"
}}

// ━┃┏┓┛┗┣┳┫┻╋


// コマンドライン入力は「行頭を > で始める」.
// コマンドライン出力は「行頭をブランクで始める」.

// 実習アイコン
// &ref(/materials/notes.png); 

// 注意アイコン
// &ref(/materials/warning.png);

// Link アイコン
// &ref(/materials/JNorth_arrow-right-sm.png); 

// OK アイコン
// &ref(/materials/OK.png);

// NG アイコン
// &ref(/materials/NG.png);

// サンプルアイコン
// &ref(/materials/Gnome-Preferences.png);

// 大文字での強調 
// CENTER:&size(24){''ほげほげ''};

// programu source 表記
// #highlighter(language=ruby,number=on,cache=on){{}}

// パイプ
// &brvbar;#contents