Part 6 -- Solving Simultaneous Equations in Python¶
In this exercise we will be solving our truss problem both as a matrix algebra problem and as a generic simultaneous equation question.
The text parts of this exercise also use a lot of LaTeX Equations. An online tool to help with make these equations can be found at the address below.
################################################
#
# Libraries
#
import numpy as np
import scipy.optimize as opt
import sympy as sym
#
################################################
Problem 1 : Our Truss Case¶
Let's explore our basic truss problem:
Consider the system below. Solve for $F_1$, $F_2$, $F_3$, $H_2$, $V_2$, and $V_3$. Correct answers in red.
Problem Setup¶
Let's start by setting up the constants.
Basic Equations¶
Let's set up the basic balance equations.
We'll have six equations. Since we have sines and cosines in the formulae we will need two sets of equations if we want to play with the symbolic and SciPy/NumPy solvers.
$$f_{1x}\left( F_1,F_3 \right) = -F_1 \cos \theta_2 + F_3 \cos \theta_3 + P_{1x}$$$$f_{1y}\left( F_1,F_3 \right) = -F_1 \sin \theta_2 - F_3 \sin \theta_3 + P_{1y}$$$$f_{2x}\left( F_1,F_2,H_2 \right) = F_1 \cos \theta_2 + F_2 + H_2 + P_{2x}$$$$f_{2y}\left( F_1,V_2 \right) = F_1 \sin \theta_2 + V_2 + P_{2y}$$$$f_{3x}\left( F_2,F_3 \right) = -F_2 - F_3 \cos \theta_3 + P_{3x}$$$$f_{3y}\left( F_3,V_3 \right) = F_3 \sin \theta_3 + V_3 + P_{3y}$$Matrix Setup¶
For the matrix setup the above equations translates into the linear algebra system
(also Fun Fact: Doing Matricies in LaTeX Markdown is even worse than doing it in MS Word's new Equation Editor even with the Web Tool!)
$$\mathbf{F} \vec{f}=-\vec{p}$$$$\begin{bmatrix} -\cos\theta_2 & 0 & \cos\theta_3 & 0 & 0 & 0 \\ -\sin\theta_2 & 0 & -\sin\theta_3 & 0 & 0 & 0 \\ \cos\theta_2 & 1 & 0 & 1 & 0 & 0 \\ \sin\theta_2 & 0 & 0 & 0 & 1 & 0 \\ 0 & -1 & -\cos\theta_3 & 0 & 0 & 0 \\ 0 & 0 & \sin\theta_3 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} F_1\\ F_2\\ F_3\\ H_2\\ V_2\\ V_3 \end{bmatrix} = \begin{bmatrix} -P_{1x}\\ -P_{1y}\\ -P_{2x}\\ -P_{2y}\\ -P_{3x}\\ -P_{3y}\ \end{bmatrix}$$Constants and Parameters¶
Let's load up our known values. Notice that handy function that comes with NumPy to convert Degrees to Radians.
################################################
#
# Constants and Parameters
#
theta_2 = np.radians(30.0)
theta_3 = np.radians(60.0)
P1x = 0.
P1y = -1000.
P2x = 0.
P2y = 0.
P3x = 0.
P3y = 0.
#
################################################
Treating the System as a System of Simutaneous Equations¶
Let's start by creating a symbolic solution system (we'll prefix these functions with a "g")
################################################
#
# SymPy Formulas for Our System
#
#
# Create our "symbolic variables"
#
F1 = sym.symbols("F_1")
F2 = sym.symbols("F_2")
F3 = sym.symbols("F_3")
H2 = sym.symbols("H_2")
V2 = sym.symbols("V_2")
V3 = sym.symbols("V_3")
def g1x(F1, F3):
return( - F1 * sym.cos(theta_2) + F3 * sym.cos(theta_3) + P1x )
def g1y(F1, F3):
return( - F1 * sym.sin(theta_2) - F3 * sym.sin(theta_3) + P1y )
def g2x(F1, F2, H2):
return( F1 * sym.cos(theta_2) + F2 + H2 + P2x )
def g2y(F1, V2):
return( F1 * sym.sin(theta_2) + V2 + P2y )
def g3x(F2, F3):
return( -F2 - F3*sym.cos(theta_3) + P3x )
def g3y(F3, V3):
return( F3*sym.cos(theta_3) + V3 + P3y )
print()
print(">> Symbolic Equation System")
print()
display(g1x(F1, F3))
display(g1y(F1, F3))
print()
display(g2x(F1, F2, H2))
display(g2y(F1, V2))
print()
display(g3x(F2, F3))
display(g3y(F3, V3))
#
################################################
SymPySympy Solution for Simultaneous Equations¶
Early versions of SymPy would recommend the generic solve() function in SymPy. However this is being depreciated in favor of other resources. The justification for some of the issues and it's replacement sub-package is called solveset.
But for legacy purposes, here is the sympy.solve function.
################################################
#
# Solving Symbolic Equations in Sympy
#
solve_eq_sympy = sym.solve([g1x(F1, F3),
g1y(F1, F3),
g2x(F1, F2, H2),
g2y(F1, V2),
g3x(F2, F3),
g3y(F3, V3)],
(F1, F2, F3, H2, V2, V3))
display(solve_eq_sympy)
#
################################################
But here are two examples of working with the solveset packages.
First we have the recommended solver for this type of system: linsolve
################################################
#
# Solving Symbolic Equations in SymPy with linsolve
#
linsolve_eq_sympy = sym.linsolve([g1x(F1, F3),
g1y(F1, F3),
g2x(F1, F2, H2),
g2y(F1, V2),
g3x(F2, F3),
g3y(F3, V3)],
(F1, F2, F3, H2, V2, V3))
print(linsolve_eq_sympy)
#
################################################
Second we have the recommended solver for a nonlinear system: nonlinsolve just to demonstate the function if we had nonlinear terms in our system for future reference.
################################################
#
# Solving Symbolic Equations in SymPy with nonlinsolve
#
nonlinsolve_eq_sympy = sym.nonlinsolve([g1x(F1, F3),
g1y(F1, F3),
g2x(F1, F2, H2),
g2y(F1, V2),
g3x(F2, F3),
g3y(F3, V3)],
(F1, F2, F3, H2, V2, V3))
print(nonlinsolve_eq_sympy)
#
################################################
################################################
#
# Solving Symbolic Equations in SymPy with nonlinsolve
#
F = sym.Matrix([[-sym.cos(theta_2), 0, sym.cos(theta_3), 0, 0, 0],
[-sym.sin(theta_2), 0, -sym.sin(theta_3), 0, 0, 0],
[ sym.cos(theta_2), 1, 0 , 1, 0, 0],
[ sym.sin(theta_2), 0, 0 , 0, 1, 0],
[ 0 , -1, -sym.cos(theta_3), 0, 0, 0],
[ 0 , 0, sym.sin(theta_3), 0, 0, 1]])
neg_p = sym.Matrix([[-P1x],
[-P1y],
[-P2x],
[-P2y],
[-P3x],
[-P3y]])
display(F)
display(neg_p)
#
################################################
If we want to do a mimalist matrix algebra solution it's rather easy. We just invert F and multiply it against -$\vec{p}$. Inversion, Determinant and other Matrix operators and SymPy. Can be found inside the Matrix object in SymPy.
################################################
#
# Solving Symbolic Linear Algebra Matrices w/ SymPy.
#
# The matrix object in SymPy comes with objects
# for basic matrix operations
#
inverted_matrix_sympy = F.inv() * neg_p
display(inverted_matrix_sympy)
#
################################################
We can also use the linsolve function: linsolve
################################################
#
# Solving Symbolic Matrix for Linear Equations in SymPy with linsolve
#
linsolve_matrix_sympy = sym.linsolve((F, neg_p),
(F1, F2, F3, H2, V2, V3))
print(linsolve_matrix_sympy)
#
################################################
Simultaneous Equations with SciPy's fslolve function¶
Now let's do the above problem with the SciPy libraries.
As with our root exercise a function that has a similar functionality to our Mathcad Solve Block construct us the fsolve function.
Let's first rewrite our equations here (we'll use NumPy trig functions) We also need to cluster them as a single system.
This first appraoch will work with nonlinear equations as well just as it does in a solve block.
################################################
#
# SciPy/NumPy Formulas for Our System
#
#
# We create a single function for all of our equations
# there will be a single argument. It will expect
# in this case an array 6 slots long which will
# hold our F1, F2, F3, H2, V2, and F3.
#
# Holding true to our earlier nomenclature, let's call this new
# function of multiple g#x's capital-G
def G(f):
#
# firsrt, let's give the input vector, f
# parameter names to which we can relate.
#
F1 = f[0]
F2 = f[1]
F3 = f[2]
H2 = f[3]
V2 = f[4]
V3 = f[5]
#
# Now we create an empty 6-long vector
#
system = np.empty(6)
#
# And we fill that fector with our equations
#
system[0] = - F1 * np.cos(theta_2) + F3 * np.cos(theta_3) + P1x
system[1] = - F1 * np.sin(theta_2) - F3 * np.sin(theta_3) + P1y
system[2] = F1 * np.cos(theta_2) + F2 + H2 + P2x
system[3] = F1 * np.sin(theta_2) + V2 + P2y
system[4] = -F2 - F3 * np.cos(theta_3) + P3x
system[5] = F3 * np.cos(theta_3) + V3 + P3y
#
# We're done so let's return the variable "system"
#
return(system)
#
################################################
And now we create a first guess of a 6-zero-long vector (or any abitrary number) and feed it into our fsolve function.
################################################
#
# SciPy's FSolve function.
#
# F1 F2,F3,H2,V2,V3
first_guess = np.array(( 0, 0, 0, 0, 0, 0))
f_root_scalar_fsolve = opt.fsolve(G, first_guess)
print()
print(">>> Fsolve")
print(">>>>> F1 = ", f_root_scalar_fsolve[0])
print(">>>>> F2 = ", f_root_scalar_fsolve[1])
print(">>>>> F3 = ", f_root_scalar_fsolve[2])
print(">>>>> H2 = ", f_root_scalar_fsolve[3])
print(">>>>> V2 = ", f_root_scalar_fsolve[4])
print(">>>>> V3 = ", f_root_scalar_fsolve[5])
print()
print(">>> Test : G(f)", G(f_root_scalar_fsolve))
#
################################################
Matrix Approach with NumPy¶
Let's construct a simple pair of matricies with NumPy's matrix object. Like the SymPy variant, this will have a number of attributes that will help you calculate determinants, inverses, etc.
Advanced warning: These are not going to be as pretty as it is when working with SymPy
################################################
#
# Solving Matrix Equations in NumPy
#
#
# Let's make a matrix!
#
F = np.matrix([[-np.cos(theta_2), 0, np.cos(theta_3), 0, 0, 0],
[-np.sin(theta_2), 0, -np.sin(theta_3), 0, 0, 0],
[ np.cos(theta_2), 1, 0 , 1, 0, 0],
[ np.sin(theta_2), 0, 0 , 0, 1, 0],
[ 0 , -1, -np.cos(theta_3), 0, 0, 0],
[ 0 , 0, np.sin(theta_3), 0, 0, 1]])
#
# And a vector using the matrix object.
#
neg_p = np.matrix([[-P1x],
[-P1y],
[-P2x],
[-P2y],
[-P3x],
[-P3y]])
#
# Do a fast print
#
print(">>>> F Matrix")
display(F)
print()
print(">>>> -p vector")
display(neg_p)
print()
#
# Now let's try to make an inverse.
#
# the [X].I operator will invert the matrix.
#
f_matrix_inverse = F.I * neg_p
print(">>>> f via Matrix Solution")
print()
display(f_matrix_inverse)
print(">>> Test : G(f)", G(f_matrix_inverse))
#
################################################
Version Information¶
################################################################
#
# Loading Version Information
#
%load_ext version_information
%version_information version_information, numpy, matplotlib, scipy, sympy
#
################################################################