Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

Untitled

By: a guest on Nov 1st, 2011  |  syntax: Python  |  size: 1.39 KB  |  views: 32  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. import numpy as np
  2.  
  3. # input_values
  4. a = np.array([0.123456, 0.234567, 0.345678])
  5. #outut_values
  6. b = np.array([0.876543, 0.765432, 0.654321])
  7.  
  8.  
  9. # A 3x3 Matrix multiplication $X a = b$ is a set of 3 equations
  10. #
  11. #   x_11 * a_1 + x_12 * a_2 + x_13 * a_3 = b_1
  12. #   x_21 * a_1 + x_22 * a_2 + x_23 * a_3 = b_2
  13. #   x_31 * a_1 + x_32 * a_2 + x_33 * a_3 = b_3
  14. #
  15. # There we have 3*3 = 9 unknown, but only 3 equations.  This doesn't work.  We
  16. # have to cheat.  We just set 6 of the 9 unkown to 1.  For example every x_ij
  17. # but the fist column
  18. #
  19. #   x_11 * a_1 + a_2 + a_3 = b_1
  20. #   x_21 * a_1 + a_2 + a_3 = b_2
  21. #   x_31 * a_1 + a_2 + a_3 = b_3
  22. #
  23. #   x_11 = b_1 - (a_2 + a_3) / a_1
  24. #   x_21 = b_2 - (a_2 + a_3) / a_1
  25. #   x_31 = b_3 - (a_2 + a_3) / a_1
  26. #
  27.  
  28. X = np.ones((3,3))
  29. X[0, 0] = (b[0] - (a[1] + a[2])) / a[0]
  30. X[1, 0] = (b[1] - (a[1] + a[2])) / a[0]
  31. X[2, 0] = (b[2] - (a[1] + a[2])) / a[0]
  32.  
  33. assert np.all( (np.dot(X, a) - b) < 1e-15 )  # approximately the same
  34.  
  35.  
  36. # DOWNSIDE: Doesn't work if a_1 is zero or almost zero.  But this is easy to
  37. # fix be picking the index i where |a_i| is the the greatest.
  38.  
  39. i_max = np.argmax(np.abs(a))
  40. i_other = list(set(range(3)).difference([i_max]))
  41.  
  42. X = np.ones((3,3))
  43. X[0, i_max] = (b[0] - np.sum(a[i_other])) / a[i_max]
  44. X[1, i_max] = (b[1] - np.sum(a[i_other])) / a[i_max]
  45. X[2, i_max] = (b[2] - np.sum(a[i_other])) / a[i_max]
  46.  
  47. assert np.all( (np.dot(X, a) - b) < 1e-15 )  # approximately the same