I am trying to solve a system of linear equations generated from a 2-D finite element problem. The equations (shown below) are not in the general form Ax=B..
4 views (last 30 days)
Show older comments
I am trying to solve a system of linear equations generated from a 2-D finite element problem. The equations (shown below) are not in the general form Ax=B. I have used syms and the solve function to obtain an answer but this strikes me as ineffecient and syms variables are annoying to deal with. Any advice on how this could be implemented without the use of syms variables would be apprecitiated.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1246542/image.png)
5 Comments
Star Strider
on 29 Dec 2022
The code you posted seems to work —
%Global Stiffness Matrix coefficients matrix
K_global=[1346153.84615385 0 -1346153.84615385 576923.076923077 0 0 0 -576923.076923077 0 0;
0 384615.384615385 384615.384615385 -384615.384615385 0 0 -384615.384615385 0 0 0;
-1346153.84615385 384615.384615385 3461538.46153846 0 -1346153.84615385 -384615.384615385 -769230.769230769 0 0 0;
576923.076923077 -384615.384615385 0 3461538.46153846 -576923.076923077 -384615.384615385 0 -2692307.69230769 0 0;
0 0 -1346153.84615385 -576923.076923077 1730769.23076923 0 0 961538.461538462 -384615.384615385 -384615.384615385;
0 0 -384615.384615385 -384615.384615385 0 1730769.23076923 961538.461538462 0 -576923.076923077 -1346153.84615385;
0 -384615.384615385 -769230.769230769 0 0 961538.461538462 2115384.61538462 0 -1346153.84615385 -576923.076923077;
-576923.076923077 0 0 -2692307.69230769 961538.461538462 0 0 3076923.07692308 -384615.384615385 -384615.384615385;
0 0 0 0 -384615.384615385 -576923.076923077 -1346153.84615385 -384615.384615385 1730769.23076923 961538.461538462;
0 0 0 0 -384615.384615385 -1346153.84615385 -576923.076923077 -384615.384615385 961538.461538462 1730769.23076923];
syms Y1 Y2 X3 Y3 X5 u1 u2 u4 v4 v5
%Force Vector Initialization
Force = [0 ; Y1 ; 0 ; Y2 ; X3 ; Y3 ; 0 ; 1000 ; X5 ; 1000];
%Disp Vector
Displacement = [u1 ; 0 ; u2 ; 0 ; 0 ; 0 ; u4 ; v4 ; 0 ; v5] ;
%Solving the system equation F=
Solved=solve(K_global*Displacement==Force,[Y1 Y2 X3 Y3 X5 u1 u2 u4 v4 v5])
vpaSolved = vpa(struct2cell(Solved))
.
Accepted Answer
Karim
on 29 Dec 2022
Edited: Karim
on 29 Dec 2022
Well to figure out the simple analytic equation you can in fact use the symolic toolbox to investigate the set of equations. Notice how you can seperate you degrees of dreedom in 2 groups: for some you know that the discplacment is zero (i.e. they are constrained) for others you know the applied forces. Hence you could split the DOF's in two groups a and b. Using the symbolic toolbox this leads to:
syms Kaa Kab Kbb Ua Ub Fa Fb real
K = [Kaa Kab; Kab' Kbb]
U = [Ua; Ub]
F = [Fa;Fb]
EOM = K*U == F
Lets say that Ub are the known displacments and Fa the known forces. Since in this case Ub = 0, we can reduce the equations to simply:
Ua = Kaa\Fa
Fb = Kab*Ua
If we apply this to your system we obtain:
% indexes for the known forces
a = [1 3 7 8 10];
% indexes for the known displacements (in this case constraints)
b = [2 4 5 6 9];
% stiffness matrix of the system
K = [1346153.84615385 0 -1346153.84615385 576923.076923077 0 0 0 -576923.076923077 0 0;
0 384615.384615385 384615.384615385 -384615.384615385 0 0 -384615.384615385 0 0 0;
-1346153.84615385 384615.384615385 3461538.46153846 0 -1346153.84615385 -384615.384615385 -769230.769230769 0 0 0;
576923.076923077 -384615.384615385 0 3461538.46153846 -576923.076923077 -384615.384615385 0 -2692307.69230769 0 0;
0 0 -1346153.84615385 -576923.076923077 1730769.23076923 0 0 961538.461538462 -384615.384615385 -384615.384615385;
0 0 -384615.384615385 -384615.384615385 0 1730769.23076923 961538.461538462 0 -576923.076923077 -1346153.84615385;
0 -384615.384615385 -769230.769230769 0 0 961538.461538462 2115384.61538462 0 -1346153.84615385 -576923.076923077;
-576923.076923077 0 0 -2692307.69230769 961538.461538462 0 0 3076923.07692308 -384615.384615385 -384615.384615385;
0 0 0 0 -384615.384615385 -576923.076923077 -1346153.84615385 -384615.384615385 1730769.23076923 961538.461538462;
0 0 0 0 -384615.384615385 -1346153.84615385 -576923.076923077 -384615.384615385 961538.461538462 1730769.23076923];
% known forces
Fa = [0;0;0;1000;1000];
% known displacements (which we actually don't need since they are zero)
Ub = [0;0;0;0;0];
% compute the unknown displacements
Ua = K(a,a) \ Fa
% compute the unkown forces
Fb = K(b,a) * Ua
% print things a bit more pretty
fprintf('Y1 = %.3f, Y2 = %.3f, X3 = %.3f, Y3 = %.3f, X5 = %.3f\n',Fb)
fprintf('u1 = %.6f, u2 = %.6f, u4 = %.6f, v4 = %.6f, v5 = %.6f\n',Ua)
Notice that these results are the same as what you obtained with the symbolic toolbox.
More Answers (1)
John D'Errico
on 29 Dec 2022
Edited: John D'Errico
on 29 Dec 2022
Ok. Now I see a real problem. I almost was going to make one up. But that is always less fun, and takes just a little more thought on my part. So here is your problem.
K_global=[1346153.84615385 0 -1346153.84615385 576923.076923077 0 0 0 -576923.076923077 0 0;
0 384615.384615385 384615.384615385 -384615.384615385 0 0 -384615.384615385 0 0 0;
-1346153.84615385 384615.384615385 3461538.46153846 0 -1346153.84615385 -384615.384615385 -769230.769230769 0 0 0;
576923.076923077 -384615.384615385 0 3461538.46153846 -576923.076923077 -384615.384615385 0 -2692307.69230769 0 0;
0 0 -1346153.84615385 -576923.076923077 1730769.23076923 0 0 961538.461538462 -384615.384615385 -384615.384615385;
0 0 -384615.384615385 -384615.384615385 0 1730769.23076923 961538.461538462 0 -576923.076923077 -1346153.84615385;
0 -384615.384615385 -769230.769230769 0 0 961538.461538462 2115384.61538462 0 -1346153.84615385 -576923.076923077;
-576923.076923077 0 0 -2692307.69230769 961538.461538462 0 0 3076923.07692308 -384615.384615385 -384615.384615385;
0 0 0 0 -384615.384615385 -576923.076923077 -1346153.84615385 -384615.384615385 1730769.23076923 961538.461538462;
0 0 0 0 -384615.384615385 -1346153.84615385 -576923.076923077 -384615.384615385 961538.461538462 1730769.23076923];
syms Y1 Y2 X3 Y3 X5 u1 u2 u4 v4 v5
%Force Vector Initialization
Force = [0 ; Y1 ; 0 ; Y2 ; X3 ; Y3 ; 0 ; 1000 ; X5 ; 1000];
%Disp Vector
Displacement = [u1 ; 0 ; u2 ; 0 ; 0 ; 0 ; u4 ; v4 ; 0 ; v5];
Solve will handle it trivially of course, but you don't want that. Picky, picky, picky. ;-)
The trick is to split the vectors into two pieces, making one of the pieces entirely constant, and the other one entirely variable. You have a square linear system of 10 unknowns.
size(K_global)
First, let me break the two vectors down into their respective components. The set of all unknowns is this:
unknowns = [Y1 Y2 X3 Y3 X5 u1 u2 u4 v4 v5].';
Can we extract only the first 5 of them? Of course. I'm doing this part symbolically, just to show you how it works.
A = zeros(10);
A(2,1) = 1;A(4,2) = 1;A(5,3) = 1;A(6,4) = 1;A(9,5) = 1;
A*unknowns
Similarly, I'll build a second 10x10 array that can extract the last 5 unknowns.
B = zeros(10);
B(1,6) = 1;B(3,7) = 1;B(7,8) = 1;B(8,9) = 1;B(10,10) = 1;
B*unknowns
Again, that was just to show you what these matrices are doing. Next, I'll build two constant vectors.
Con1 = [0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 1000 ; 0 ; 1000];
Con2 = [0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0]; % All zeros, but still...
Now, let me convince you that these combine as you have them.
Con1 + A*unknowns
You should see this worked.
Now, what were you trying to solve? I'll expand the problem you had.
K_global*Con1 + K_global*A*unknowns = Con2 + B*unknowns
Swap some terms around...
K_global*A*unknowns - B*unknowns = Con2 - K_global*Con1
Do you see where this is heading? I'll do one more operation.
(K_global*A - B)*unknowns = Con2 - K_global*Con1
Finally, lets solve the problem, using NO symbolic algebra.
format long g
unk_num = (K_global*A - B)\(Con2 - K_global*Con1)
So all 10 unknowns, solved using simple linear algebra. Just a call to backslash. The only step of any importance was to know how to construct the matrices A and B.
Does it agree with solve?
sol = solve(K_global*Force == Displacement)
double(sol.Y1)
double(sol.v4)
That is about as close as backslash can resolve the solution in double precision arithmetic, because the numerical condition number of the linear system is moderately large. That means you should expect to see only about 8 or 9 significant digits in the solution, below the maximum element of the solution. And since we have Y1 == 1000, I would expect to start seeing problems around 1e-6 in the result.
cond(K_global*A - B)
You could probably resolve SOME of that by rescaling the variables, IF it was important to you.
0 Comments
See Also
Categories
Find more on Equation Solving in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!