## 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

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

Once again, we divide out , rearrange, and are left with

In order to get this into the form necessary to use the Lambert W function, we multiply both sides by .

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.

This means the eigenvalues of the system can be found as

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

Well it turns out that, in general, the matrices and do not commute. In other words, in general,

(it’s true, run some numbers!!)

This is due to the fact that,for a general case,

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.

## Useful Soil Coefficients

Previously we discussed mechanical analysis and the method of sieve analysis. This allowed us to develop the particle size distribution curve. With this curve we can find some useful parameters. They are:

- The Uniformity Coefficient
- The Coefficient of Gradation

Both of these parameters are used in soil classification. The Uniformity Coefficient is defined as:

The Coefficient of Gradation is defined as:

where

These diameters are obtained from the particle size distribution curve by going across from the percent finer, coming to the curve, and turning 90 degrees down to the abscissa to locate the diameter as seen below.

**Related Posts:**

## Sieve Analysis

The first method of mechanical analysis is sieve analysis. It is pretty straight forward and involves straining a soil sample through a series of sieves. Generally the sieves are stacked as seen in the above diagram with a US Standard Sieve No. 4 at the top (This sieve size has a 4.75 mm opening) and a Sieve No. 200 at the bottom (with a 0.075 mm opening). The stack usually consists of 5 to 7 sieves and has a pan at the bottom. As you can see, the opening in the sieve gets smaller as you go down. Therefore larger particles will be captured while smaller particles will pass through. This allows us to measure how much soil of each size was in the original soil sample.

A few important steps:

- the soil sample needs to be oven-dried and all clumps broke.
- the weight of the empty sieves needs to be obtained.

Once this is done, the soil sample is placed into the top sieve, the lid is placed on the stack and the stack is placed in a mechanical shaker. The stack will then be shaken vigorously for approximately 5 minutes to ensure proper distribution of the particles.

The next step is to calculate what is called the *percent passing* or the *percent finer*. This is simply the percentage of the entire soil sample that passed through a particular sieve. For example, the percent passing the No. 200 sieve would be the mass of the soil in the bottom pan (since that is the soil that passed through the No. 200 sieve) divided by the total mass of the original soil sample. So in general we can write:

where

- %P is the percent passing
- i represents the ith sieve
- n represents the number of sieves
- m is the mass

Once the percent passing is determined for each sieve, a plot is made of sieve or particle size versus percent passing as seen below. The percent passing is on the ordinate and the sieve or particle size is on the abscissa. Also notice the abscissa is reverse logarithmic.

In future posts we will look more carefully at what information this distribution curve can provide, but for now we are highlighting the process. If you have any questions, please let me know in the comments.

**Related Posts:**

## 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

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

Functions of matrices can generally be evaluated using the following relationship

where is a matrix, is a matrix of eigenvectors, and is the Jordan canonical form of the matrix . 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 with eigenvalues and and let’s assume that is repeated once and is repeated twice. This means that, since there are two repeated eigenvalues, there are two Jordan blocks. Then, can be written as

Each Jordan block accounts for a single repeated eigenvalue. Now, the expression is

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

where is the duplicity of the eigenvalue. If there are no repeated eigenvalues, the matrix 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 can be written as

Differentiating the defining equation for and solving for (implicit differentiation), we obtain the following expression for the derivative of

Further derivatives of can be taken and a general formulation of the th derivative of 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

,

where and are matrices, the solution, , is valid only when and commute; i.e. when . The proof is quite simple and satisfactory. First, assume that . It follows that

since will always commute with its own matrix exponential, . 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.

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.

In this case 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

Thus, the original equation becomes

If we assume a solution of the form , the equation becomes

Collecting terms and dividing out

Clearly cannot equal zero, thus the portion of the equation in the parentheses must equal zero. Thus, after re-arranging,

Thus, using the Lambert W function,

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

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

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

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:**

## 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:

This type of equation is associated with a system where the output is equivalent to the input delayed by a small time constant and multiplied by a system constant . 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.

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

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

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

The solution to this equation is:

In other words,

where is the Lambert W function of and 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 or ), 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. ).

Returning to the Laplace transform equation,

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

Multiplying both sides by yields:

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:

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

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

where is a vector and is a matrix. Stay tuned for more!

## Mesh Circuit Analysis

Recently I had a reader comment with a request for an example of Kirchoff’s Voltage Law (KVL) with Ohm’s Law. KVL and Ohm’s Law are both used in the circuit analysis methodology called Mesh Circuit Analysis. Consequently I thought I would analyze a small circuit and demonstrate the general principles of this technique that can be used on any circuit.

The circuit we will use below has a 10 volt source and 4 resistors configured as seen.

**STEP 1:**

The first step is to label the circuit as I have done in red ink. This includes the current flowing clockwise through each loop (or mesh) and the voltage drops across each passive element. (In this case all the passive elements are resistors.) Next we will use the picture to write some equations.

**STEP 2:**

If you look at the diagram at the first loop and start at the top, then travel around the loop going clockwise and add up the voltage drops you get the first equation below. Notice the voltage drop across the source is negative because it is not a drop it is an increase.Also notice the sum must equal zero.

