Solve a system of symbolic variables

Hi,

I have been trying hard at this for a while now. Does anybody know if there is some way to solve symbolic variables in a matrix, resulting in a double format?

I want to solve for the unknowns of F and d. The amount of unknowns depends on the values specified in the input:

% This program calculates any straight equal-length element analysis
% for a bar with two fixed ends

syms Freact1 Freact2;

A = input('Cross-sectional area (mm^2): ');
E = input('Elastic Modulus (MPa): ');
L = input('Length of the system (mm): ');

num_ele = input('Enter the number of elements to be analyzed: ');
num_nodes = num_ele + 1;

k = (E*A)/(L/num_ele);

F = zeros(num_nodes,1);
d = zeros(num_nodes,1);

d(1) = 0; d(num_nodes) = 0;
F1 = vpa(Freact1); Fend = vpa(Freact2);
F(1) = F1; F(num_nodes) = Fend;
disp(' ');

for i = 2:num_ele
    str = sprintf('Node %d: ', i);
    disp(str);
    F(i) = input('Enter the force at the node (N): ');
    disp(' ');
end

glbl_stiff = 0*diag(num_nodes,num_nodes-1) + 2*eye(num_nodes);
glbl_stiff(1,1) = 1; glbl_stiff(num_nodes,num_nodes) = 1;

for j = 1:num_ele
    glbl_stiff(j,j+1) = -1;
    glbl_stiff(j+1,j) = -1;
end

glbl_stiff = k*glbl_stiff;

F = glbl_stiff*d
d = F\glbl_stiff

and the error that always results:

The following error occurred converting from sym to
double:
Error using mupadmex
Error in MuPAD command: DOUBLE cannot convert the
input expression into a double array.

If the input expression contains a symbolic
variable, use the VPA function instead.

Error in FEM_IA1 (line 20)
F(1) = F1; F(num_nodes) = Fend;

It's confusing because I did use the VPA function. It just seems like nothing's working.. any ideas would greatly be appreciated.

Thanks, Ian

Answer 1

[Star Strider]

The variables Freact1 and Freact2 in that line are symbolic variables. You haven't assigned any numeric values to them.

[Ian Wood]

yes i'm trying to solve for the symbolic variables in my program, but to my knowledge it is not achievable

[Star Strider]

I don't see Freadt1 or Freact2 in any equation anywhere. What are you solving for?

[Ian Wood]

they're in matrix F, and I want to solve for them

[Star Strider]

Even in symbolic variables, array references have to be numeric. You don't need vpa there since there's nothing to evaluate — you simply need to define Freact1 and Freact2 as numbers in order to use them the way you want to, since you've already defined F as a numeric vector.

Add — I also just noticed that you define F three times — once in the line that is throwing the error, once in your i for loop, and again at the end where you defined it in terms of glbl_stiff*d.

I'm lost!

[Ian Wood]

well the first and last elements of the matrix are supposed to be the ones i'm solving for, and the other F's, whatever the size, are from user inputs. If I initialize the variables as numbers, I will not be able to solve for them later because they already have a value.. and that's where i'm lost.

[Star Strider]

Please post the equation — and any necessary additional information — that you are using in your problem. A Wikipedia or other link would work if it adequately describes what you are doing.

[Ian Wood]

The equation is Hooke's Law: F = k*x where F, k and x are matrices of size depending on the input specified in the number of elements. F is a force column matrix, and the same goes for the displacement matrix. k is a stiffness square matrix that contains all known values. I have provided a link to show what i'm doing (please see example 2.1 on page 42):

http://www.scribd.com/doc/62642917/A-First-Course-in-the-Finite-Element-Method

[Star Strider]

Thank you for the reference. I finally figured out what you want to do, even though my experiences with the finite element method are not recent. I apologise for the delay, but it took me a few minutes to read the reference and put your code into context.

I made some changes, so I'm appending my version of your code here. See if it does what you want:

