¼ø¶È»ñÎÁ/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"