# Category Archives: Mathematics

the study of quantity, shape, space, and change in a rigorously logical and abstract framework

## Systems of Delay Differential Equations and the Lambert W Function Continued…

Now we reach the end of this particular journey down the road of DDE systems. As it turns out the problem of determining the stability of a DDE system with both a time dependent term and a delay term is very difficult. Once again consider the DDE system

$\dot{\bf x}\left(t\right) = {\bf Ax}\left(t\right) + {\bf A}_d {\bf x}\left(t-\tau\right)$

The method that we will use to solve this system is to use the Laplace transform with zero initial conditions.

$s{\bf X}\left(s\right) = {\bf AX}\left(s\right) + {\bf A}_d {\bf X}\left(s\right)e^{-s\tau}$

Once again, we divide out ${\bf X}\left(s\right)$, rearrange, and are left with

$\left(s{\bf I} - {\bf A}\right)e^{s\tau} = {\bf A}_d$

In order to get this into the form necessary to use the Lambert W function, we multiply both sides by $\tau e^{-{\bf A}\tau}$.

$\left(s{\bf I} - {\bf A}\right)e^{s\tau}\tau e^{-{\bf A}\tau} = {\bf A}_d \tau e^{-{\bf A}\tau}$

$\left(s{\bf I} - {\bf A}\right)\tau e^{\left(s{\bf I} - {\bf A}\right)\tau} = {\bf A}_d \tau e^{-{\bf A}\tau}$

With our current knowledge of the Lambert W function, we have this problem pegged, right? We can use the Lambert W function to solve this equation easily.

$s{\bf I} = \dfrac{1}{\tau}W_{k}\left({\bf A}_d \tau e^{-{\bf A}\tau}\right) + {\bf A}$

This means the eigenvalues of the system can be found as

$\mbox{det}\left(s{\bf I} - \left[\dfrac{1}{\tau}W_{k}\left({\bf A}_d \tau e^{-{\bf A}\tau}\right) + {\bf A}\right]\right) = 0$

However, there is a problem with this conclusion. Remember our discussion on how matrices must commute (i.e. ${\bf XY} = {\bf YX}$) in order for them to satisfy the original exponential equation leading to the Lambert W function? Look again at this equation.

$\left(s{\bf I} - {\bf A}\right)\tau e^{\left(s{\bf I} - {\bf A}\right)\tau} = {\bf A}_d \tau e^{-{\bf A}\tau}$

Well it turns out that, in general, the matrices $\left(s{\bf I} - {\bf A}\right)\tau$ and ${\bf A}_d \tau e^{-{\bf A}\tau}$ do not commute. In other words, in general,

$\left[\left(s{\bf I} - {\bf A}\right)\tau\right]\left[{\bf A}_d \tau e^{-{\bf A}\tau}\right] \neq \left[{\bf A}_d \tau e^{-{\bf A}\tau}\right]\left[\left(s{\bf I} - {\bf A}\right)\tau\right]$ (it’s true, run some numbers!!)

This is due to the fact that,for a general case, ${\bf AA}_d \neq {\bf A}_d{\bf A}$

As is turns out, there is a proposed method to deal with this problem which involves introducing an unknown matrix into the equation that forces commutivity. Thus, we will explore this in the next, and final, post.

## Systems of Delay Differential Equations and the Lambert W Function

The evaluation of the Lambert W function for matrices must be understood in order to move forward in our understanding of solving DDE systems. For the DDE system given as

$\dot{\bf x}\left(t\right) = {\bf A}_d{\bf x}\left(t-\tau\right)$

The eigenvalues (in this case we will use $s$ to denote the eigenvalues instead of the more common $\lambda$) of the system can be found as:

$\mbox{det}\left(s{\bf I} - \dfrac{1}{\tau}W_{k}\left({\bf A}_{d}\tau\right)\right) = 0$

Functions of matrices can generally be evaluated using the following relationship

$F\left({\bf A}\right) = {\bf V}F\left({\bf J}\right){\bf V}^{-1}$

where ${\bf A}$ is a matrix, ${\bf V}$ is a matrix of eigenvectors, and ${\bf J}$ is the Jordan canonical form of the matrix ${\bf A}$. The Jordan canonical form of a matrix contains all the eigenvalues of the matrix and accounts for any repeated eigenvalues. As a quick example, consider a matrix ${\bf A}$ with eigenvalues $\lambda_{1}$ and $\lambda_{2}$ and let’s assume that $\lambda_{1}$ is repeated once and $\lambda_{2}$ is repeated twice. This means that, since there are two repeated eigenvalues, there are two Jordan blocks. Then, ${\bf J}$ can be written as

${\bf J} = \left[\begin{array}{ccccc} \lambda_{1} & 1 & 0 & 0 & 0\\ 0 & \lambda_{1} & 0 & 0 & 0\\ 0 & 0 & \lambda_{2} & 1 & 0\\ 0 & 0 & 0 & \lambda_{2} & 1\\ 0 & 0 & 0 & 0 & \lambda_{2} \end{array}\right] = \left[\begin{array}{cc} {\bf J}_{1} & {\bf 0}\\ {\bf 0} & {\bf J}_{2}\end{array}\right]$

${\bf J}_{1} = \left[\begin{array}{cc} \lambda_{1} & 1\\ 0 & \lambda_{1}\end{array}\right]$

${\bf J}_{2} = \left[\begin{array}{ccc} \lambda_{2} & 1 & 0\\ 0 & \lambda_{2} & 1\\ 0 & 0 & \lambda_{2}\end{array}\right]$

Each Jordan block accounts for a single repeated eigenvalue. Now, the expression $F\left({\bf J}\right)$ is