The next 2 lines are simply writing the unknown voltage drops in terms of Ohm’s Law. The voltage drop across the resistor is equal to the current through the resistor times the resistance. Notice for resistor 2 the currents are opposing one another. In other words current 1 runs from the top of the circuit toward the bottom of the circuit through resistor 2. Current 2 runs from the bottom to the top. That is why the current component in the second Ohm’s Law equation is i1 minus i2.

Next these two Ohm’s Law equations are substituted in the original equation which gives the 4th line below.

From the 4th line the rest is just algebra. Then do the exact same process for loop 2 as seen above. This gives 2 equations with 2 unknowns. Of course these can be solved by any linear algebra method you prefer.

STEP 3:

Below I solved it using simple substitution. This solves for the unknown currents.

STEP 4:

Finally, substitute these found currents back into the Ohm’s Law equations we developed earlier and you find the unknown voltage drops. The last diagram shows the final solution.

I hope this answers my reader’s question and has been helpful. Please let me know if anyone needs more help.

## Newton’s Second Law

From my education and experience, I would venture to say that the vast majority of engineering mechanics can have its origins traced back to Newton’s Second Law. It is amazing how much this relatively simple little law can be developed and tell you a great deal about how the world works. In a previous post on kinematics, we discussed the relationship between position, velocity, and acceleration. There we said that those relationships could be developed independent of the forces that created the various motions. Well the study of kinetics is the study of the forces which cause motion and how that all works and it of course is based on Newton’s Second Law. The form I would like to state Newton’s Second Law in is this:

The summation symbol is to show that it is the net force that creates motion or the sum of the forces. For instance if I have a ball on the table and I push in one direction with a certain force and also simultaneously push in the exact opposite direction with the same force, the ball will not move. (It may deform, but no motion.) Now I applied two forces, but nothing happened. That is what is meant by the net force or sum of forces. There must be an unbalanced total force to create motion. The differential on the right-hand side is the rate of change of linear momentum with respect to time. Now I know that some, not well acquainted with calculus may say that this doesn’t look that simple, but I mean it only has one term on the left and one on the right. A equals B. That’s it. The motion of almost all things (at least all the everyday things you see. Some variation lies with subatomic particle motion, and things moving at close to the speed of light, but that’s another topic.) Is summed up in this?! Amazing!

Now if we make a few modifications we can manipulate the equation into the form that most people are familiar with. Now linear momentum is the product of the mass and velocity of a body.

Well if we assume that the mass of the body is constant and does not change with time as the body moves, and realizing that the time rate of change of velocity is the definition of the instantaneous acceleration, you have:

This is typical form you see. Finally we will divide by the mass to change it to this form:

What you can understand from this is that as the sum of the unbalanced forces increases, the acceleration increases, and as the mass of a body increases, the acceleration decreases. This is of course all of our real world experience. The harder something is pushed, the more it speeds up. The bigger it is, the harder it is to speed up.

Already you can see that we can extrapolate a lot of understanding from this little relationship. I plan to build on this in future posts to show some of its application to the different areas of engineering.

Also, here are a few links to some additional helps. First is a physics lecture on Newton’s Laws provided by MIT OCW. Second is a interactive graphical demonstration of Newton’s Second Law provided by Wolfram Demonstrations. Finally, a Khan Academy video.

