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

26 vues (au cours des 30 derniers jours)
Simon Cook le 29 Déc 2022
Modifié(e) : John D'Errico le 29 Déc 2022
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.
5 commentairesAfficher 3 commentaires plus anciensMasquer 3 commentaires plus anciens
Star Strider le 29 Déc 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])
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 le 29 Déc 2022
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.

Connectez-vous pour commenter.

Réponse acceptée

Karim le 29 Déc 2022
Modifié(e) : Karim le 29 Déc 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]
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.
1 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Simon Cook le 29 Déc 2022
Thanks

Connectez-vous pour commenter.

Plus de réponses (1)

John D'Errico le 29 Déc 2022
Modifié(e) : John D'Errico le 29 Déc 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)
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.
0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

Catégories

En savoir plus sur Equation Solving dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by