Ȼ/03

Top / Ȼ / 03

ϢΩ켡ε

LU ʬ

ݤФˡǤªľġʤ٤ʷľȡLU ʬȤΤˤʤ롥

ץץ

  1. # LU decomposition test program.
  2. # Solve a problem ax=b where a is a matrix and b is a vector.
  3.  
  4. include Math
  5.  
  6. # size of the problem
  7. n = 4
  8.  
  9. # matrix
  10. a = [
  11.   [3.0, 3.0, -5.0, -6.0],
  12.   [1.0, 2.0, -3.0, -1.0],
  13.   [2.0, 3.0, -5.0, -3.0],
  14.   [-1.0, 0, 0, 1.0 ]
  15. ]
  16.  
  17. # vector
  18. b = [1.0, 2.0, 3.0, 4.0]
  19.  
  20. # LU decomposition without pivot operation.
  21. # Note that we don't need the vector b and computation cost order is n^3 because it is triple loop.
  22. for k in (1..(n-1))
  23.   for i in (k+1..n)
  24.     a[i-1][k-1] /= a[k-1][k-1]
  25.     for j in (k+1..n)
  26.       a[i-1][j-1] -= a[i-1][k-1]*a[k-1][j-1]
  27.     end
  28.   end
  29. end
  30.  
  31. # We obtain y=L^(-1)b by a forward substitution.
  32. # Its loop is double and it means its computation cost order is n^2.
  33. y = [0.0, 0.0, 0.0, 0.0]
  34. for i in (1..n)
  35.   y[i-1] = b[i-1]
  36.   if i < n
  37.     for j in ((i+1)..n)
  38.       b[j-1] -= a[j-1][i-1]*b[i-1]
  39.     end
  40.   end
  41. end
  42.  
  43. # Solve the solution x=U^(-1)y by a backward substitution.
  44. # Its computation cost order is also n^2.
  45. x = [0.0, 0.0, 0.0, 0.0]
  46. n.downto(1) { |i|
  47.   x[i-1] = y[i-1]/a[i-1][i-1]
  48.   if i > 2
  49.     (i-1).downto(1){ |j|
  50.     y[j-1] -= (a[j-1][i-1]/a[i-1][i-1])*y[i-1]
  51.     }
  52.   end
  53. }
  54.  
  55. # result
  56. p x
  57.  
  58. # Sample result
  59. # [-4.999999999999999, 8.881784197001252e-16, -1.9999999999999993, -1.0]

ǵ޹߲ˡ

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 ƤˡǤ롥
warning.png S(a) äȤ᤯Ȥ S(a) θۤεǡ׻Ƥߤ b - Aa Ǥ롥

ʤߤˡ֤⤦ΤҤ CGˡǤȤߤʤ롥