## 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 Num_Load_Types; global Loads; global Num_Load_Combos; global Load_Factors; global Ultimate_Loads; 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 //modify the menus of 'Main_Window' delmenu(Main_Window.figure_id, gettext("&File")); delmenu(Main_Window.figure_id, gettext("&Tools")); delmenu(Main_Window.figure_id, gettext("&Edit")); delmenu(Main_Window.figure_id, gettext("&?")); toolbar(Main_Window.figure_id, "off"); // create the "File" menu on the menu bar file_menu = uimenu(Main_Window,'label', '&File'); //create two submenus in the "File" menu new_submenu=uimenu(file_menu,'label', '&New'); close_submenu=uimenu(file_menu,'label', '&Close', 'callback', "Close_Program()"); //create submenus of "New" new_select1=uimenu(new_submenu,'label', '&Beam','callback',"Beam_Analysis()"); new_select2=uimenu(new_submenu,'label', '&Frame','callback',"Frame_Analysis()"); new_select3=uimenu(new_submenu,'label', '&Truss','callback',"Truss_Analysis()"); endfunction //======================================================================================================================== //MENU-ROUTINES //======================================================================================================================== 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(); Load_Controls(); 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(); Load_Controls(); 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(); Load_Controls(); endfunction //======================================================================================================================== //SUB-ROUTINES //======================================================================================================================== function Clear_Figures() scf(2); clf(); scf(1); f=gcf(); C=f.children; for i=1:length(C) c=C(i); if c.type~="uimenu" delete(c); end end endfunction //======================================================================================================================== //======================================================================================================================== function Load_Controls() 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()"; //==================================================================================================================== //4th Line: applied loads //creates text in the main window Loads_Text=uicontrol(gcf(),"style","text"); Loads_Text.Position = [margin_x,main_window_h-margin_y-3*text_h-3*spacing_y,text_w,text_h]; Loads_Text.BackgroundColor = [1,1,1]; //white Loads_Text.String = "Number of Load Types" //creates edit box for entering data in the main window Num_Loads_Box=uicontrol(gcf(),"style","edit"); Num_Loads_Box.Position = [margin_x+text_w+spacing_x,main_window_h-margin_y-3*text_h-3*spacing_y,text_h,text_h]; Num_Loads_Box.String = "0"; Num_Loads_Box.Tag = "NLTB" //create a pushbutton in the main window Loads_Button=uicontrol(gcf(),"style","pushbutton"); Loads_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y-3*text_h-3*spacing_y,button_w,text_h]; Loads_Button.String = "Loads"; Loads_Button.Callback = "Get_Loads()"; //========================================================================================================================= //5th Line: load factors and load combinations //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]; Combos_Button.String = "Load Factors"; Combos_Button.Callback = "Get_Load_Combos()"; //========================================================================================================================= //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 Menu = list('',1,['Singularly','By Sets','From File']); Node_Entry_Method = x_choices('Select Entry Method',list(Menu)); //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 not then adjust 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 Menu = list('',1,['Singularly','By Sets','From File']); Member_Entry_Method = x_choices('Select Entry Method',list(Menu)); //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 //============================================================================================================================ //============================================================================================================================ function Get_Loads() //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_Load_Types; global Loads; global Num_Nodes; //Extract Num_Load_Types from Edit Box in Main_Window Num_Load_Types = evstr(get(findobj("Tag","NLTB"),"String")); Loads = zeros(3*Num_Nodes,Num_Load_Types); for k=1:Num_Load_Types //initializes dynamically sized matrix for loads entry and associated //labels Row_Labels=zeros(Num_Nodes); Loads_Entry=zeros(Num_Nodes,3); for i=1:Num_Nodes Row_Labels(i) = i; end //displays a matrix dialog box for user inputs Loads_Entry =evstr(x_mdialog('Enter the Loads Applied to each Node for Load Type '+string(k),string(Row_Labels),['Horizontal Load','Vertical Load','Applied Moment'],string(Loads_Entry))); for i=1:Num_Nodes //for each node for j=1:3 //for each applied load Loads(3*(i-1)+j,k)=Loads_Entry(i,j); end end end //print out the loads matrix print(6,Loads); endfunction //===================================================================================================================== //===================================================================================================================== function Get_Load_Combos() //obtains the number of load combinations from the user and provides an x_mdialog //box for the user to input the load factors for each load type for each load //combination global Num_Load_Combos; global Load_Factors; global Num_Load_Types; global Loads; global Num_Nodes; global Axes_Size; global Node_Coord; global Ultimate_Loads; global Num_Members; global Connect_Mat; //Extract Num_Load_Combos from Edit Box in Main_Window Num_Load_Combos = evstr(get(findobj("Tag","NLCB"),"String")); //initializes dynamically sized matrix for members and associated //labels Load_Factors = zeros(Num_Load_Types,Num_Load_Combos); for i=1:Num_Load_Types Row_Labels(i)= 'Load Type '+string(i); end for i=1:Num_Load_Combos 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' Load_Factors =evstr(x_mdialog('Enter the Load Factors',string(Row_Labels),string(Column_Labels),string(Load_Factors))); //Apply the load combinations to the loads and find the ultimate load at each node Combined_Loads = Loads*Load_Factors; Ultimate_Loads = zeros(3*Num_Nodes,1); for i=1:(3*Num_Nodes) if max(Combined_Loads(i,:))==max(abs(Combined_Loads(i,:))) then Ultimate_Loads(i,1) = max(Combined_Loads(i,:)); else Ultimate_Loads(i,1) = -max(abs(Combined_Loads(i,:))); 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 if Ultimate_Loads(3*(i-1)+j,1)~=0 then select j case 1 then //horizontal forces X(1)=Node_Coord(i,1); Y(1)=Node_Coord(i,2); X(2)=Node_Coord(i,1)+unit_size*(Ultimate_Loads(3*(i-1)+j,1)/abs(Ultimate_Loads(3*(i-1)+j,1))); 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); Y(2)=Node_Coord(i,2)+unit_size*(Ultimate_Loads(3*(i-1)+j,1)/abs(Ultimate_Loads(3*(i-1)+j,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; if Ultimate_Loads(3*(i-1)+j,1)<0 then a1 = -80*64; a2 = 160*64; elseif Ultimate_Loads(3*(i-1)+j,1)>0 then 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 //print out the load factors matrix and final ultimate loads vector print(6,Load_Factors); print(6,Ultimate_Loads); 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 Ultimate_Loads; global Displacements; global Num_Load_Types; 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; ultimate_loads = Ultimate_Loads; //delete rows and columns corresponding to given parameters and boundary conditions k(Eliminate,:)=[]; k(:,Eliminate)=[]; Index(Eliminate,:)=[]; ultimate_loads(Eliminate,:)=[]; //calculate the unknown displacements Unknown_Displacements=inv(k)*ultimate_loads; //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