$F\left({\bf J}\right) = \left[\begin{array}{ccccc} F\left(\lambda_{1}\right) & F'\left(\lambda_{1}\right) & 0 & 0 & 0\\ 0 & F\left(\lambda_{1}\right) & 0 & 0 & 0\\ 0 & 0 & F\left(\lambda_{2}\right) & F'\left(\lambda_{2}\right) & \frac{1}{2}F''\left(\lambda_{2}\right)\\ 0 & 0 & 0 & F\left(\lambda_{2}\right) & F'\left(\lambda_{2}\right)\\ 0 & 0 & 0 & 0 & F\left(\lambda_{2}\right) \end{array}\right]$

Thus, in general, the function of a single Jordan block (for a single repeated eigenvalue $\lambda$) can be written as

$F\left({\bf J}\right) = \left[\begin{array}{ccccc} F\left(\lambda\right) & F'\left(\lambda\right) & \frac{1}{2}F''\left(\lambda\right) & \cdots & \frac{1}{\left(n-1\right)!}F^{\left(n-1\right)}\left(\lambda\right)\\ 0 & F\left(\lambda\right) & F'\left(\lambda\right) & \cdots & \frac{1}{\left(n-2\right)!}F^{\left(n-2\right)}\left(\lambda\right)\\ \vdots & \ddots & \ddots & \ddots & \vdots\\ {} & {} & {} & F\left(\lambda\right) & F'\left(\lambda\right)\\ 0 & \cdots & \cdots & \cdots & F\left(\lambda\right) \end{array}\right]$

where $n$ is the duplicity of the eigenvalue. If there are no repeated eigenvalues, the matrix ${\bf J}$ becomes simple diagonal matrix with the eigenvalues on the diagonal (and zeros everywhere else). Therefore, since the Lambert W function is a function (same as sin or cos), the Lambert W of ${\bf A}$ can be written as

$W\left({\bf A}\right) = {\bf V}W\left({\bf J}\right){\bf V}^{-1}$

Differentiating the defining equation $x = W\left(x\right)e^{W\left(x\right)}$ for $W$ and solving for $W'$ (implicit differentiation), we obtain the following expression for the derivative of $W$

$W'\left(x\right) = \dfrac{1}{\left(1 + W\left(x\right)\right)e^{W\left(x\right)}}$
$W'\left(x\right) = \dfrac{W\left(x\right)}{x\left(1 + W\left(x\right)\right)}, \mbox{ for } x\neq 0$

Further derivatives of $W$ can be taken and a general formulation of the $n$th derivative of $W$ can be found in the paper “On the Lambert W Function” by Corless, et al. There is a Matlab function that will evaluate the Lambert W function for matrices (lambertwm), however you must search for it as it was written by a couple of doctoral students at Michigan.

One quick note before we move on to the current problems encountered with DDE systems. For the equation

${\bf Y}e^{\bf Y} = {\bf X}$,

where ${\bf X}$ and ${\bf Y}$ are matrices, the solution, ${\bf Y} = W_{k}\left({\bf X}\right)$, is valid only when ${\bf X}$ and ${\bf Y}$ commute; i.e. when ${\bf XY} = {\bf YX}$. The proof is quite simple and satisfactory. First, assume that ${\bf Y}e^{\bf Y} = {\bf X}$. It follows that

${\bf XY} = {\bf Y}e^{\bf Y}{\bf Y} = {\bf Y}^2e^{\bf Y} = {\bf YX}$

since ${\bf Y}$ will always commute with its own matrix exponential, $e^{\bf Y}$. Thus, it can be seen that the two matrices satisfy the original exponential equation if and only if they commute, and thus satisfying the Lambert W function logically follows.

The next type of system that we will encounter is a DDE system with both a delay term and a non-delay term.

${\bf x}\left(t\right) = {\bf Ax}\left(t\right) + {\bf A}_d{\bf x}\left(t-\tau\right)$

More on this to come. Stay tuned!

## Delay Differential Equations and the Lambert W Function Continued…

In order to get a more firm grasp on the Lambert W function and how it can be used to solve delay differential equations, let’s consider a system with a scalar state and a delay on the input to the system.

$\dot{x}\left(t\right) = ax\left(t\right) + a_{d}u\left(t-\tau\right)$

In this case $u\left(t\right)$ is the system input. In traditional differential equation language, this is a non-homogeneous delay differential equation (similar to the first post but this time there is also a non-delay term). So, let’s have the state be fed back into the system (state-feedback system), so that

$u\left(t\right) = x\left(t\right)$

Thus, the original equation becomes

$\dot{x}\left(t\right) = ax\left(t\right) + a_{d}x\left(t-\tau\right)$

If we assume a solution of the form $x\left(t\right) = Ce^{st}$, the equation becomes

$Cse^{st} = aCe^{st} + a_{d}Ce^{s\left(t-\tau\right)}$

Collecting terms and dividing out $C$

$e^{st}\left(s - a - a_{d}e^{-s\tau}\right) = 0$

Clearly $e^{st}$ cannot equal zero, thus the portion of the equation in the parentheses must equal zero. Thus, after re-arranging,

$s-a = a_{d}e^{-s\tau}$

$\left(s-a\right)e^{s\tau} = a_{d}$

$\left(s-a\right)\tau e^{\left(s-a\right)\tau} = a_{d}\tau e^{-a\tau}$

Thus, using the Lambert W function,

$\left(s-a\right)\tau = W_{k}\left(a_{d}\tau e^{-a\tau}\right)$

Stability of the differential equation is thus determined by using the principal branch of the Lambert W function.

$s = \dfrac{1}{\tau}W_{0}\left(a_{d}\tau e^{-a\tau}\right) + a$

The Lambert W function for a scalar can be evaluated in Matlab with the following commands.

lambertw(k,expression), or in this case,

$s = \dfrac{1}{\tau}\mbox{lambertw}\left(0,a_{d}\tau e^{-a\tau}\right) + a$

This process for a system of delay differential equations is the same as that for a scalar delay differential equation (as shown in the first Lambert W post). However, the difference lies in the evaluation of the Lambert W function for a matrix as opposed to a scalar. More on the DDE systems in the next post.

## Least Squares Linear Regression

Least-squares regression is a methodology for finding the equation of a best fit line through a set of data points. It also provides a means of describing how well the data correlates to a linear relationship. An example of data with a general linear trend is seen in the above graph. First, we will go over the derivation of the formulas from theory and then I have also appended at the end of this post Scilab code for implementation of the algorithm.

The equation of a line through a data point can be written as:

The value of any data points that are not directly on the line but are in the proximity of the line can be given by:

Where e is the vertical error between the y-value given by the line and the actual y-value of the data. The goal would be to come up with a line which minimizes this error. In least-squares regression, this is accomplished  by minimizing the sum of the squares of the errors. The sum of the squares of the errors is given by:

In order to minimize this value, the minimum finding techniques of differential calculus will be used. First take the derivative with respect to the slope.

Then with respect to the y-intercept yields:

Which can be substituted in the previous equation to solve for the slope.

The y-intercept is then:

It can be seen that these last two formulas only require knowledge about the data point coordinates and the number of points and the equation for the least squares linear regression line can be found.

Finally, below is the Scilab code implementation.

//the linear regression function takes x-values and
//y-values of data in the column vectors 'X' and 'Y' and finds
//the best fit line through the data points. It returns
//the slope and y-intercept of the line as well as the
//coefficient of determination ('r_sq').

//the function call for this should be of the form:
//'[m,b,r2]=Linear_Regression(x,y)'
function [slope, y_int, r_sq]=Linear_Regression(X, Y)
//determine the number of data points
n=size(X,'r');

//initialize each summation
sum_x=0;
sum_y=0;
sum_xy=0;
sum_x_sq=0;
sum_y_sq=0;

//calculate each sum required to find the slope, y-intercept and r_sq
for i=1:n
sum_x=sum_x+X(i);
sum_y=sum_y+Y(i);
sum_xy=sum_xy+X(i)*Y(i);
sum_x_sq=sum_x_sq+X(i)*X(i);
sum_y_sq=sum_y_sq+Y(i)*Y(i);
end

//determine the average x and y values for the
//y-intercept calculation
x_bar=sum_x/n;
y_bar=sum_y/n;

//calculate the slope, y-intercept and r_sq and return the results
slope=(n*sum_xy-sum_x*sum_y)/(n*sum_x_sq-sum_x^2);
y_int=y_bar-slope*x_bar;
r_sq=((n*sum_xy-sum_x*sum_y)/(sqrt(n*sum_x_sq-sum_x^2)*sqrt(n*sum_y_sq-sum_y^2)))^2;

//determine the appropriate axes size for plotting the data and
//linear regression line
axes_size=[min(X)-0.1*(max(X)-min(X)),min(Y)-0.1*(max(Y)-min(Y)),max(X)+0.1*(max(X)-min(X)),max(Y)+0.1*(max(Y)-min(Y))];

//plot the provided data
plot2d(X,Y,style=-4,rect=axes_size);

//plot the calculated regression line
plot2d(X,(slope*X+y_int));
endfunction

I hope this proves helpful. Let me know in the comments if you have any questions.

Related Posts:

Maximizing and Minimizing

The Bisection Method Using Scilab

Structural Finite Element Analysis Software Installation

## Delay Differential Equations and the Lambert W Function

Delay differential equations are natural to study since most systems involve a delay from the input to the output (accelerators, computers, etc.). On a more nerdy level they involve some pretty interesting math so let’s take a look.

Consider the scalar, linear, pure delay differential equation:

$\dot{x}\left(t\right) = a_{d}x\left(t-\tau\right)$

This type of equation is associated with a system where the output is equivalent to the input delayed by a small time constant $\tau$ and multiplied by a system constant $a_d$. In studying the stability of differential equations, we want to know whether they will behave in one of two ways. Either the solution to the differential equation approaches infinity as time approaches infinity (unstable) or the solution approaches some constant as time approaches infinity (stable). To do this, let’s take the Laplace transform of the DDE.

$sX(s) - x(0) = a_{d}e^{-s\tau}X(s)$

Collecting the terms and assuming the initial condition to be zero, we arrive at the following equation.

$X(s)\left(s - a_{d}e^{-s\tau}\right) = 0$

To avoid a trivial solution, only the part of the equation in the parentheses can equal zero. Therefore,

$s - a_{d}e^{-s\tau} = 0$

Let’s introduce the Lambert W function. This function satisfies the equation:

$Ye^Y = X$

The solution to this equation is:

$Y_k = W_{k}\left(X\right)$

In other words,

$W_{k}\left(X\right) e^{W_{k}\left(X\right)} = X$

where $W_{k}\left(X\right)$ is the Lambert W function of $X$ and $k$ is the branch number. This equation is called a transcendental equation since there are infinite values that satisfy this equation. The Lambert W function has infinite branches (similar to $\mbox{tan}^{-1}$ or $\mbox{sin}^{-1}$), meaning there are infinite values that will satisfy the equation above. Additionally, it can be shown that the maximum values yielded by the Lambert W function are given by the principal branch (i.e. $k=0$).

Returning to the Laplace transform equation,

$s - a_{d}e^{-s\tau}=0$

The roots of this equation determine the stability of the DDE. Therefore, we solve for the values that make this equation zero.

$s = a_{d}e^{-s\tau}$

Multiplying both sides by $\tau e^{s\tau}$ yields:

$\underbrace{s\tau}_Y e\underbrace{^{s\tau}}_Y = \underbrace{a_{d}\tau}_X$

Clearly, this equation is a candidate for using the Lambert W function. So, for this simple function, it can be seen that the roots for the DDE are given as:

$s_{k} = \dfrac{1}{\tau}W_{k}\left(a_{d}\tau\right)$

There are infinite values that satisfy this equation; however, since the maximum values for the Lambert W function are given by the principal branch, the only branch that need be evaluated is the principal branch. In other words, if the maximum value for the DDE roots is negative, then all the rest of the values are guaranteed negative. Therefore, for stability

$s_{0} = \dfrac{1}{\tau}W_{0}\left(a_{d}\tau\right) < 0$

This is an important result for understanding delay systems. However, the study of delays in systems of differential equations is much more difficult and remains an open problem in the field of dynamics and control systems. For example consider the DDE system

${\bf x}\left(t\right) = {\bf A}_{d}{\bf x}\left(t-\tau\right)$

where ${\bf x}\left(\cdot\right)$ is a vector and ${\bf A}_{d}$ is a matrix. Stay tuned for more!

## Structural Finite Element Analysis Software Installation

My engineering capstone design team and I recently completed our design and implementation of a structural finite element analysis (SFEA) tool. It allows a user to analyze any two-dimensional structure (beam, frame, truss). The user provides the geometry of the structure along with element parameters (cross-sectional area, modulus of elasticity, etc.) and the loads applied. The program will then calculate the nodal displacements, the reactions at the supports, and the forces acting on each element. The tool is built on the free, open-source platform called Scilab. My desire here is to make the code available to everyone for use or for further modification and development, and to provide instructions on its installation.

In order to run the application it will of course be necessary to have Scilab. Click here and you will be directed to the Scilab site where you can download and install the program. Once you have done this, open Scilab. This will bring up the main console. In the top menu, select ‘Applications->Scinotes’.

This will open the Scinotes editor in a new window. In this window you should paste the SFEA code from the end of this post (sorry for the length)

Once you have pasted the code into the Scinotes editor, save it as ‘SFEA.sci’ in the ‘C:\Program_Files\scilabx.x.x\modules\scinotes\macros’ folder and execute the code from the execute menu or the control button.

Now return to the main console and at the command prompt enter the function name, ‘–>SFEA’ , in all caps. When you hit enter, the program will run and open up two new windows. The left window will provide the controls for user entry while the right window will be used to plot the structure graphically.

To learn how to utilize the program, I have written a user’s guide.

If you have any questions or problems, please let me know in the comments section.

————————————————————————————————————————————————————————–

function SFEA()
//=======================================================================================
//DECLARATIONS OF GLOBAL VARIABLES AND CONSTANTS
//=======================================================================================
clearglobal();

global main_window_x;   //constants
global main_window_y;
global main_window_w;
global main_window_h;
global plot_window_x;
global plot_window_y;
global plot_window_w;
global plot_window_h;
global margin_x;
global margin_y;
global spacing_x;
global spacing_y;
global text_w;
global text_h;
global button_w;

global Truss_Switch;   //variables
global Beam_Switch;
global Node_Entry_Switch;
global Member_Entry_Switch;
global Reactions_Entry_Switch;
global Num_Nodes;
global Node_Coord;
global Axes_Size;
global Node_Type;
global Num_Members;
global Connect_Mat;
global Local_k;
global Transf_Mat;
global Displacements;
global Structure_Forces;

//==============================================================================================
//GLOBAL CONSTANT DEFINITIONS
//==============================================================================================

main_window_x      =0;
main_window_y      =0;
main_window_w      =335; //width
main_window_h      =600; //height

plot_window_x      =main_window_x+main_window_w+18;
plot_window_y      =main_window_y;
plot_window_w      =640;
plot_window_h      =main_window_h;

margin_x           =10;
margin_y           =15;

spacing_x          =10;
spacing_y          =10;

text_w             =145;
text_h             =25;
button_w           =120;

//==================================================================================================
//CREATE THE PRIMARY WINDOWS AND MENUS
//==================================================================================================
//create the figure 'Main_Window'
Main_Window=figure(1);
Main_Window.Figure_name = "Structural Finite Element Analysis";
Main_Window.Position = [main_window_x,main_window_y,main_window_w,main_window_h];
Main_Window.Background = 8; //white

//create the figure 'Plot_Window'
Plot_Window=figure(2);
Plot_Window.Figure_name = "Structure Plot"
Plot_Window.Position = [plot_window_x,plot_window_y,plot_window_w,plot_window_h];
Plot_Window.Background = 8; //white

toolbar(Main_Window.figure_id, "off");

endfunction
//========================================================================================================================
//========================================================================================================================
function Close_Program()
scf(2);
close(); //close the plot window
scf(1);
close(); //close the main window
endfunction
//========================================================================================================================
function Beam_Analysis()
global Truss_Switch;
global Beam_Switch;
global Node_Entry_Switch;
global Member_Entry_Switch;
global Reactions_Entry_Switch;
Truss_Switch = 0;
Beam_Switch = 1;
Node_Entry_Switch = 0;
Member_Entry_Switch = 0;
Reactions_Entry_Switch = 0;
Clear_Figures();
endfunction
//========================================================================================================================
function Frame_Analysis()
global Truss_Switch;
global Beam_Switch;
global Node_Entry_Switch;
global Member_Entry_Switch;
global Reactions_Entry_Switch;
Truss_Switch = 0;
Beam_Switch = 0;
Node_Entry_Switch = 0;
Member_Entry_Switch = 0;
Reactions_Entry_Switch = 0;
Clear_Figures();
endfunction
//========================================================================================================================
function Truss_Analysis()
global Truss_Switch;
global Beam_Switch;
global Node_Entry_Switch;
global Member_Entry_Switch;
global Reactions_Entry_Switch;
Truss_Switch = 1;
Beam_Switch = 0;
Node_Entry_Switch = 0;
Member_Entry_Switch = 0;
Reactions_Entry_Switch = 0;
Clear_Figures();
endfunction
//========================================================================================================================
//SUB-ROUTINES
//========================================================================================================================
function Clear_Figures()
scf(2);
clf();
scf(1);
f=gcf();
C=f.children;
for i=1:length(C)
c=C(i);
delete(c);
end
end
endfunction
//========================================================================================================================
//========================================================================================================================
global main_window_x;
global main_window_y;
global main_window_w;
global main_window_h;
global plot_window_x;
global plot_window_y;
global plot_window_w;
global plot_window_h;
global margin_x;
global margin_y;
global spacing_x;
global spacing_y;
global text_w;
global text_h;
global button_w;

//================================================================================================
scf(1); //set current figure to 'Main_Window'
//================================================================================================
//frames:
//inputs frame
Frame1=uicontrol(gcf(),"style","frame");
Frame1.Position = [(margin_x/2),435,320,185];
Frame1.HorizontalAlignment = "center";
Frame1.BackgroundColor = [1,1,1]; //white

//results frame
Frame2=uicontrol(gcf(),"style","frame");
Frame2.Position = [(margin_x/2),300,320,130];
Frame2.HorizontalAlignment = "center";
Frame2.BackgroundColor = [1,1,1]; //white
//==================================================================================================
//1st Line: number of nodes and their coordinates
//creates text in the main window
Node_Text=uicontrol(gcf(),"style","text");
Node_Text.Position = [margin_x,main_window_h-margin_y,text_w,text_h];
Node_Text.BackgroundColor = [1,1,1]; //white
Node_Text.String = "Number of Nodes";

//creates edit box for entering data in the main window
Num_Node_Box=uicontrol(gcf(),"style","edit");
Num_Node_Box.Position = [margin_x+text_w+spacing_x,main_window_h-margin_y,text_h,text_h];
Num_Node_Box.String = "0";
Num_Node_Box.Tag = "NNB"

//create a pushbutton in the main window
Node_Button=uicontrol(gcf(),"style","pushbutton");
Node_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y,button_w,text_h];
Node_Button.String = "Node Coordinates";
Node_Button.Callback = "Get_Node_Coord()";
//====================================================================================================
//2nd Line: node reactions
//creates text in the main window
Reactions_Text=uicontrol(gcf(),"style","text");
Reactions_Text.Position = [margin_x,main_window_h-margin_y-text_h-spacing_y,text_w,text_h];
Reactions_Text.BackgroundColor = [1,1,1]; //white
Reactions_Text.String = "Node Reactions";

//create a pushbutton in the main window
Reactions_Button=uicontrol(gcf(),"style","pushbutton");
Reactions_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y-text_h-spacing_y,button_w,text_h];
Reactions_Button.String = "Select Reactions";
Reactions_Button.Callback = "Get_Node_Reactions()";
//====================================================================================================
//3rd Line: number of members and their connectivity
//creates text in the main window
Member_Text=uicontrol(gcf(),"style","text");
Member_Text.Position = [margin_x,main_window_h-margin_y-2*text_h-2*spacing_y,text_w,text_h];
Member_Text.BackgroundColor = [1,1,1]; //white
Member_Text.String = "Number of Members";