% This program calculates any straight equal-length element analysis
% for a bar with two fixed ends
syms Freact1 Freact2;
A = input('Cross-sectional area (mm^2): ');
E = input('Elastic Modulus (MPa): ');
L = input('Length of the system (mm): ');
num_ele = input('Enter the number of elements to be analyzed: ');
num_nodes = num_ele + 1;
k = (E*A)/(L/num_ele);
F = zeros(num_nodes,1);
d = zeros(num_nodes,1);         % You calculate ‘d’ later, so you don't need to define it (unless you want to preallocate it)
d(1) = 0; d(num_nodes) = 0;
% F1 = vpa(Freact1); Fend = vpa(Freact2);   % Not necessary
% F(1) = F1; F(num_nodes) = Fend;           % Not necessary
disp(' ');
for ki = 2:num_ele
    str = sprintf('Node %d: ', ki);
    disp(str);
    F(ki) = input('Enter the force at the node (N): ');
    disp(' ');
    fprintf(1,'\n\tF(%d) = %f\n', ki, F(ki))                        % DIAGNOSTIC
end
fprintf(1,['\nF = ' repmat('\n\t\t%f\n',1,size(F,1)) '\n\n'], F)    % DIAGNOSTIC
glbl_stiff = 0*diag(num_nodes,num_nodes-1) + 2*eye(num_nodes);
glbl_stiff(1,1) = 1; 
glbl_stiff(num_nodes,num_nodes) = 1;
for j = 1:num_ele
    glbl_stiff(j,j+1) = -1;
    glbl_stiff(j+1,j) = -1;
end
glbl_stiff = k*glbl_stiff;
% F = glbl_stiff*d              % PLEASE don't redefine ‘F’ here!
% d = F\glbl_stiff              % Correct idea, but wrong order —
%                                 ‘glbl_stiff’ PRE-multiplies ‘d’ so the correct 
%                                 way to calculate ‘F’ is to do the equivalent 
%                                 of: F = inv(glbl_stiff)*d.  
%                                 PLEASE do not use ‘inv’ in practice!  

% d = glbl_stiff\F              % Correct idea, but ...
d = lsqr(glbl_stiff,F)          % ... since matrix ‘glbl_stiff’ is sparse, you need to use sparse techinques with it

You don't need to use the Symbolic Math Toolbox in this situation. You can do everything you need to without it. You can comment or delete the lines I labeled % DIAGNOSTIC since they were for my benefit.

I'll keep this open for a while to be sure you and I are converging on a solution. If this works as you want it to, please follow up with a comment so I'll know it worked.

[Ian Wood]

That seemed to do it the proper way, but what if I have any unknown F values in my vector? How could I then solve for d?

[Star Strider]

My approach to unknown F values in your vector would be to reduce the size of glbl_stiff to accomodate them. Two connected springs with an unknown force at the node would probably have to be lumped as a single spring. You would have to reduce the order of your system in that situation, unless you wanted to do the solution symbolically, accounting for the unknown F-values as symbolic variables. I can't think of any other way to approach that situation. (My background is EE not ME.)

[Ian Wood]

Yeah usually when solving these by hand I would cancel out the rows and respective columns of all the boundary conditions (basically displacements and forces are zero). I have to apply some code to my program that eliminates such rows and columns, and obtain a smaller matrix that contains what's left.

[Star Strider]

See the documentation for sparse matrices. There are a number of functions — particularly find (link at the end of the sparse page) that can help you find and eliminate the zero rows and columns.

Did this solve your problem?

[Ian Wood]

I need to experiment with these functions a little bit before I can confirm. Thanks a lot for your help to date.

[Star Strider]

I'll keep this open until we're happy we've converged on a solution then.

Answer 2

[Matt Fig]

Use symbolic arrays instead:

F = sym(zeros(num_nodes,1));
d = sym(zeros(num_nodes,1));

[Ian Wood]

Never really thought about doing this, thanks. I'm trying another method for the time being, but if that doesn't work i'll revert to this.

[Phan HaNhut]

if memory of computer is slow, code can not end? (matlab: "busy") And there is not any result?

I have to add more RAM?

[Walter Roberson]

If your memory is limited, then you can speed up operations by turning off virtual memory (or configuring it to be size 0). Swapping memory to disk is very very slow, and my practical experience on MS Windows systems is that once you swap enough program memory to disk then you cannot make any progress because you run into "thrashing" (part of memory you need for the calculation is swapped out in order to bring in something else you need, but then that has to get swapped out in order to bring the first back in in order to proceed, but then that needs... etc.) Turning off virtual memory would result in the calculation failing cleanly with complaints about insufficient memory instead of running for days getting nowhere.