ץץ

  1. # Steepest descent method to solve Ax = b where A: matrix, b: vector.
  2. include Math
  3.  
  4. ##### contants #####
  5. # Allowed error which is evaluated by norm.
  6. Eps = 1.0e-08
  7.  
  8. # Matrix
  9. A = [
  10.   [7.0, 1.0, 2.0, 3.0],
  11.   [1.0, 10.0, 3.0, 4.0],
  12.   [2.0, 3.0, 13.0, 6.0],
  13.   [3.0, 4.0, 6.0, 16.0 ]
  14. ]
  15.  
  16. # Vector on the right hand side
  17. B = [1.0, 2.0, 3.0, 4.0]
  18.  
  19. # size of the problem
  20. N = A.size
  21.  
  22. ##### Functions ######
  23. ### Function: inner product
  24. def inner_product(x,y)
  25.   z = 0.0
  26.   for i in (1..N)
  27.     z += x[i-1]*y[i-1]
  28.   end
  29.   return z
  30. end
  31.  
  32. ### Function: scalar multiply
  33. def vc_coeff(c,b)
  34.   z = Array.new(N).collect{ 0.0 }
  35.   for i in (1..N)
  36.     z[i-1] = c * b[i-1]
  37.   end
  38.   return z
  39. end
  40.  
  41. ### Function: sum of vectors
  42. def vc_sum(a,b)
  43.   c = Array.new(N).collect{ 0.0 }
  44.   for i in (1..N)
  45.     c[i-1] = a[i-1]+b[i-1]
  46.   end
  47.   return c
  48. end
  49.  
  50. ### Function: matrix * vector
  51. def mtrx_vc_product(a,b)
  52.   c = Array.new(N).collect{ 0.0 }
  53.   for i in (1..N)
  54.     for j in (1..N)
  55.       c[i-1] += a[i-1][j-1]*b[j-1]
  56.     end
  57.   end
  58.  
  59.   return c
  60. end
  61.  
  62. ### Function: A * vector
  63. def a_vc(b)
  64.   return mtrx_vc_product(A,b)
  65. end
  66.  
  67. ### Function: residual = b - Ax
  68. def residual(x)
  69.   return vc_sum(B, vc_coeff(-1.0, a_vc(x)))
  70. end
  71.  
  72. ### Function: The norm to use in the Jacobi method.
  73. def norm(x)
  74.   v = 0.0
  75.   for i in (1..N)
  76.     v += x[i-1]**2
  77.   end
  78.   return sqrt(v)
  79. end
  80.  
  81. ### Function: print results
  82. def print_result(x)
  83.   print "x = ["
  84.   for i in (1..N-1)
  85.     printf("%.6f",x[i-1])
  86.     print ", "
  87.   end
  88.   printf("%.6f",x[N-1])
  89.   print "], "
  90.   print "residual= ",norm(residual(x)),"\n"
  91. end
  92.  
  93. ### Function: Steepest descent method iteration
  94. def steepest_descent_iteration(x)
  95.   r = residual(x)
  96.   w = inner_product(r,r) / inner_product(r, a_vc(r))
  97.  
  98.   return vc_sum(x, vc_coeff(w,r))
  99. end
  100.  
  101. ##### main program #####
  102. # initial values
  103. x = Array.new(N).collect{ 1.0 }
  104. print_result(x)
  105.  
  106. # loop to obtain a good x.
  107. loop_times = 0
  108. loop do
  109.   if (norm(residual(x)) < Eps) then
  110.       break
  111.     else
  112.       x = steepest_descent_iteration(x)
  113.       print_result(x)
  114.       loop_times += 1
  115.     end
  116. end
  117.  
  118. # result
  119. print "\n\nYou obtain a numerical solution after ",loop_times," loops.\n\n"
  120. print_result(x)
  121. print "\n"

CG ˡ

S(a) Ǥ᤯˲ǵ޹߲ˡĥꤹǤäƼ«٤ΤȸäɤΤǡŪˤȹͤʤΤ CG ˡǤ롥 Ūˤϡ a ƤפФβ A-ľ򤹤褦ˤΤ CG ˡǤ롥
warning.png ǵ޹߲ˡϡβФβɬľ򤷤Ƥ롥
warning.png ĤΥ٥ȥ x,y A-ľǤȤϡx^t A y = 0 Ǥ뤳ȡ

ơΤ褦˲ A-ľˤʤ褦ˤȡָŪˤϡA μ n Ȥȥ٥ȥ a n 뤫ޤǤˤʤ餺 a = x ȤʤȤɤ⤿餵롥