//creates edit box for entering data in the main window
Num_Member_Box=uicontrol(gcf(),"style","edit");
Num_Member_Box.Position = [margin_x+text_w+spacing_x,main_window_h-margin_y-2*text_h-2*spacing_y,text_h,text_h];
Num_Member_Box.String = "0";
Num_Member_Box.Tag = "NMB"

//create a pushbutton in the main window
Member_Button=uicontrol(gcf(),"style","pushbutton");
Member_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y-2*text_h-2*spacing_y,button_w,text_h];
Member_Button.String = "Member Data";
Member_Button.Callback = "Get_Member_Data()";
//====================================================================================================================
//creates text in the main window

//creates edit box for entering data in the main window

//create a pushbutton in the main window
//=========================================================================================================================
//creates text in the main window
Combos_Text=uicontrol(gcf(),"style","text");
Combos_Text.Position = [margin_x,main_window_h-margin_y-4*text_h-4*spacing_y,text_w,text_h];
Combos_Text.BackgroundColor = [1,1,1]; //white
Combos_Text.String = "Number of Load Combinations"

//creates edit box for entering data in the main window
Num_Combos_Box=uicontrol(gcf(),"style","edit");
Num_Combos_Box.Position = [margin_x+text_w+spacing_x,main_window_h-margin_y-4*text_h-4*spacing_y,text_h,text_h];
Num_Combos_Box.String = "0";
Num_Combos_Box.Tag = "NLCB"

