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

[John D'Errico]

Sigh. By posting a picture of the equations, you make it impossible for someone to paste those coefficients into MATLAB. And while I would have answered your question quickly if I could have just pasted the coefficient array into MATLAB. there is no way I'd touch it as is. (The problem IS simple to solve.)

Is there a good reason why you want to make it more difficult to get help?

[Simon Cook]

Just wanted a general idea of where to look/ strategies to employ. Appreciate your willing to put it in matlab but I hope if sent in the right direction could get there myself

[Simon Cook]

This is a coefficients matrix that is solved using syms (this a slightly different problem to the one above.)

%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])

[Star Strider]

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])

Solved = struct with fields:
    Y1: -4036192649312651740175608838454866923095997416655529980658973193831841978026000/186241460818286368935039332151175406902707246982877619366601958753940483110093
    Y2: -206422424064849525457419958683837935439876853243717067371307579023151209804542000/186241460818286368935039332151175406902707246982877619366601958753940483110093
    X3: -28253348545189310078497994185890367976918203945359452972264142116095063643083200/186241460818286368935039332151175406902707246982877619366601958753940483110093
    Y3: -162024304922410839458520061544453002820777764275394239208805460147645336382137550/186241460818286368935039332151175406902707246982877619366601958753940483110093
    X5: 28253348545188961761750420455665936982195277584804555066055997059445158833591900/186241460818286368935039332151175406902707246982877619366601958753940483110093
    u1: 87850616007040499973187133460274492118576438093866386504140615856422912000/186241460818286368935039332151175406902707246982877619366601958753940483110093
    u2: 46923622543009790497369096096733171226688088298909721848612942018877849600/186241460818286368935039332151175406902707246982877619366601958753940483110093
    u4: 57417723431222675248432621591712101404645514500036458967762703694692352000/186241460818286368935039332151175406902707246982877619366601958753940483110093
    v4: 95496318082738580877112272529081300843527340761672154555777907753274572800/186241460818286368935039332151175406902707246982877619366601958753940483110093
    v5: 147966822523803895160735314580194653165712234516574564147511532345674956800/186241460818286368935039332151175406902707246982877619366601958753940483110093

vpaSolved = vpa(struct2cell(Solved))

vpaSolved = 

.

[Simon Cook]

Thanks for the vpa I was using double(). The code does run however, syms is not the ideal way to solve these linear equations. I was wondering if there is a way to rearrange to Ax=B quickly or some other method that doesnt use syms.

Answer 1

[Karim]

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]

K = 

U = [Ua; Ub]

U = 

F = [Fa;Fb]

F = 

EOM = K*U == F

EOM = 

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

Ua = 

Fb = Kab*Ua

Fb = 

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

Ua = 5×1

1.0e-03 * 0.4717 0.2520 0.3083 0.5128 0.7945

% compute the unkown forces
Fb = K(b,a) * Ua

Fb = 5×1

1.0e+03 * -0.0217 -1.1084 -0.1517 -0.8700 0.1517

% print things a bit more pretty
fprintf('Y1 = %.3f, Y2 = %.3f, X3 = %.3f, Y3 = %.3f, X5 = %.3f\n',Fb)

Y1 = -21.672, Y2 = -1108.359, X3 = -151.703, Y3 = -869.969, X5 = 151.703

fprintf('u1 = %.6f, u2 = %.6f, u4 = %.6f, v4 = %.6f, v5 = %.6f\n',Ua)

u1 = 0.000472, u2 = 0.000252, u4 = 0.000308, v4 = 0.000513, v5 = 0.000794

Notice that these results are the same as what you obtained with the symbolic toolbox.

[Simon Cook]

Thanks

Answer 2

[John D'Errico]

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)

ans = 1×2

10 10

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

ans = 

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

ans = 

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

ans = 

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)

unk_num = 10×1

1.0e+00 * 1000 1000 5.60195357711227e-13 1000 1.33427350502589e-12 3.53635636054996e-07 -1.89978583443388e-06 1.23322581330808e-06 3.23712254484621e-06 -8.54193592985598e-06

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)

sol = struct with fields:
    X3: 3466651189609543863176773926060344186418570778125/7966333403627170179660754211064126078839339744564465923719758
    X5: 4818645153557259204023352484473860068332890305125/3983166701813585089830377105532063039419669872282232961859879
    Y1: 3983166701813587077377059148336330793718643961888609764554462750/3983166701813585089830377105532063039419669872282232961859879
    Y2: 3983166701813587077377059148336330793718643961888609764554462750/3983166701813585089830377105532063039419669872282232961859879
    Y3: 7966333403627196780430882481572685299706294927378502533059764875/7966333403627170179660754211064126078839339744564465923719758
    u1: 615609230671669893905685472111789920493274937236008461642099375/2138446339850691481166836481950918929214730115071085931630174454939648
    u2: -57418684212725072514345251300944223098258114879174370279584040625/34215141437611063698669383711214702867435681841137374906082791279034368
    u4: 23784468260989039514541426270923944338654460357023725060965647375/17107570718805531849334691855607351433717840920568687453041395639517184
    v4: 6984858221910320136314736392159253994917095392194937398191940625/2138446339850691481166836481950918929214730115071085931630174454939648
    v5: -287017066430363424749745068598638581291345059494632017489045870375/34215141437611063698669383711214702867435681841137374906082791279034368

double(sol.Y1)

ans = 
                      1000

double(sol.v4)

ans = 
      3.26632382199406e-06

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)

ans = 
          10355164.0277055

You could probably resolve SOME of that by rescaling the variables, IF it was important to you.