ץץ

  1. # CG method to solve Ax = b where A: matrix, b: vector.
  2. include Math
  3.  
  4. ##### contants #####
  5. # Allowed error which is evaluated by norm.
  6. Eps = 1.0e-08
  7.  
  8. # Matrix
  9. A = [
  10.   [7.0, 1.0, 2.0, 3.0],
  11.   [1.0, 10.0, 3.0, 4.0],
  12.   [2.0, 3.0, 13.0, 6.0],
  13.   [3.0, 4.0, 6.0, 16.0 ]
  14. ]
  15.  
  16. # Vector on the right hand side
  17. B = [1.0, 2.0, 3.0, 4.0]
  18.  
  19. # size of the problem
  20. N = A.size
  21.  
  22. ##### Functions ######
  23. ### Function: inner product
  24. def inner_product(x,y)
  25.   z = 0.0
  26.   for i in (1..N)
  27.     z += x[i-1]*y[i-1]
  28.   end
  29.   return z
  30. end
  31.  
  32. ### Function: scalar multiply
  33. def vc_coeff(c,b)
  34.   z = Array.new(N).collect{ 0.0 }
  35.   for i in (1..N)
  36.     z[i-1] = c * b[i-1]
  37.   end
  38.   return z
  39. end
  40.  
  41. ### Function: sum of vectors
  42. def vc_sum(a,b)
  43.   c = Array.new(N).collect{ 0.0 }
  44.   for i in (1..N)
  45.     c[i-1] = a[i-1]+b[i-1]
  46.   end
  47.   return c
  48. end
  49.  
  50. ### Function: matrix * vector
  51. def mtrx_vc_product(a,b)
  52.   c = Array.new(N).collect{ 0.0 }
  53.   for i in (1..N)
  54.     for j in (1..N)
  55.       c[i-1] += a[i-1][j-1]*b[j-1]
  56.     end
  57.   end
  58.  
  59.   return c
  60. end
  61.  
  62. ### Function: A * vector
  63. def a_vc(b)
  64.   return mtrx_vc_product(A,b)
  65. end
  66.  
  67. ### Function: residual = b - Ax
  68. def residual(x)
  69.   return vc_sum(B, vc_coeff(-1.0, a_vc(x)))
  70. end
  71.  
  72. ### Function: The norm to use in the Jacobi method.
  73. def norm(x)
  74.   v = 0.0
  75.   for i in (1..N)
  76.     v += x[i-1]**2
  77.   end
  78.   return sqrt(v)
  79. end
  80.  
  81. ### Function: print results
  82. def print_result(x)
  83.   print "x = ["
  84.   for i in (1..N-1)
  85.     printf("%.6f",x[i-1])
  86.     print ", "
  87.   end
  88.   printf("%.6f",x[N-1])
  89.   print "], "
  90.   print "residual= ",norm(residual(x)),"\n"
  91. end
  92.  
  93. ### Function: CG iteration
  94. def cg_iteration(x,r,p,r_R,q)
  95.   q_Q = inner_product(p,q)
  96.   w = r_R/q_Q
  97.   x_new = vc_sum(x, vc_coeff(w,p))
  98.   r_new = vc_sum(r, vc_coeff((-1.0*w),q))
  99.  
  100.   r_R_new = inner_product(r_new, r_new)
  101.   beta = r_R_new / r_R
  102.   p_new = vc_sum(r_new, vc_coeff(beta,p))
  103.   q_new = a_vc(p_new)
  104.  
  105.   return [x_new, r_new, p_new, r_R_new, q_new]
  106. end
  107.  
  108. ##### main program #####
  109. # initial values
  110. x = Array.new(N).collect{ 1.0 }
  111. r = residual(x)
  112. p = vc_coeff(1.0, r)
  113. r_R = inner_product(r,r)
  114. q = a_vc(p)
  115.  
  116. print_result(x)
  117.  
  118. # loop to obtain a good x.
  119. loop_times = 0
  120. loop do
  121.   if (norm(residual(x)) < Eps) then
  122.       break
  123.     else
  124.       x,r,p,r_R,q = cg_iteration(x,r,p,r_R,q)
  125.       print_result(x)
  126.       loop_times += 1
  127.     end
  128. end
  129.  
  130. # result
  131. print "\n\nYou obtain a numerical solution after ",loop_times," loops.\n\n"
  132. print_result(x)
  133. print "\n"