//create a pushbutton in the main window
Combos_Button=uicontrol(gcf(),"style","pushbutton");
Combos_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y-4*text_h-4*spacing_y,button_w,text_h];
//=========================================================================================================================
//6th Line: find displacements
//create a pushbutton in the main window
Displacements_Button=uicontrol(gcf(),"style","pushbutton");
Displacements_Button.Position = [margin_x+10*spacing_x,main_window_h-margin_y-5.5*text_h-5.5*spacing_y,button_w,text_h];
Displacements_Button.String = "Find Displacements";
Displacements_Button.Callback = "Find_Displacements()";
//=========================================================================================================================
//7th Line: find reactions
//create a pushbutton in the main window
Reactions_Button=uicontrol(gcf(),"style","pushbutton");
Reactions_Button.Position = [margin_x+10*spacing_x,main_window_h-margin_y-6.5*text_h-6.5*spacing_y,button_w,text_h];
Reactions_Button.String = "Find Reactions";
Reactions_Button.Callback = "Find_Reactions()";
//=========================================================================================================================
//8th Line: find member forces
//create a pushbutton in the main window
Reactions_Button=uicontrol(gcf(),"style","pushbutton");
Reactions_Button.Position = [margin_x+10*spacing_x,main_window_h-margin_y-7.5*text_h-7.5*spacing_y,button_w,text_h];
Reactions_Button.String = "Find Member Forces";
Reactions_Button.Callback = "Find_Forces()";
endfunction
//=============================================================================================================================
//Data Entry and Callback Functions
//=============================================================================================================================
function Get_Node_Coord()
//Extracts the number of nodes submitted by the user, obtains the node coordinates, and plots the nodes in the plot window
global Num_Nodes;
global Node_Entry_Switch;
global Node_Coord;
global Axes_Size;

//Extract Num_Nodes from Edit Box in Main_Window
Num_Nodes = evstr(get(findobj("Tag","NNB"), "String"));

//Used to initialize the node coordinate matrix to zero during the first iteration
//but maintain the inputted node coordinate entries during future callbacks
if Node_Entry_Switch == 0 then
Node_Coord=zeros(Num_Nodes,2);
Node_Entry_Switch = 1;
end

//Let the user select the Node Coordinate data entry method

//execute the users selection
select Node_Entry_Method
case 1 then //======================================================================================================
//Enter Nodes Singularly:

//initializes dynamically sized matrix for coordinate labels
Row_Labels=zeros(Num_Nodes,1);
for i=1:Num_Nodes
Row_Labels(i,1) = i;
end

//displays a matrix dialog box for user inputs and saves the coordinates to the matrix 'Node_Coord'
Node_Coord=evstr(x_mdialog('Enter the Coordinates of the Nodes',string(Row_Labels),['X','Y'],string(Node_Coord)));

case 2 then //===========================================================================================================
//Enter Nodes via Sets:

//Obtain the number of sets to be entered from the user
Num_Sets = evstr(x_dialog('How many sets do you want to enter?','0'));

//initializes dynamically sized matrix for set labels
Row_Labels=zeros(Num_Sets);
for i=1:Num_Sets
Row_Labels(i) = i;
end

Set_Data = zeros(Num_Sets,7); //initialize the Set_Data matrix

//displays a matrix dialog box for user inputs of the set data
Set_Data=evstr(x_mdialog('Enter the Set Data',string(Row_Labels),['Node Number Start','Node Number Increment','X Start','Y Start','X Increment','Y Increment','Number of Nodes in Set'],string(Set_Data)));

//take the data entered into the 'Set_Data' matrix and determine the coordinate matrix it produces
for i=1:Num_Sets //for each set
for j=0:(Set_Data(i,7)-1) //for each node in the set
Node_Coord(Set_Data(i,1)+j*Set_Data(i,2),1)=Set_Data(i,3)+j*Set_Data(i,5);
Node_Coord(Set_Data(i,1)+j*Set_Data(i,2),2)=Set_Data(i,4)+j*Set_Data(i,6);
end
end

case 3 then //================================================================================================================
//Enter Nodes from a File:

File_Directory = uigetfile(["*.txt"],"SCI/modules/scinotes/macros/","Choose a file");
fi=mopen(File_Directory,'r');
Node_Coord=evstr(mgetl(fi,Num_Nodes));
mclose(fi);
end //========================================================================================================================

//plots the nodes in the plot window
Axes_Size = zeros(4);
X=Node_Coord;
Y=Node_Coord;
X(:,2)=[];
Y(:,1)=[];
scf(2);
clf(2);

//determine the axes size for plotting
//initial estimates
Axes_Size(1)=min(Node_Coord(:,1))-((max(Node_Coord(:,1))-min(Node_Coord(:,1)))/10);
Axes_Size(2)=min(Node_Coord(:,2))-((max(Node_Coord(:,2))-min(Node_Coord(:,2)))/10);
Axes_Size(3)=max(Node_Coord(:,1))+((max(Node_Coord(:,1))-min(Node_Coord(:,1)))/10);
Axes_Size(4)=max(Node_Coord(:,2))+((max(Node_Coord(:,2))-min(Node_Coord(:,2)))/10);

//check for good proportionality
width=Axes_Size(3)-Axes_Size(1);
height=Axes_Size(4)-Axes_Size(2);
if height==0 then
aspect_ratio=3;
else
aspect_ratio=width/height;
end

if aspect_ratio<0.5 then
new_width=height/2;
Axes_Size(1)=Axes_Size(1)-((new_width-width)/2);
Axes_Size(3)=Axes_Size(3)+((new_width-width)/2);
elseif aspect_ratio>2 then
new_height=width/2;
Axes_Size(2)=Axes_Size(2)-((new_height-height)/2);
Axes_Size(4)=Axes_Size(4)+((new_height-height)/2);
end

