Mathematik für Ingenieure mit Python

Mathematische Grundlagen
www.grund-wissen.de/mathematik
WinPython
Alles was man für Mathematik mit Python unter Windows braucht.
winpython.sourceforge.net
Download WinPython

Python-Code
vollständiger Python-Code: fachwerk_2d_a3_solve.py
Was hat's gebracht?
Wir haben gesehn, dass der Code zum Aufstellen des Gleichungssystems Gemeinsamkeiten mit dem Code zum Erzeugen der Abbildung für die Aufgabenstellung hat (fachwer_2d_a1_plot.py). In beiden Skripten wird das Fachwerk durch die Arrays nodes und bars beschrieben. Die Abbildung definiert die Aufgabe für den Menschen und das Gleichungssystem definiert die Aufgabe für den Rechner.
Wie geht's weiter?
Wenn wir ein kleines Fachwerk automatisiert rechnen können, dann können wir auch große Fachwerke rechnen.
Fragen & Antworten
Python Forum:
Mathematik für Ingenieure mit Python: Fachwerk

Einfaches zweidimensionales Fachwerk, mit Numpy automatisiert gerechnet

Das Gleichungssystem mit Python und Numpy aufstellen und lösen

Der Hauptaufwand bei Fachwerken ist nicht das Lösen des Gleichungssystems, sondern das Aufstellen des Gleichungssystems. Deshalb soll hier nun das Aufstellen des Gleichungssystems mit Python automatisiert werden.
import numpy as np

###############################################################################
def unitVector(v):
    """
    Richtungsvektor berechnen.
    """
    v = np.asanyarray(v)
    return v/np.sqrt((v**2).sum())

###############################################################################
# Das Fachwerk
# index: nodeIndex
# value: (x,y)
nodes = np.array([[6,7],  # A
                  [2,7],  # B
                  [2,4],  # C
                  [0,4],  # D
                  [2,0],  # E
                  [0,0]], # F
                 dtype=np.float64)

# index: barIndex
# value: (node1, node2)
bars = np.array([[0,1],  # S1
                 [0,2],  # S2
                 [1,2],  # S3
                 [1,3],  # S4
                 [2,3],  # S5
                 [2,4],  # S6
                 [3,4],  # S7
                 [3,5],  # S8
                 [4,5]], # S9
                dtype=np.int32)

# index: nodeIndex, value: connected to (barIndex, nodeIndex)
nodeToBar = [[] for i in range(nodes.shape[0])]
for barIndex, (n1,n2) in enumerate(bars):
    nodeToBar[n1].append( (barIndex, n2))
    nodeToBar[n2].append( (barIndex, n1))
    

###############################################################################
# Externe Kraefte (F0 und Lagerkraefte)
# index: forceIndex
# value: (node, direction, value), None heisst unbekannte Lagerkraft)
externalForces = [[0, [0,-1], 1.0],  # F0
                  [4, [1, 0], None], # F_Ex
                  [4, [0, 1], None], # F_Ey
                  [5, [0, 1], None]] # F_Fy
                  

###############################################################################
# Das Gleichungssystem: A*x = b
# Fast alle Koeffizienten des Gleichungssystems sind 0,
# deshalb werden A und b mit 0 vorbelegt.
A = np.zeros((12,12), dtype=np.float64)
b = np.zeros((12,),   dtype=np.float64)

#Fuer die Ausgabe und als Doku: welcher Index gehoert zu welcher Groesse?
xNames = ['S1','S2','S3','S4','S5','S6','S7','S8','S9','F_Ex','F_Ey','F_Fy']

for n1, barIndexes in enumerate(nodeToBar):
    for barIndex, n2 in barIndexes:
        nx, ny = unitVector(nodes[n2]-nodes[n1])
        A[2*n1,   barIndex] = nx
        A[2*n1+1, barIndex] = ny

forceIndex = bars.shape[0]
for nodeIndex, direction, value in externalForces:
    nx, ny = unitVector(direction)
    
    # unbekannte Lagerkraft
    if value is None: 
        A[2*nodeIndex,   forceIndex] = nx
        A[2*nodeIndex+1, forceIndex] = ny
        forceIndex += 1
        
    # gegebene externe Kraft
    else:             
        b[2*nodeIndex]   = -nx*value
        b[2*nodeIndex+1] = -ny*value
        
# Loesen des Gleichungssystems mit Numpy
x = np.linalg.solve(A, b)

# Ausgabe der Loesung
for xName, xValue in zip(xNames, x):
    print('%4s = % f' % (xName, xValue))

Ausgabe

  S1 =  1.333333
  S2 = -1.666667
  S3 = -2.000000
  S4 =  2.403701
  S5 = -1.333333
  S6 = -3.000000
  S7 =  0.000000
  S8 =  2.000000
  S9 =  0.000000
F_Ex =  0.000000
F_Ey =  3.000000
F_Fy = -2.000000