Expanding physical memory is usually good. You may need to switch to a 64 bit operating system (with the extra memory) to make real progress.

Answer 3

[Ian Wood]

OK, so working through this I came up with another code, and I provided some comments to explain what I did. The code changed quite a bit, but the desired end result is still the same.

Still don't know how to solve for the unknowns (F and d), but I've figured out how to eliminate the entries that are not needed to solve the equations.

Here's my updated code:

 % This program calculates any series-element analysis for a bar

 A = input('Cross-sectional area (mm^2): ');
 E = input('Elastic Modulus (MPa): ');
 L = input('Length of the system (mm): ');

 num_ele = input('Enter the number of elements to be analyzed: ');
 num_nodes = num_ele + 1;

 % initialize the element length vector

 ele_length = zeros(num_ele,1);

 % input the length of each element, and if the sum is not equal to
 % specified input, quit the program

 for n=1:num_ele
     disp(' ');
     str = sprintf('Element %d: ', n);
     disp(str);
     ele_length(n) = input('Enter the length of the element (mm): ');
 end 

 if sum(ele_length) ~= L
     error('error: lengths are not equal');
 end

 k = (E*A)./(ele_length);    % stiffness constant for each element length

 % initialize force, displacement, and degree of freedom vectors

 F = zeros(num_nodes,1);
 d = zeros(num_nodes,1);
 DOF = zeros(num_nodes,1);

 % prompt for nodal forces and degrees of freedom for each node. A degree of
 % freedom of zero means a fixed location, and thus displacement of zero.

 for m = 1:num_nodes
     disp(' ');
     str = sprintf('Node %d: ', m);
     disp(str);

     F(m) = input('Enter the force at the node in newtons: ');
     DOF = input('Enter the degrees of freedom (0 or 1): ');

     if DOF == 0
         d(m) = 0;

     % if there is a degree of freedom, there exists displacement

     elseif DOF == 1
         d(m) = input('Enter the displacement at the node in mm: ');      

     else
         error('error: invalid input');
     end
 end

 % obtain the global stiffness matrix

 glbl_stiff = 0*diag(num_nodes,num_nodes-1) + 2*eye(num_nodes);
 glbl_stiff(1,1) = 1; glbl_stiff(num_nodes,num_nodes) = 1;

 for j = 1:num_ele
     glbl_stiff(j,j+1) = -1;
     glbl_stiff(j+1,j) = -1;
 end

 glbl_stiff = sum(k)*glbl_stiff;

 ind = find(d==0);   % find the indices where the elements are zero

 % remove elements corresponding to boundary condition

 glbl_stiff(ind,:) = [];
 glbl_stiff(:,ind) = [];
 F(ind) = [];
 d(ind) = [];

 % d = lsqr(glbl_stiff,F)
 % F = glbl_stiff*d

[Star Strider]

Maybe I'm missing something here, but you're inputing F and d to generate glbl_stiff.

When you're finished, there's nothing left to solve for!

If you later want to input different F or d vectors with your defined glbl_stiff and solve for the other vector, you can then do it with lsqr.

NOTE: ‘Answers’ was down for a long time yesterday afternoon and evening, so I didn't check back in until this morning.

[Ian Wood]

Well I do want to input some F and d values, but not all values. Certain values will be desired. glbl_stiff will always have completely known values. As an example for a two-element analysis, F could be [F1; 1000; F3] and d [0 d2 0]. I would need to find F1, F3, and d2 with the use of glbl_stiff.

This is why I thought that the use of symbolic variables could be implemented, but apparently there is no need to use them to solve this system (according to my professor).

[Star Strider]

Scanning the book you referred me to earlier (the discussions on pages 46-7), it seems that you would solve independently for the forces with zero (or other defined) displacements and then for the displacements with known forces as separate procedures.

So you would for instance solve for d2 as a linear equation with a force of 1000, then use sparse matrix techniques to solve for F1 and F3 with displacements = 0. That seems to be what the book suggests, and that's certainly how I would do it. It doesn't seem mathematically possible to do both simultaneously, since that would mean having unknowns on both sides of the equation.