//plot points
plot2d(X,Y,style=-4,rect=Axes_Size);

//print out the node coordinate matrix
print(6,Node_Coord);
endfunction
//===========================================================================================================
//===========================================================================================================
function Get_Node_Reactions()
//Presents a toggle dialog box for the user to select a type for each node and then
//re-plots the nodes in the plot window with representative symbols
global Reactions_Entry_Switch;
global Num_Nodes;
global Node_Type;
global Node_Coord;
global Axes_Size;

//Used to initialize the node type matrix to 1s during the first iteration
//but maintain the inputted node type entries during future iterations
if Reactions_Entry_Switch == 0 then
for j=1:Num_Nodes
Node_Type(j) = 1;
end
Reactions_Entry_Switch = 1;

end

//generates a toggle selection dialog box for node type selection
//current options are: Free, X-Roller, Y-Roller, Pin, Fixed, X-Spring, or Y-String
//Free is the default value
L=list();
for i=1:Num_Nodes
s='Node'+string(i);
l  = list(s,Node_Type(i),['Free','X-Roller','Y-Roller','Pin','Fixed']);
L(i)=l;
end
Node_Type = x_choices('Select Node Type',L);

//re-plot the nodes as their respective node types
X=Node_Coord;
Y=Node_Coord;
X(:,2)=[];
Y(:,1)=[];
scf(2);
clf(2);
for i=1:Num_Nodes
select Node_Type(i)
case 1 then
s=-4; //Free
case 2 then
s=-3; //X-Roller
case 3 then
s=-9; //Y-Roller
case 4 then
s=-6; //Pin
case 5 then
s=-11; //Fixed
end
plot2d(X(i),Y(i),style=s,rect=Axes_Size);
end

//print out the node type matrix
print(6,Node_Type);

endfunction
//==============================================================================================================
//==============================================================================================================
function Get_Member_Data()
//Extracts the number of members from the user input, presents an x_mdialog box for the
//user to input the starting and ending nodes of each member as well as the cross-secional
//area, modulus of elasticity and moment of inertia, and plots the members
//connecting the nodes in the plot window
global Num_Members;
global Member_Entry_Switch;
global Connect_Mat;
global Node_Coord;
global Axes_Size;

//Extract Num_Members from Edit Box in Main_Window
Num_Members = evstr(get(findobj("Tag","NMB"), "String"));

//Used to initialize the member data matrix to zero during the first iteration
//but maintain the inputted member data entries during future iterations
if Member_Entry_Switch == 0 then
Connect_Mat=zeros(Num_Members,5);
Member_Entry_Switch = 1;
end

//Let the user select the Member data entry method

//execute the users selection
select Member_Entry_Method
case 1 then //======================================================================================================
//enter Member Data Singularly

//initializes dynamically sized matrix for member labels
Row_Labels=zeros(Num_Members);
for i=1:Num_Members
Row_Labels(i) = i;
end

//displays a matrix dialog box for user inputs and saves the member connections to the matrix 'Connect_Mat'
Connect_Mat=evstr(x_mdialog('Enter the Members Information',string(Row_Labels),['Starting Node','Ending Node','Cross-Sectional Area (A)','Modulus of Elasticity (E)','Moment of Inertia (I)'],string(Connect_Mat)));
case 2 then //=========================================================================================================
//enter Member Data by sets

//Obtain the number of sets to be entered from the user
Num_Sets = evstr(x_dialog('How many sets do you want to enter?','0'));

//initializes dynamically sized matrix for set labels
Row_Labels=zeros(Num_Sets);
for i=1:Num_Sets
Row_Labels(i) = i;
end

Set_Data = zeros(Num_Sets,9); //initialize the Set_Data matrix

//displays a matrix dialog box for user inputs of the set data
Set_Data=evstr(x_mdialog('Enter the Set Data',string(Row_Labels),['Member Number Start','Member Number Increment','First Start Node','First End Node','Node Increment','Number of Members in Set','A','E','I'],string(Set_Data)));

//take the data entered into the 'Set_Data' matrix and determine the connectivity matrix it produces
for i=1:Num_Sets //for each set
for j=0:(Set_Data(i,6)-1) //for each member in the set
Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),1)=Set_Data(i,3)+j*Set_Data(i,5); //extracts member starting nodes
Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),2)=Set_Data(i,4)+j*Set_Data(i,5); //extracts member ending nodes
Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),3)=Set_Data(i,7); //extracts member cross-sectional area
Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),4)=Set_Data(i,8); //extracts member modulus of elasticity
Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),5)=Set_Data(i,9); //extracts member moment of inertia
end
end
case 3 then //=========================================================================================================
//enter Member Data from a file

File_Directory = uigetfile(["*.txt"],"SCI/modules/scinotes/macros/","Choose a file");
fi=mopen(File_Directory,'r');
Connect_Mat=evstr(mgetl(fi,Num_Members));
mclose(fi);
end //=================================================================================================================

//Draws the members between the nodes with a line
for i=1:Num_Members
for j=1:2
X(j)=Node_Coord(Connect_Mat(i,j),1);
Y(j)=Node_Coord(Connect_Mat(i,j),2);
end
plot2d(X,Y,style=1,rect=Axes_Size);
end

//print out the connectivity matrix
print(6,Connect_Mat);

endfunction
//============================================================================================================================
//============================================================================================================================
//obtains the number of load types (i.e. dead, live, seismic, etc.) and presents the user
//with an x_mdialog box to enter the horizontal and vertical loads as well as applied
//moments at each node for each load type
global Num_Nodes;

//Extract Num_Load_Types from Edit Box in Main_Window

//initializes dynamically sized matrix for loads entry and associated
//labels
Row_Labels=zeros(Num_Nodes);
for i=1:Num_Nodes
Row_Labels(i) = i;
end

//displays a matrix dialog box for user inputs

for i=1:Num_Nodes //for each node
for j=1:3 //for each applied load
end
end
end

endfunction
//=====================================================================================================================
//=====================================================================================================================
//obtains the number of load combinations from the user and provides an x_mdialog
//combination
global Num_Nodes;
global Axes_Size;
global Node_Coord;
global Num_Members;
global Connect_Mat;

//Extract Num_Load_Combos from Edit Box in Main_Window

//initializes dynamically sized matrix for members and associated
//labels

end
Column_Labels(i) = 'Combination '+string(i);
end

//displays a matrix dialog box for user inputs and saves the member connections to the matrix 'Connect_Mat'

for i=1:(3*Num_Nodes)
else
end
end

//define a unit length for graphing the forces
L=zeros(Num_Members,1);
for i=1:Num_Members
Start_Node=Connect_Mat(i,1); //starting node number
End_Node=Connect_Mat(i,2); //ending node number
X1=Node_Coord(Start_Node,1); //starting node coordinates
Y1=Node_Coord(Start_Node,2);
X2=Node_Coord(End_Node,1); //ending node coordinates
Y2=Node_Coord(End_Node,2);
L(i)=(((Y2-Y1)^2)+((X2-X1)^2))^(1/2); //length of the member
end
unit_size=min(L)/8;

//graph the ultimate applied forces
for i=1:Num_Nodes
for j=1:3
select j
case 1 then //horizontal forces
X(1)=Node_Coord(i,1);
Y(1)=Node_Coord(i,2);
Y(2)=Node_Coord(i,2);
case 2 then //vertical forces
X(1)=Node_Coord(i,1);
Y(1)=Node_Coord(i,2);
X(2)=Node_Coord(i,1);
case 3 then //moments
x = Node_Coord(i,1)-unit_size/2;
y = Node_Coord(i,2)+unit_size/2;
w = unit_size;
h = unit_size;
a1 = -80*64;
a2 = 160*64;
a1 = 100*64;
a2 = 160*64;
end
//plot moment arcs
xarcs([x;y;w;h;a1;a2],[5]); //5 - red
end
//plot force lines
plot2d(X,Y,style=5,rect=Axes_Size); //5 - red
end
end
end

endfunction
//=====================================================================================================================
//Find Results Subroutines
//=====================================================================================================================
function Find_Displacements()
//this function calculates the local and global stiffness matrices as well as
//the transformation matrices for converting between the two. Finally, it finds
//the node displacements and plots the deformed structure
global Num_Nodes;
global Num_Members;
global Connect_Mat;
global Node_Coord;
global Local_k;
global Transf_Mat;
global Truss_Switch;
global Beam_Switch;
global Node_Type;
global Displacements;
global Axes_Size;
global Structure_Forces;

