Solving Structural Analysis Problems with Linear Algebra (Matrix Method) in R and GNU Octave
Irucka Embry, E.I.T. (EcoC2S)
2023-11-27
Replicate the R code
Install the packages
Note: If you wish to replicate the R code below, then you will need to copy and paste the following commands in R first to make sure you have all the packages and their dependencies installed:
install.packages(c("install.load", "import", "ramify", "data.table", "pander", "iemisc"))
# install the packages and their dependencies
Load the packages
Note: If you wish to replicate the R code below, then you will need to copy and paste the following commands in R first to make sure you have all the packages loaded:
install.load::load_package("ramify", "data.table", "pander", "iemisc")
# load needed packages using the load_package function from the install.load
# package (it is assumed that you have already installed these packages)
# set the pander options
panderOptions("table.continues", "")
panderOptions("table.caption.prefix", "")
panderOptions("missing", "")
import::from(pracma, mldivide)
# import mldivide from the pracma package
Example Problem 3.2 (Hibbeler) solved using R
# Solves the system of equations AX = B, where X is a column vector of unknown
# forces, of the form X [FAG FAB FGB FGF FBF FBC]'
# Coefficient matrix
Aa <- mat("sind(30), 0, 0, 0, 0, 0; 0, 1, 0, 0, 0, 0; 0, 0, sind(60), 0, 0, 0; 0, 0, 0, 1, 0, 0; 0, 0, 0, 0, sind(60), 0; 0, 0, 0, 0, 0, 1",
eval = TRUE)
pander(Aa)| 0.5 | 0 | 0 | 0 | 0 | 0 |
| 0 | 1 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0.866 | 0 | 0 | 0 |
| 0 | 0 | 0 | 1 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0.866 | 0 |
| 0 | 0 | 0 | 0 | 0 | 1 |
# Right-hand side, external forces
Bb <- mat("4; 8 * cosd(30); 3 * cosd(30); 8 - 3 * sind(30) - 3.00 * cosd(60); 3.00 * sind(30); -1.73 * cosd(60) - 3.00 * cosd(30) + 6.928",
eval = TRUE)
pander(Bb)| 4 |
| 6.928 |
| 2.598 |
| 5 |
| 1.5 |
| 3.465 |
# Solve for the unknown forces
Xx <- solve(Aa, Bb)
# Use base R's built-in solve function for solving the system of equations
pander(Xx)| 8 |
| 6.928 |
| 3 |
| 5 |
| 1.732 |
| 3.465 |
# or
Xxx <- mldivide(Aa, Bb)
# Use pracma's mldivide (which is the same as the GNU Octave built-in \
# operator for solving the system of equations)
pander(Xxx)| 8 |
| 6.928 |
| 3 |
| 5 |
| 1.732 |
| 3.465 |
# Display the results in a data.table
Xdisplay <- data.table(Forces = c("FAG", "FAB", "FGB", "FGF", "FBF", "FBC"), Value = Xxx,
TC = c("C", "T", "C", "C", "T", "T"))
# set the column names
setnames(Xdisplay, c("Members", "Value (kN)", "Tension or Compression"))
pander(Xdisplay)| Members | Value (kN) | Tension or Compression |
|---|---|---|
| FAG | 8 | C |
| FAB | 6.928 | T |
| FGB | 3 | C |
| FGF | 5 | C |
| FBF | 1.732 | T |
| FBC | 3.465 | T |
Use GNU Octave to obtain the answer for Example Problem 3.2 (Hibbeler)
% Solution method (Bech)
% Coefficient matrix
A = [0, 0, 0, 0, 0, 6; -1, 0, 0, 0, 0, 1; 0, -1, 0, 0, 0, 0; 0, -8, 4, 0, 0, 0; -1, 0, 0, 1, 0, 0; 0, -1, 1, 0, -1, 0];
% Right-hand side, external forces
B = [300 * 2 + 300 * 6; 0; -300 - 300; -300 * 6 - 300 * 2; 0; -300 - 300];
% Solve for the unknown forces
X = A \ B; % Use GNU Octave's built-in slash operator for solving the system of equations
Ax = X(1);
Ay = X(2);
F_BD = X(3);
Cx = X(4);
Cy = X(5);
E = X(6);
disp('Solution:')
fprintf('\n');
disp('Ax = ')
fprintf('%9.0f\t', Ax);
fprintf('\n');
disp('Ay = ')
fprintf('%9.0f\t', Ay);
fprintf('\n');
disp('F_BD = ')
fprintf('%9.0f\t', F_BD);
fprintf('\n');
disp('Cx = ')
fprintf('%9.0f\t', Cx);
fprintf('\n');
disp('Cy = ')
fprintf('%9.0f\t', Cy);
fprintf('\n');
disp('E = ')
fprintf('%9.0f\t', E);Solution:
Ax = 400
Ay = 600
F_BD = 600
Cx = 400
Cy = 600
E = 400
Problem 7.17 (Olia 206-207) solved using R
# Solution method (Bech)
# Solves the system of equations AX = B, where X is a column vector of unknown
# forces, of the form X [Ax Ay F_BD Cx Cy E]'
# Coefficient matrix
A <- matrix(c(0, 0, 0, 0, 0, 6, -1, 0, 0, 0, 0, 1, 0, -1, 0, 0, 0, 0, 0, -8, 4, 0,
0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 1, 0, -1, 0), nrow = 6, byrow = TRUE)
pander(A)| 0 | 0 | 0 | 0 | 0 | 6 |
| -1 | 0 | 0 | 0 | 0 | 1 |
| 0 | -1 | 0 | 0 | 0 | 0 |
| 0 | -8 | 4 | 0 | 0 | 0 |
| -1 | 0 | 0 | 1 | 0 | 0 |
| 0 | -1 | 1 | 0 | -1 | 0 |
# Right-hand side, external forces
B <- mat("300 * 2 + 300 * 6; 0; -300 - 300; -300 * 6 - 300 * 2; 0; -300 - 300", eval = TRUE)
pander(B)| 2400 |
| 0 |
| -600 |
| -2400 |
| 0 |
| -600 |
# Solve for the unknown forces
X <- solve(A, B)
# Use base R's built-in solve function for solving the system of equations
pander(X)| 400 |
| 600 |
| 600 |
| 400 |
| 600 |
| 400 |
# or
X <- mldivide(A, B)
# Use pracma's mldivide (which is the same as the GNU Octave built-in \
# operator for solving the system of equations)
pander(X)| 400 |
| 600 |
| 600 |
| 400 |
| 600 |
| 400 |
# Display the results in a data.table
Xdisplay <- data.table(Forces = c("Ax", "Ay", "F_BD", "Cx", "Cy", "E"), Value = X)
# set the column names
setnames(Xdisplay, c("Members", "Value (N)"))
pander(Xdisplay)| Members | Value (N) |
|---|---|
| Ax | 400 |
| Ay | 600 |
| F_BD | 600 |
| Cx | 400 |
| Cy | 600 |
| E | 400 |
Use GNU Octave to obtain the answer for Olia 206-207
% Solution method (Bech)
% Coefficient matrix
A = [0, 0, 0, 0, 0, 6; -1, 0, 0, 0, 0, 1; 0, -1, 0, 0, 0, 0; 0, -8, 4, 0, 0, 0; -1, 0, 0, 1, 0, 0; 0, -1, 1, 0, -1, 0];
% Right-hand side, external forces
B = [300 * 2 + 300 * 6; 0; -300 - 300; -300 * 6 - 300 * 2; 0; -300 - 300];
% Solve for the unknown forces
X = A \ B; % Use GNU Octave's built-in slash operator for solving the system of equations
Ax = X(1);
Ay = X(2);
F_BD = X(3);
Cx = X(4);
Cy = X(5);
E = X(6);
disp('Solution:')
fprintf('\n');
disp('Ax = ')
fprintf('%9.0f\t', Ax);
fprintf('\n');
disp('Ay = ')
fprintf('%9.0f\t', Ay);
fprintf('\n');
disp('F_BD = ')
fprintf('%9.0f\t', F_BD);
fprintf('\n');
disp('Cx = ')
fprintf('%9.0f\t', Cx);
fprintf('\n');
disp('Cy = ')
fprintf('%9.0f\t', Cy);
fprintf('\n');
disp('E = ')
fprintf('%9.0f\t', E);Solution:
Ax = 400
Ay = 600
F_BD = 600
Cx = 400
Cy = 600
E = 400
Example 4.1 (Prakash) solved using R
# Solution method (Bech)
# Solves the system of equations AX = B, where X is a column vector of unknown
# forces, of the form X [Ax Ay Bx]'
# gravitational acceleration due to gravity (Wikimedia)
g <- 9.80665 # m / s^2
# Coefficient matrix
A <- matrix(c(1, 0, 1, 0, 1, 0, 0, 0, 1.5), nrow = 3, byrow = TRUE)
pander(A)| 1 | 0 | 1 |
| 0 | 1 | 0 |
| 0 | 0 | 1.5 |
# Right-hand side, external forces
B <- mat("0; 1000 * g + 2400 * g; 1000 * g * 2 + 2400 * g * 6", eval = TRUE) # Nm
pander(B)| 0 |
| 33343 |
| 160829 |
# convert from Nm to kN-m
B <- B/1000 # kN-m
pander(B)| 0 |
| 33.34 |
| 160.8 |
# Solve for the unknown forces
X <- solve(A, B)
# Use base R's built-in solve function for solving the system of equations
pander(X)| -107.2 |
| 33.34 |
| 107.2 |
# or
X <- mldivide(A, B)
# Use pracma's mldivide (which is the same as the GNU Octave built-in \
# operator for solving the system of equations)
pander(X)| -107.2 |
| 33.34 |
| 107.2 |
# Display the results in a data.table
Xdisplay <- data.table(Forces = c("Ax", "Ay", "Bx"), Value = X)
# set the column names
setnames(Xdisplay, c("Members", "Value (kN)"))
pander(Xdisplay)| Members | Value (kN) |
|---|---|
| Ax | -107.2 |
| Ay | 33.34 |
| Bx | 107.2 |
Use GNU Octave to obtain the answer for Example 4.1 (Prakash)
% Solution method (Bech)
% gravitational acceleration due to gravity (Wikimedia)
g = 9.80665; # m / s^2
% Coefficient matrix
A = [1, 0, 1; 0, 1, 0; 0, 0, 1.5];
% Right-hand side, external forces
B = [0; 1000 * g + 2400 * g; 1000 * g * 2 + 2400 * g * 6];
% convert from Nm to kN-m
B = B / 1000;
% Solve for the unknown forces
X = A \ B; % Use GNU Octave's built-in slash operator for solving the system of equations
Ax = X(1);
Ay = X(2);
Bx = X(3);
disp('Solution:')
fprintf('\n');
disp('Ax = ')
fprintf('%9.2f\t', Ax);
fprintf('\n');
disp('Ay = ')
fprintf('%9.2f\t', Ay);
fprintf('\n');
disp('Bx = ')
fprintf('%9.2f\t', Bx);
fprintf('\n');Solution:
Ax = -107.22
Ay = 33.34
Bx = 107.22
Works Cited
Arun Prakash, CEE 297: Basic Mechanics-I: Chapter 4: Equilibrium of Rigid Bodies Statics, Purdue University School of Civil Engineering, https://engineering.purdue.edu/~aprakas/CE297/CE297-Ch4.pdf.
Steven C. Chapra, Applied Numerical Methods with MATLAB for Engineers and Scientists, 2nd Edition, Boston, Massachusetts: McGraw-Hill, 2008, pages 49.
Edited by Michael Møller Bech, Morten Lykkegaard Christensen, Lars Diekhöner, Christian Frier, Olav Geil, Erik Lund, Peter Nielsen, Thomas Garm Pedersen, Bo Rosbjerg, Applications of Mathematics in Engineering and Science, School of Engineering and Science, Aalborg University, 2012, https://first.math.aau.dk/dan/ressources/applications/?file=applications-ses-eng.pdf.
Masoud Olia, Ph.D., P.E. and Contributing Authors, Barron’s Fundamentals of Engineering Exam, Hauppauge, New York: Barron’s Educational Series, Inc., 2015, pages 206-207.
Russell C. Hibbeler, Structural Analysis, Hauppauge, New York: Pearson, 2012, page 96.
Wikimedia Foundation, Inc. Wikipedia, 5 May 2016, Gravitational acceleration, https://en.wikipedia.org/wiki/Gravitational_acceleration
EcoC2S Links
EcoC2S Home
About EcoC2S
EcoC2S Services
1 Stop Shop
Products
EcoC2S Media
EcoC2S Resources
R Trainings and Resources provided by EcoC2S (Irucka Embry, EIT)
Donations accepted with Liberapay
If you would like to contribute to the continued development of Irucka Embry’s R packages and/or Irucka Embry’s R Examples, please feel free to donate via the link below:
Please feel free to review Irucka Embry (iaembry)’s profile on Liberapay.
Copyright and License
All R and GNU Octave code written by Irucka Embry is distributed under the GPL-3 (or later) license, see the GNU General Public License (GPL) page.
All written content originally created by Irucka Embry is copyrighted under the Creative Commons Attribution-ShareAlike 4.0 International License. All other written content retains the copyright of the original author(s).
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.