计算一个144×144的矩阵，使用多个循环。

huangapple
go评论95阅读模式
英文:

Calculating a matrix of 144x144 using multiple loops

问题

I have a problem calculating the kuu_s matrix, which is a size of 144x144; the program is only calculating 121x121 elements leaving the first row and columns to 0, and the last 22x22 elements also resulting in 0.

  1. import numpy as np
  2. import scipy.integrate as spi
  4. # Define the integrand function
  5. def f0(x,n):
  6.     if n == 1:
  7.         return (1-x)
  8.     elif n == 2:
  9.         return x
  10.     else:
  11.         return np.sin((n - 2) * np.pi * x)
  13. def f1(x,n):
  14.     if n == 1:
  15.         return -1
  16.     elif n == 2:
  17.         return 1
  18.     else:
  19.         return np.cos((n - 2)*np.pi*x)*(n-2)*np.pi
  21. def f00(x, m, n):
  22.         return f0(x,m)*f0(x,n)
  24. def f11(x, m, n):
  25.     return f1(x,m)*f1(x,n)
  27. def f10(x, m, n):
  28.         return f1(x,m)*f0(x,n)
  29. def f01(x, m, n):
  30.         return f0(x,m)*f1(x,n)
  31. # Define the dimensions of the tensor
  32. NMHT_0 = 13 # Updated to 13
  34. # Create the tensor
  35. I00 = I11 = I10 = I01 = np.zeros((NMHT_0, NMHT_0))
  36. # Perform the integration
  37. for m in range(1, NMHT_0 ):
  38.     for n in range(1, NMHT_0 ):
  39.         result1, error = spi.quad_vec(f00, 0, 1, args=(m, n)) 
  40.         result2, error2 = spi.quad_vec(f11, 0, 1, args=(m, n)) 
  41.         result3, error3 = spi.quad_vec(f10, 0, 1, args=(m, n)) 
  42.         result4, error4 = spi.quad_vec(f01, 0, 1, args=(m, n))
  44.         I00[m, n] = result1 
  45.         I11[m, n] = result2  
  46.         I10[m, n] = result3 
  47.         I01[m, n] = result4
  48. a = 1
  49. b = 1
  50. A = np.random.rand(5,5)
  52. Kuu_s = np.zeros((NMHT_0 ** 2, NMHT_0 ** 2))
  54. k = 0
  55. for i in range(1, NMHT_0):
  56.      for j in range(1, NMHT_0):
  57.         q = 0
  58.         k += 1
  59.         for m in range(1, NMHT_0):
  60.             for n in range(1, NMHT_0):
  61.                 q += 1
  62.                 Kuu_s[k, q] += (
  63.                             (A[1, 1] * (1 / a ** 2) * (I11[i, m] * I00[j, n])) +
  64.                             (A[1, 3] * (1 / (a * b)) * (I10[i, m] * I01[j, n] + I01[i, m] * I10[j, n])) +
  65.                             (A[3, 3] * (1 / b ** 2) * (I00[j, n] * I11[i, m]))
  66.                     )
  68. np.savetxt("matrixkuu2.docx", Kuu_s, delimiter='  ', fmt='%1.3e')

You should update NMHT_0 to 13 instead of 12 to ensure the program calculates all the 144x144 elements.

英文:

I have a problem calculating the kuu_s matrix ,which is a size of 144x144; the program is only calculating 121x121 elements leaving the first row and columns to 0, and the last 22x22 elements also resulting in 0.

  1. import numpy as np
  2. import scipy.integrate as spi
  3. # Define the integrand function
  4. def f0(x,n):
  5. if n == 1:
  6. return (1-x)
  7. elif n == 2:
  8. return x
  9. else:
  10. return np.sin((n - 2) * np.pi * x)
  11. def f1(x,n):
  12. if n == 1:
  13. return -1
  14. elif n == 2:
  15. return 1
  16. else:
  17. return np.cos((n - 2)*np.pi*x)*(n-2)*np.pi
  18. def f00(x, m, n):
  19. return f0(x,m)*f0(x,n)
  20. def f11(x, m, n):
  21. return f1(x,m)*f1(x,n)
  22. def f10(x, m, n):
  23. return f1(x,m)*f0(x,n)
  24. def f01(x, m, n):
  25. return f0(x,m)*f1(x,n)
  26. # Define the dimensions of the tensor
  27. NMHT_0 = 12
  28. # Create the tensor
  29. I00 = I11 = I10 = I01 = np.zeros((NMHT_0, NMHT_0))
  30. # Perform the integration
  31. for m in range(1, NMHT_0 ):
  32. for n in range(1, NMHT_0 ):
  33. result1, error = spi.quad_vec(f00, 0, 1, args=(m, n)) result2, error2 = spi.quad_vec(f11, 0, 1, args=(m, n)) result3, error3 = spi.quad_vec(f10, 0, 1, args=(m, n)) result4, error4 = spi.quad_vec(f01, 0, 1, args=(m, n))
  34. I00[m, n] = result1 I11[m, n] = result2  I10[m, n] = result3 I01[m, n] = result4
  35. a = 1
  36. b = 1
  37. A = np.random.rand(5,5)
  38. Kuu_s = np.zeros((NMHT_0 ** 2, NMHT_0 ** 2))
  39. k = 0
  40. for i in range(1, NMHT_0):
  41. for j in range(1, NMHT_0):
  42. q = 0
  43. k += 1
  44. for m in range(1, NMHT_0):
  45. for n in range(1, NMHT_0):
  46. q += 1
  47. Kuu_s[k, q] += (
  48. (A[1, 1] * (1 / a ** 2) * (I11[i, m] * I00[j, n])) +
  49. (A[1, 3] * (1 / (a * b)) * (I10[i, m] * I01[j, n] + I01[i, m] * I10[j, n])) +
  50. (A[3, 3] * (1 / b ** 2) * (I00[j, n] * I11[i, m]))
  51. )
  52. np.savetxt("matrixkuu2.docx",Kuu_s,delimiter='  ',fmt='%1.3e')

the following pic is from mathcad
计算一个144×144的矩阵，使用多个循环。

I've tried changing the limits of m n i j ,but I get the error of out-of-bound index.
How could I fix the issue so that the program calculates all the 144x144 elements and adds kuu_s to each loop

答案1

得分: 2

我已修复了您的代码中的间距问题，以及评论中提出的问题。以下是结果。它确实生成一个144x144的输出数组。

  1. import numpy as np
  2. import scipy.integrate as spi
  4. # 定义被积函数
  5. def f0(x, n):
  6.     if n == 1:
  7.         return (1 - x)
  8.     elif n == 2:
  9.         return x
  10.     else:
  11.         return np.sin((n - 2) * np.pi * x)
  13. def f1(x, n):
  14.     if n == 1:
  15.         return -1
  16.     elif n == 2:
  17.         return 1
  18.     else:
  19.         return np.cos((n - 2) * np.pi * x) * (n - 2) * np.pi
  21. def f00(x, m, n):
  22.     return f0(x, m) * f0(x, n)
  24. def f11(x, m, n):
  25.     return f1(x, m) * f1(x, n)
  27. def f10(x, m, n):
  28.     return f1(x, m) * f0(x, n)
  30. def f01(x, m, n):
  31.     return f0(x, m) * f1(x, n)
  33. # 定义张量的维度
  34. NMHT_0 = 12
  36. # 创建张量
  37. I00 = np.zeros((NMHT_0, NMHT_0))
  38. I11 = np.zeros((NMHT_0, NMHT_0))
  39. I10 = np.zeros((NMHT_0, NMHT_0))
  40. I01 = np.zeros((NMHT_0, NMHT_0))
  42. # 执行积分
  43. for m in range(NMHT_0):
  44.     for n in range(NMHT_0):
  45.         result1, error1 = spi.quad_vec(f00, 0, 1, args=(m, n))
  46.         result2, error2 = spi.quad_vec(f11, 0, 1, args=(m, n))
  47.         result3, error3 = spi.quad_vec(f10, 0, 1, args=(m, n))
  48.         result4, error4 = spi.quad_vec(f01, 0, 1, args=(m, n))
  50.         I00[m, n] = result1
  51.         I11[m, n] = result2
  52.         I10[m, n] = result3
  53.         I01[m, n] = result4
  55. a = 1
  56. b = 1
  57. A = np.random.rand(5, 5)
  59. Kuu_s = np.zeros((NMHT_0 ** 2, NMHT_0 ** 2))
  61. for i in range(NMHT_0):
  62.     for j in range(NMHT_0):
  63.         k = i * NMHT_0 + j
  64.         for m in range(NMHT_0):
  65.             for n in range(NMHT_0):
  66.                 q = m * NMHT_0 + n
  67.                 Kuu_s[k, q] += (
  68.                     (A[1, 1] * (1 / a ** 2) * (I11[i, m] * I00[j, n])) +
  69.                     (A[1, 3] * (1 / (a * b)) * (I10[i, m] * I01[j, n] + I01[i, m] * I10[j, n])) +
  70.                     (A[3, 3] * (1 / b ** 2) * (I00[j, n] * I11[i, m]))
  71.                 )
  73. np.savetxt("matrixkuu2.txt", Kuu_s, delimiter='  ', fmt='%1.3e')