//global structure stiffness matrix
K=zeros(3*Num_Nodes,3*Num_Nodes);

for i=1:Num_Members
//obtain the parameters for the member stiffness matrices
Start_Node=Connect_Mat(i,1); //starting node number
End_Node=Connect_Mat(i,2); //ending node number
AE=Connect_Mat(i,3)*Connect_Mat(i,4); //the product of A and E
EI=Connect_Mat(i,4)*Connect_Mat(i,5); //the product of E and I
X1=Node_Coord(Start_Node,1); //starting node coordinates
Y1=Node_Coord(Start_Node,2);
X2=Node_Coord(End_Node,1); //ending node coordinates
Y2=Node_Coord(End_Node,2);
L=(((Y2-Y1)^2)+((X2-X1)^2))^(1/2); //length of the member

//find the member directional sines and cosines
if X2==X1 then //the member is vertical
if Y1<Y2 then //the member points up
phi = (%pi/2); //angle of the member with the horizontal
elseif Y1>Y2 then //the member points down
phi = ((3*%pi)/2);
end
elseif Y2==Y1 then //the member is horizontal
if X1<X2 then //the member points to the right
phi = 0;
elseif X1>X2 then //the member points to the left
phi = %pi;
end
else //the member is at some intermediate angle
phi=atan((Y2-Y1)/(X2-X1));
end
C=cos(phi);
S=sin(phi);
//===================================================================================================
//STIFFNESS MATRIX
//===================================================================================================
//compute the stiffness matrix for each member in terms of the local
//coordinate axis
Local_k(1:3,1:6,i)=[(AE/L),0,0,(-AE/L),0,0;0,(12*EI/L^3),(6*EI/L^2),0,(-12*EI/L^3),(6*EI/L^2);0,(6*EI/L^2),(4*EI/L),0,(-6*EI/L^2),(2*EI/L);];
Local_k(4:6,1:6,i)=[(-AE/L),0,0,(AE/L),0,0;0,(-12*EI/L^3),(-6*EI/L^2),0,(12*EI/L^3),(-6*EI/L^2);0,(6*EI/L^2),(2*EI/L),0,(-6*EI/L^2),(4*EI/L);];

//compute the transformation matrix for converting forces and displacements
//from the local coordinate axis to the global structure coordinate system
Transf_Mat(1:3,1:6,i)=[C,S,0,0,0,0;-S,C,0,0,0,0;0,0,1,0,0,0;];
Transf_Mat(4:6,1:6,i)=[0,0,0,C,S,0;0,0,0,-S,C,0;0,0,0,0,0,1;];

//stiffness matrix for each member in terms of global coordinate system
Global_k(:,:,i)=Transf_Mat(:,:,i)'*Local_k(:,:,i)*Transf_Mat(:,:,i); //k*=T'kT

//add the global stiffness matrix for each member into the global structure stiffness matrix,K
K((3*Start_Node-2):(3*Start_Node),(3*Start_Node-2):(3*Start_Node))= K((3*Start_Node-2):(3*Start_Node),(3*Start_Node-2):(3*Start_Node))+Global_k(1:3,1:3,i); //adds submatrix global_k11 into K
K((3*Start_Node-2):(3*Start_Node),(3*End_Node-2):(3*End_Node))= K((3*Start_Node-2):(3*Start_Node),(3*End_Node-2):(3*End_Node))+Global_k(1:3,4:6,i); //adds submatrix global_k12 into K
K((3*End_Node-2):(3*End_Node),(3*Start_Node-2):(3*Start_Node))= K((3*End_Node-2):(3*End_Node),(3*Start_Node-2):(3*Start_Node))+Global_k(4:6,1:3,i); //adds submatrix global_k21 into K
K((3*End_Node-2):(3*End_Node),(3*End_Node-2):(3*End_Node))= K((3*End_Node-2):(3*End_Node),(3*End_Node-2):(3*End_Node))+Global_k(4:6,4:6,i); //adds submatrix global_k22 into K
//======================================================================================================
end

//create an index matrix to keep up with what index of rows are left after deletion due to
//boundary conditions
Index=zeros(3*Num_Nodes,1);
for i=1:(3*Num_Nodes)
Index(i,1) = i;
end

Eliminate = [];
Counter = 1;
//if the structure is a truss there are no moments
if Truss_Switch == 1 then
for i = 1:Num_Nodes
Eliminate(Counter)=3*i;
Counter = Counter+1;
end
end

//if the structure is a beam there are no axial forces
if Beam_Switch == 1 then
for i = 1:Num_Nodes
Eliminate(Counter)=3*i-2;
Counter = Counter+1;
end
end

//extract boundary conditions from the Node_Type selections
for i=1:Num_Nodes
select Node_Type(i)
case 1 then
u=1;//Free
v=1;
theta=1;
case 2 then
u=0; //X-Roller
v=1;
theta=1;
case 3 then
u=1; //Y-Roller
v=0;
theta=1;
case 4 then
u=0; //Pin
v=0;
theta=1;
case 5 then
u=0; //Fixed
v=0;
theta=0;
end

//Determine where the boundary conditions give zero displacement
//these rows/cols will be deleted to generate the final
//invertible matrix
if u==0 then
Eliminate(Counter)=3*i-2;
Counter=Counter+1;
end
if v==0 then
Eliminate(Counter)=3*i-1;
Counter=Counter+1;
end
if theta==0 then
Eliminate(Counter)=3*i;
Counter=Counter+1;
end
end

//before removing rows and columns, preserve the full vectors and matrices
k = K;

//delete rows and columns corresponding to given parameters and boundary conditions
k(Eliminate,:)=[];
k(:,Eliminate)=[];
Index(Eliminate,:)=[];

//calculate the unknown displacements

//Develop the full displacement matrix
Displacements = zeros(3*Num_Nodes,1);
for i=1:size(Index,"r") //for each element in the index
//copy the rows from the calculated Unknown_Displacements matrix into the correct rows
//of the full Displacements matrix
Displacements(Index(i,1),:) = Unknown_Displacements(i,:);
end

//plot the deflected shape
Def_Node_Coord = zeros(Num_Nodes,2);
for i=1:Num_Nodes
Def_Node_Coord(i,1) = Node_Coord(i,1)+Displacements(3*(i-1)+1,1); //x-coordinates
Def_Node_Coord(i,2) = Node_Coord(i,2)+Displacements(3*(i-1)+2,1); //y-coordinates
end
for i=1:Num_Members
for j=1:2
X(j)=Def_Node_Coord(Connect_Mat(i,j),1);
Y(j)=Def_Node_Coord(Connect_Mat(i,j),2);
end
plot2d(X,Y,style=2,rect=Axes_Size);
end

//Calculate the global force vector for finding the reactions
Structure_Forces = K*Displacements;

//print out the global structure displacements vector
print(6,Displacements);

endfunction
//=====================================================================================================
//=====================================================================================================
function Find_Reactions()
//this function determines the reactions at the structure supports
global Num_Nodes;
global Node_Type;
global Structure_Forces;

//Extract the reactions from the Structure_Forces matrix
Reactions = zeros(3*Num_Nodes,1);
for i=1:Num_Nodes
if Node_Type(i)==1 then
Reactions((3*i-2:3*i),1)=0;
else
Reactions((3*i-2:3*i),1)=Structure_Forces((3*i-2:3*i),1);
end
end

//print out the reactions vector
print(6,Reactions);
endfunction
//=====================================================================================================
//=====================================================================================================
function Find_Forces()
//this function uses the global displacements to determine the forces on each
//member of the structure in terms of that member's local coordinate axis
global Num_Members;
global Connect_Mat;
global Displacements;
global Local_k;
global Transf_Mat;

//Adjust the forces to local coordinates
Local_Member_Forces = zeros(6,Num_Members); //the i-th column vector represents the forces on the i-th member
//in terms of the local coordinate axis
Member_Displacements=zeros(6,1); //for extracted each member's displacements in terms of global coordinates

for i=1:Num_Members
Start_Node=Connect_Mat(i,1); //starting node number
End_Node=Connect_Mat(i,2); //ending node number

//extract the member's displacements from the structure displacements matrix
Member_Displacements(1:3,1)=Displacements(3*Start_Node-2:3*Start_Node);
Member_Displacements(4:6,1)=Displacements(3*End_Node-2:3*End_Node);

