I am trying to modify the code from Oliver Natt's book "Physik mit Python" to calculate a more complex structure. However, when I add points to the code, I receive the error message "numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square" when attempting to solve the equation system. The code is written in Python and uses the numpy library to solve the system. Can anyone help me understand why I am getting this error and how to fix it?
This is the code which causes the Error:
import numpy as np
#Set the scaling factor for the force vectors [m/N].
force_scale = 0.005
#Number of dimensions for the problem (2 or 3).
dim = 3
#Set the positions of the points [m].
points = np.array([
[0, 0, 0], [1.5, 0, 0], [0, 1.5, 0], [0, 0, 2], [1, 0, 2],
[1, 1, 2], [0, 1, 2], [1.5, 1.5, 0]
])
#Create a list of the indices of the support points.
support_indices = [0, 1, 2, 7]
#Each bar connects exactly two points. We store the indices
#of the corresponding points in an array.
bars = np.array([
[0, 3], [1, 4], [2, 6], [3, 4], [4, 5], [5, 6], [6, 3],
[3, 5], [0, 4], [0, 6], [0, 5], [1, 5], [7, 5]
])
#Set the external force acting on each point [N].
#For the support points, we initially set this force to 0.
#This will be calculated later.
external_force = np.array([
[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, -98.1], [0, 0, -98.1],
[0, 0, -98.1], [0, 0, -98.1], [0, 0, 0]
])
#Define the number of points, bars, etc.
n_points = points.shape[0]
n_bars = bars.shape[0]
n_support = len(support_indices)
n_nodes = n_points - n_support
n_equations = n_nodes * dim
#Create a list of the indices of the nodes.
node_indices = list(set(range(n_points)) - set(support_indices))
def unit_vector(point_index, bar_index):
"""Returns the unit vector pointing from the point with index
point_index along the bar with index bar_index. """
i1, i2 = bars[bar_index]
if point_index == i1:
vec = points[i2] - points[i1]
else:
vec = points[i1] - points[i2]
return vec / np.linalg.norm(vec)
# Set up the equation system for the forces.
A = np.zeros((n_equations, n_bars))
for i, bar in enumerate(bars):
for k in np.intersect1d(bar, node_indices):
n = node_indices.index(k)
A[n * dim:(n + 1) * dim, i] = unit_vector(k, i)
#Solve the equation system A @ F = -external_force for the forces F.
b = -external_force[node_indices].reshape(-1)
forces = np.linalg.solve(A, b)
#Calculate the external forces.
for i, bar in enumerate(bars):
for k in np.intersect1d(bar, support_indices):
external_force[k] -= forces[i] * unit_vector(k, i)
If I use the following arrays, then the program has no problems
points = np.array([
[0, 0, 0], [1.5, 0, 0], [0, 1.5, 0], [0, 0, 2], [1, 0, 2],
[1, 1, 2], [0, 1, 2]
])
support_indices = [0, 1, 2]
bars = np.array([
[0, 3], [1, 4], [2, 6], [3, 4], [4, 5], [5, 6], [6, 3],
[3, 5], [0, 4], [0, 6], [0, 5], [1, 5]
])
external_force = np.array([
[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, -98.1], [0, 0, -98.1],
[0, 0, -98.1], [0, 0, -98.1]
])
This is the traceback I received:
Traceback (most recent call last):
File "C:\...\main.py", line 62, in <module>
forces = np.linalg.solve(A, b)
File "<__array_function__ internals>", line 200, in solve
File "C:\...\venv\lib\site-packages\numpy\linalg\linalg.py", line 373, in solve
_assert_stacked_square(a)
File "C:\...\venv\lib\site-packages\numpy\linalg\linalg.py", line 190, in _assert_stacked_square
raise LinAlgError('Last 2 dimensions of the array must be square')
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square
If I copy and paste the original code from Oliver Natt's book "Physik mit Python" and modify the points and bars to calculate a more complex structure, I get a different traceback error compared to the code above:
Traceback (most recent call last):
File "C:\...\original.py", line 66, in <module>
A[n * dim:(n + 1) * dim, i] = einheitsvektor(k, i)
IndexError: index 12 is out of bounds for axis 1 with size 12
I have tried several steps to debug the code, including checking for syntax errors, re-arranging the code, and verifying the mathematical calculations. However, I am still encountering errors and suspect that there may be a mathematical issue that I am currently unable to identify.
Additionally, it is quite surprising to me that when I try to write the same code in two different programming languages, I get completely different Tracebacks. This leads me to believe that the issue might not just be a simple syntax error, but perhaps there is something deeper going on that I am currently unable to see.