英文:

I've fixed the spacing problems with your code, plus the issues raised in the comments. This is the result. It does produce a 144x144 output array.

  1. import numpy as np
  2. import scipy.integrate as spi
  3. # Define the integrand function
  4. def f0(x,n):
  5. if n == 1:
  6. return (1-x)
  7. elif n == 2:
  8. return x
  9. else:
  10. return np.sin((n - 2) * np.pi * x)
  11. def f1(x,n):
  12. if n == 1:
  13. return -1
  14. elif n == 2:
  15. return 1
  16. else:
  17. return np.cos((n - 2)*np.pi*x)*(n-2)*np.pi
  18. def f00(x, m, n):
  19. return f0(x,m)*f0(x,n)
  20. def f11(x, m, n):
  21. return f1(x,m)*f1(x,n)
  22. def f10(x, m, n):
  23. return f1(x,m)*f0(x,n)
  24. def f01(x, m, n):
  25. return f0(x,m)*f1(x,n)
  26. # Define the dimensions of the tensor
  27. NMHT_0 = 12
  28. # Create the tensor
  29. I00 = np.zeros((NMHT_0, NMHT_0))
  30. I11 = np.zeros((NMHT_0, NMHT_0))
  31. I10 = np.zeros((NMHT_0, NMHT_0))
  32. I01 = np.zeros((NMHT_0, NMHT_0))
  33. # Perform the integration
  34. for m in range(NMHT_0):
  35. for n in range(NMHT_0):
  36. result1, error1 = spi.quad_vec(f00, 0, 1, args=(m, n))
  37. result2, error2 = spi.quad_vec(f11, 0, 1, args=(m, n))
  38. result3, error3 = spi.quad_vec(f10, 0, 1, args=(m, n))
  39. result4, error4 = spi.quad_vec(f01, 0, 1, args=(m, n))
  40. I00[m, n] = result1 
  41. I11[m, n] = result2  
  42. I10[m, n] = result3 
  43. I01[m, n] = result4
  44. a = 1
  45. b = 1
  46. A = np.random.rand(5,5)
  47. Kuu_s = np.zeros((NMHT_0 ** 2, NMHT_0 ** 2))
  48. for i in range(NMHT_0):
  49. for j in range(NMHT_0):
  50. k = i * NMHT_0 + j
  51. for m in range(NMHT_0):
  52. for n in range(NMHT_0):
  53. q = m * NMHT_0 + n
  54. Kuu_s[k, q] += (
  55. (A[1, 1] * (1 / a ** 2) * (I11[i, m] * I00[j, n])) +
  56. (A[1, 3] * (1 / (a * b)) * (I10[i, m] * I01[j, n] + I01[i, m] * I10[j, n])) +
  57. (A[3, 3] * (1 / b ** 2) * (I00[j, n] * I11[i, m]))
  58. )
  59. np.savetxt("matrixkuu2.txt",Kuu_s,delimiter='  ',fmt='%1.3e')

答案2

得分: 1

在创建Kuu_s时，您的循环从1开始，但应该从0开始，因为Python是从零开始索引的。第二个问题是关于起始k和q值的问题。将它们从0开始并在循环开始时递增意味着它们首次使用时会等于1，但（如我已经提到的）Python是从零开始索引的，所以您会错过第一行和列。您有两个选择：

  1. 将k和q从0开始并在循环结束时递增它们。
  2. 将k和q从-1开始并在循环开始时递增它们。