Local_Member_Forces(:,i)=Local_k(:,:,i)*Transf_Mat(:,:,i)*Member_Displacements(:,1); //f=kTD

end

//print out the local member forces matrix
print(6,Local_Member_Forces);
endfunction

## Structural Finite Element Analysis User’s Guide

Now that you have downloaded and installed Scilab, the Structural Finite Element Analysis (SFEA) program, and got it up and running, how do you use it. Well that’s what we want to go through at this point. What you should see on your desktop are two blank screens. The one on the left will provide the user input controls while the one of the right will graphically depict the structure being analyzed. If this is not the case then you should read here first.

In the left screen there is a ‘File’ menu. If you click it you will see that you can select ‘New’ or ‘Close’. ‘Close’ is pretty obvious. ‘New’ will give you a sub-menu of three choices: Beam, Frame, or Truss. To make clear the difference, in this program an element of a beam supports transverse and moment forces, a frame supports axial, transverse and moment forces, and a truss element supports axial and transverse forces.

To give an example of why this might be important, consider a simple cantilever beam with a non-orthogonal loading. (In other words, the load does not meet the beam’s longitudinal axis at 90 degrees) Now although it is a ‘beam’ due to its form and positioning, since the load will produce an axial force in the structure, it is not a ‘beam’ with regard to the methodology of this program. If you were to analyze this structure, you would select ‘Frame’.

Once you make this selection, the left window will populate the controls. You will notice that there are two frames. The top frame walks the user through the step-by-step acquisition of the necessary data. The bottom frame contains the controls for obtaining the results of node displacements, support reactions, and elemental or member forces.

Entering Node Data

Now that the controls are up, the user simply starts at the top left and goes across and down just like reading text. Start by entering the number of nodes in the text box at the top and then click the ‘Node Coordinates’ button. This brings up a dialog box that asks the user if they want to enter the nodes ‘Singularly’, ‘By Sets’, or ‘From File’. ‘Singularly’ means manually enter the node coordinates one at a time. ‘By Sets’ means the user will enter some criteria about the nodes that is repetitive and the computer will determine the nodes for you. This is a helpful function when you are entering a large number of nodes that follow some consistent pattern with regard to their geometric arrangement. ‘From File’ means the user can open and import the coordinates from a user created text file that contains all the node coordinates.

Singularly

If the user selects this method, a new dialog box will open to let the user input the node coordinates.

The number of rows in the dialog box will dynamically correspond to the number of nodes the user entered would be in the structure.

By Sets

If the user selects this method, another dialog will appear to ask the user how many sets they have to enter. This is a simple integer input.

After entering the number of sets, another dialog box will be displayed to obtain the set data from the user.

The rows in the dialog box correspond to the number of sets. For each set the user must enter the node number to start with, how the node     number for each subsequent node should increment from the beginning node, the beginning x-y coordinate for the first node, the x and y distance each subsequent node will increment from the previous one, and how many nodes in the set.

For example if I was analyzing a simple cantilever beam and I wanted to divide the beam into 20 separate elements, then I would have 21 nodes. Since I would space the nodes out evenly along the length of the beam, I could enter those nodes by sets without having to enter them one at a time. If the length of the entire beam is 120 inches (10 feet), then I would have 20 elements each 6 inches long. So in entering the data by sets in the dialog above I would input: node number start = 1, node number increment = 1, x start = 0, y start = 0, x increment = 6, y increment 0, number of nodes in set = 21. The program would generate 21 nodes starting at the x-y coordinate (0,0) and at every 6 inches horizontally for the full length of the beam.

From File

If the user selects this method, a standard Windows dialog box will open where you can search in the file tree structure to find the text file containing the node coordinates data. The text file should be formatted with one node coordinate per line starting with the x-coordinate followed by a space followed by the y-coordinate.

Which ever method the user selects, once the data has been provided to the program, it will output the node coordinate matrix to the Scilab console. The user can verify the correct coordinates numerically there. Additionally, the nodes will be plotted graphically in the second SFEA window. This allows the user to verify the nodes visually.

Entering Support Reactions Data

Once the nodes have been entered, the user can set the boundary conditions for the structure by pushing the ‘Select Reactions’ button. A dialog box will then open and allow the user to select the support type at each node.

Notice the default node type is free, so the user only has to define those nodes which are supports. An ‘X-Roller’ allows rotation and movement in the y-direction but no movement in the x-direction. Likewise a ‘Y-Roller’ allows rotation and movement in the x-direction but no movement in the y-direction. A pin allows rotations but no x or y movement. Finally a fixed support allows no movement at all.

Once the user defines the node types, the node type matrix will be displayed in the Scilab console and the structure plot will be modified with different symbology for different node types as seen below. The symbols are from left to right, free, x-roller, y-roller, pin, and fixed.

Entering Member or Element Data

Once the nodes and reactions are defined the user can input the members that run between the nodes. In much the same way as the entry of the nodes, the user will enter the number of members in the textbox in the SFEA first window and then push the ‘Member Data’ button. As with the nodes, a dialog box will permit the user to input the data ‘Singularly’, ‘By Sets’, or ‘From File’.

Singularly

If the user selects this method, a dialog box will appear to allow the user to input the data for each member. The data required are the start node number, the end node number, the member’s cross-sectional area (A), the member’s modulus of elasticity (E), and the member’s moment of inertia (I).

By Sets

In the same way as with node entry, the ‘By Sets’ method of member entry allows the user to systematically enter a large number of members with repetitive characteristics. With regard to member entry, each member of a set must have the same A, E, and I.

From File

As before, this allows the user to load the member data from a text file on the computer. The format of the file should be one member per line with each line having: start node number-space-end node number-space-A-space-E-space-I

Once the member data is entered, this data is output to the Scilab console as the connectivity matrix. Also the members will be plotted in the SFEA plot window.

Where each ‘x’ is a load type and each ‘A’ is a load factor. A load factor is something that is typically used in Load and Resistance Factor Design (LRFD) as a form of scaling factor or safety factor. So for example we might have:

Where D stands for dead load and L stands for live load. In many structural applications and building codes there are sets of multiple load combinations and corresponding load factors. Each of these combinations must be calculated and the worst case scenario is then selected for design. The SFEA program provides the capacity for this type of analysis.

Once the user enters the number of load types and pushes the ‘Loads’ button,  a dialog box that allows the entry of the different loads for each type will appear.

Each row in the dialog box corresponds to each node in the structure. Consequently you can enter a horizontal, vertical, and moment load at each node. After pressing ‘OK’ another identical dialog box will appear for the next load type.

The loads matrix will be displayed in the Scilab console after entry. Each column of the matrix is a different load type. The rows follow the following pattern. Each node has 3 potential loads therefore the first 3 rows correspond to the 3 loads on node 1. The next 3 rows correspond to node 2 and so on. Therefore there should be 3N rows where N is the number of nodes. Additionally, the order of the loads in each set of 3 rows is the same as it is entered in the dialog box: horizontal, vertical, moment.

As discussed above, the number of load combinations and load factors are often listed in manuals or code books. Once the number of combinations is entered and the user pushes the ‘Load Factors’ button, the following dialog box appears.

Here the user enters the load factors that correspond to each load type for each combination.

Once this data is entered the load factors will be displayed to the Scilab console as well as the ultimate loads. The ultimate loads are the maximum loads at each node calculated from all the given loads, load combinations, and load factors. Also, the ultimate loads will be depicted in the plot window on the structure with red lines or curves. The line goes from the node in the direction of the force. A semicircular curve around the left side of a node is a positive moment according to the right-hand rule while a curve around the right side of a node is negative.

Find Displacements

At this point the structure is fully defined and graphically displayed in the plot window. Now when the user pushes the ‘Find Displacements’ button, the displacements vector will be displayed in the Scilab console and the new displaced shape of the structure will be plotted in blue. The displacements vector follows the same convention as discussed above with the loads matrix. The first 3 rows correspond to the horizontal, vertical, and angular displacements of node 1. The second 3 for node 2 and so on.

Find Reactions

When the user presses this button next, the reactions vector is displayed in the Scilab console. Again, as with the loads and displacements you have the horizontal, vertical, and moment reactions for each node starting with node 1 and going down.

Find Member Forces

Finally when the user pushes this button, the local member forces matrix is displayed in the Scilab console. Each column corresponds to a member. Column 1 for member 1 and so on. The first 3 rows are the horizontal, vertical, and moment forces on the member at the starting node while rows 4,5, and 6 are the forces on the member at the ending node. It is important to understand that these forces are with reference to the member’s local axis. The local axis for each member has its origin at the starting node with the horizontal or x-axis directed toward the ending node. The vertical or y-axis is then appropriately orthogonal to the local x-axis.

I hope this user’s guide has been helpful and provided sufficient information or answers to your questions. However if you are still stuck or want further information just let me know in the comments.

## The Bisection Method using Scilab

Recently, I have really been enjoying learning to use Scilab which is an open-source mathematics software that works alot like Matlab. Specifically we have run into some design problems in water and wastewater engineering that involved solving for variables that are described implicitly by an equation and cannot be solved through typical algebraic means. For example:

This cannot be solved directly for x. So how can we find the value of x that will satisfy the equation? Well this is where the numerical analysis technique of the bisection method comes in. The bisection method is a bounded or bracketed root-finding method. In other words, it will locate the root of an equation provided you give it the interval in which a root is located.

The search for the root is accomplished by the algorithm by dividing the interval in half and determining if the root is in one half or the other. Then it divides that interval in half and so on until it gets within a specified relative error. Below I have entered my code for this algorithm. If you save this file as Bisection_Method.sci in a scilab/modules/…/macros folder and execute it in the SciNotes editor, you can then use it by entering a command in the console like:

root = Bisection_Method(‘x-sinx-3’,0,6.5,0.5)

The first passed parameter is the function equation itself solved such that f(x) = 0. The second parameter is the lower end of the interval while the third is the upper end. The final parameter is the percent relative error that the algorithm much reach with the solution before it terminates and gives a result.

Please feel free to download Scilab for free here and copy this code and use it. Let me know if you have any troubles or if you find a problem with the code.

//The Bisection Method is a bracketed root locating method. Therefore
//the user must input the lower and upper bounds on the interval where
//the root is to be located.
//func_expression - must be a string containing the function expression
//with the variable written as 'x'.
//es - is the relative percent error criterion that must be met for the
//method to terminate.
function [root]=Bisection_Method(func_expression, x_lower, x_upper, es)
x = x_upper;
fu = evstr(func_expression); //evaluate the function at x_upper
x = x_lower;
fl = evstr(func_expression); //evaluate the function at x_lower
//***********************************************************
//test to see if there is a root in the interval
if (fu*fl >= 0) then
root = 'no root in interval given';
else
//there is a root in the interval
exact_solution = 'notfound';
ea = 100;
xr_new =(x_upper + x_lower)/2;
//*******************************************************
//iteratively progress toward to root until it is found or it
//is approximated with a relative error less than required
while ea > es & exact_solution == 'notfound',
x = xr_new;
fr = evstr(func_expression);
x = x_lower;
fl = evstr(func_expression);
//***************************************************
if(fl*fr < 0) then
x_upper = xr_new;
elseif(fl*fr > 0) then
x_lower = xr_new;
elseif(fl*fr == 0) then
root = xr_new;
exact_solution = 'found';
end //of if statement
//***************************************************
//calculate the approximate relative error for the iteration
xr_old = xr_new;
xr_new =(x_upper + x_lower)/2;
ea = abs((xr_new - xr_old)/xr_new)*100;
end //of while loop
//*******************************************************
//if error criterion has been met but the exact answer has not
//been found, set the root to the adequate approximation.
if (exact_solution == 'notfound') then
root = xr_new;
end //of if statement
//*******************************************************
end //of if-else statement
//***********************************************************
endfunction

## Calculus Applied to Manufacturing Design

Suppose you are in the can manufacturing industry. Your job is to create cans according to the specifications of canned goods suppliers who will put their products in your cans and distribute them for sale. Now suppose one of those canned good suppliers comes to you and says, “We need a can that can contain a 300 mL volume of our product. Joe Shmoe says he can provide us these cans for \$X. What can you do?” If you want the canned good suppliers business you will have to underbid Joe. (Provided you can’t sale them on some other facet of your company, i.e. service, quality, etc.) How will you do it? Where can you cut cost? What is the most inexpensive can you can make to satisfy the spec? In the real world this is a very complex question because it involves labor, facility, process, and material expenses along with several others. For the sake of this exercise however, we will use the principle of the derivative to optimize our can design in order to minimize material cost.

Now let’s assume you are set up to produce a standard cylindrical can with a top and a botttom. (I realize the canned goods supplier will have to add their product before the top is sealed and I realize that actual cans have a small lip on the top and bottom, but for the sake of simplicity…) To minimize the cost of the cans we want to minimize the amount of material used. But what shape cylinder will contain the required 300 mL but use minimal material? A short fat can or a tall thin can? How “tall” or how “short”?

First, realize the total material for the can can be determined by finding the area of the metal. For a cylinder this is the area of the bottom plus the area of the top plus the area of the side.

We can rewrite the area of the top and bottom as:

If you imagine removing the top and bottom, cutting the side along a vertical line, and laying it out flat you can see that the area of the side is the area of a rectangle with the length of one side equal to the height of the can and the length of the other side equal to the circumference of the top or bottom. So its area is:

Putting this all together, the total area of a can is given by:

Now the volume of the cylinder is given by:

If we rearrange this expression for the height, h we get:

Now we substitute this into our expression for the area to eliminate the h and make the whole thing in terms of 1 variable. (r – because the volume is given to us in the spec)

This equation tells us how much area a can has of a particular volume with any radius. With this in hand we can now use the derivative to optimize (maximize or minimize) the can design and find out which radius will give us the smallest area and thus the most inexpensive (as far as material costs go) cylindrical can that can contain the prescribed volume. Taking the derivative of the area equation using the power rule we have:

Now we set this equal to zero to find the minimum area. Solving for r we get:

Knowing that the volume is 300 mL our radius is:

Putting this back in our equation for h we get a height of:

These are the dimensions for our can of minimum material. Below I have tabulated and graphed some other values for the radius and associated heights that yield the necessary volume and show that the resulting surface area is larger than our optimized design.

One other beautiful thing about this technique is that I have a general formula for an optimized cylindrical can for any volume. I can just put in whatever the specification says and get my design.

Although this is greatly simplified, I hope this provides another demonstration of the real world applicability of these techniques.

## Basic Derivative Rules

Previously we have discussed how a derivative is found using the limit, why the derivative is useful, and how to find it. Now to attempt to simplify or shortcut the derivative finding process, we will look at some general rules that can be used instead of going through the entire limit step-by-step every time.

The first derivative we will look at is the derivative of a linear function in slope-intercept form.

Where m is the slope and b is the y-intercept. Using the limit definition for the derivative we have:

Therefore the derivative of any linear function is, as we would expect, equal to the slope of that line. The next derivative we will look at is a special case of the linear function. If the slope is zero then the linear function is a constant function.

As above, the derivative of a linear function is its slope therefore the derivative of a constant function is also equal to the slope which is zero.

The next rule we can look at is a single independent variable raised to any power. Hence this rule is often called the Power Rule.

Applying the limit process we obtain:

Notice how we used the Binomial Theorem to expand the first term in the numerator. At this point I want us to notice a certain pattern in the coefficients that can be seen in Pascal’s Triangle. The first coefficient will always be 1 (c1=1). This can be seen by looking at the left edge of the triangle. Every first number is 1. Secondly if you look at the diagonal row parallel to, directly adjacent to, and to the right of this first row of 1’s, you notice that the numbers count sequentially as you descend the horizontal rows of the triangle: 1, 2,3,4,5,..etc. These numbers correspond to the second coefficient in the binomial expansion. If we label the horizontal row of the very first 1 at the top of Pascal’s Triangle zero and then incrementally as we go down the triangle, the second row will provide the coefficients for the binomial squared, the third row for the binomial cubed, the forth row for the binomial to the fourth power, and so on such that the nth row gives the coefficients for the binomial to the nth power and the second coefficient will be n also (c2=n). From here we have:

So you can see that the derivative of any variable to an exponent is equal to that same variable times the original exponent to the power of one less than the original exponent. For example:

In summary we can use the following rules to make our lives easier with regard to finding derivatives.

Next time we’ll look at some more helpful derivative rules.