我选择使用第二个选项，因为我认为这种方式的代码更清晰，尽管它可能更难阅读。

  1. import numpy as np
  2. from scipy.integrate import quad_vec
  4. # 省略部分代码...
  6. NMHT_0 = 12
  8. M, N = np.meshgrid(np.arange(0, NMHT_0), np.arange(0, NMHT_0))
  9. I00, _ = quad_vec(f00, 0, 1, args=(M, N))
  10. I11, _ = quad_vec(f11, 0, 1, args=(M, N))
  11. I10, _ = quad_vec(f10, 0, 1, args=(M, N))
  12. I01, _ = quad_vec(f01, 0, 1, args=(M, N))
  14. a = 1
  15. b = 1
  16. A = np.random.rand(5, 5)
  18. Kuu_s = np.zeros((NMHT_0**2, NMHT_0**2))
  20. k = -1
  21. for i in range(NMHT_0):
  22.     for j in range(NMHT_0):
  23.         k += 1
  24.         q = -1
  25.         for m in range(NMHT_0):
  26.             for n in range(NMHT_0):
  27.                 q += 1
  28.                 Kuu_s[k, q] = (
  29.                     A[1, 1]*I11[i, m]*I00[j, n]/a**2
  30.                     + A[3, 3]*I00[i, m]*I11[j, n]/b**2
  31.                     + A[1, 3]*(I10[i, m]*I01[j, n]
  32.                                + I01[i, m]*I10[j, n])/a/b
  33.                 )
  35. np.savetxt("matrixkuu2.txt", Kuu_s, delimiter="  ", fmt="%1.3e")

您还会注意到，我取消了最初用于积分的双重for循环，而是选择正确使用quad_vec来计算积分。可能有一种方法可以摆脱末尾的大型嵌套循环，但我没有费心去做。只有在NHMT_0变得非常大时才会成为问题（即会很慢）。

英文:

When creating Kuu_s, you start your loops at 1, when they should start at 0 because Python is zero-indexed. The second issue is starting your k and q values. Starting them at 0 and incrementing at the start of the loop means they will first be used when they equal 1, but (as I already mentioned) Pyhton is zero-indexed, so you're missing the first row and column. You have two options:

  1. Start k and q at 0 and increment them at the end of the loop.
  2. Start k and q at -1 and increment them at the start of the loop.

I chose to use the second option because I think the code is cleaner that way, though it is arguably less readable.

  1. import numpy as np
  2. from scipy.integrate import quad_vec
  5. def f0(x, N):
  6.     res = np.sin((N - 2)*np.pi*x)
  7.     res[N == 1] = (1-x)
  8.     res[N == 2] = x
  10.     return res
  13. def f1(x, N):
  14.     res = np.cos((N - 2)*np.pi*x)*(N-2)*np.pi
  15.     res[N == 1] = -1
  16.     res[N == 2] = 1
  18.     return res
  21. def f00(x, M, N):
  22.     return f0(x, M)*f0(x, N)
  25. def f11(x, M, N):
  26.     return f1(x, M)*f1(x, N)
  29. def f10(x, M, N):
  30.     return f1(x, M)*f0(x, N)
  33. def f01(x, M, N):
  34.     return f0(x, M)*f1(x, N)
  37. NMHT_0 = 12
  39. M, N = np.meshgrid(np.arange(0, NMHT_0), np.arange(0, NMHT_0))
  40. I00, _ = quad_vec(f00, 0, 1, args=(M, N))
  41. I11, _ = quad_vec(f11, 0, 1, args=(M, N))
  42. I10, _ = quad_vec(f10, 0, 1, args=(M, N))
  43. I01, _ = quad_vec(f01, 0, 1, args=(M, N))
  46. a = 1
  47. b = 1
  48. A = np.random.rand(5, 5)
  50. Kuu_s = np.zeros((NMHT_0**2, NMHT_0**2))
  52. k = -1
  53. for i in range(NMHT_0):
  54.     for j in range(NMHT_0):
  55.         k += 1
  56.         q = -1
  57.         for m in range(NMHT_0):
  58.             for n in range(NMHT_0):
  59.                 q += 1
  60.                 Kuu_s[k, q] = (
  61.                     A[1, 1]*I11[i, m]*I00[j, n]/a**2
  62.                     + A[3, 3]*I00[i, m]*I11[j, n]/b**2
  63.                     + A[1, 3]*(I10[i, m]*I01[j, n]
  64.                                + I01[i, m]*I10[j, n])/a/b
  65.                 )
  67. np.savetxt("matrixkuu2.txt", Kuu_s, delimiter="  ", fmt="%1.3e")

You'll also note that I got rid of the initial double for loops used for integration in favor of properly using quad_vec for evaluating the integrals. There might be a way to get rid of the large nested loop at the end, but I didn't bother. It will only be a problem (i.e. it will be slow) if NHMT_0 gets very large.

huangapple
  • 本文由 发表于 2023年6月26日 14:24:27
  • 转载请务必保留本文链接：https://go.coder-hub.com/76554006.html
  • numpy
  • python
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定