input-matrix4

有限元程序

/* ================================================================ */
/* Test linear equation solvers and eigenvalue/eigenvector routines */
/* ================================================================ */

print "==============================================\n";
print "*** Linear matrix equation solver [A]{x} = {b}\n";
print "==============================================\n\n";

A = [3,  4, 5, 7,  9,  10;
     1,  6, 7, 9, 20, 100;
     11,26,47, 9,  2,  13;
     12, 6,87, 1,  3,  33;
     11, 0, 1, 6,  31,  5;
     81, 71,2, 2,   1,  54];

A = ColumnUnits(A, [N/m, N/m, N, N/m, N/m, N/rad]);
A = RowUnits(A, [m], [3]);
A = RowUnits(A, [m], [6]);

b = [20 N; 100 N;  10 N*m;
      3 N;   2 N; 200 N*m];
PrintMatrix(A, b);

print "\n";
print "*** Solve by matrix inverse (then check results)\n";

Y = Inverse(A);
X = Y*b;
Residual = A*X-b;
PrintMatrix(X, Residual);

print "\n";
print "*** Solve by matrix solver (then check results)\n";

X = Solve(A, b);
Residual = A*X-b;
PrintMatrix(X, Residual);

print "\n";
print "*** Solve by Matrix LU Decomposition and Backsubstitution (then check results)\n";

LU = Decompose(A);
X  = Substitution(LU, b);
Residual = A*X-b;
PrintMatrix(X, Residual );

print "\n";
print "======================================\n";
print "*** Solve matrix equations (Example 2)\n";
print "======================================\n\n";

/* Structure calculation */

E    = 200 MPa;
I    = 0.4 m^4;
Area = 0.1 m^2;
L    = 5.0 m;

/* Setup and print 3x3 stiffness matrix */

stiff = Matrix([3,3]);
stiff = ColumnUnits( stiff, [N/m, N/m, N/rad]);
stiff = RowUnits( stiff,     [m], [3] );

stiff[1][1] =  E*Area/L;
stiff[2][2] = 12*E*I/L^3;
stiff[2][3] = -6*E*I/L/L;
stiff[3][2] = -6*E*I/L/L;
stiff[3][3] =  4*E*I/L;

PrintMatrix(stiff);

/* Setup and print 3x1 force vector  */

force = [ 0.0 N ; -5 N; 0.0 N*m ];
PrintMatrix(force);

/* Compute and print displacement matrix */

displ = Solve(stiff, force );
PrintMatrix(displ);
PrintMatrixCast(displ, [ cm; deg ] );

print "\n";
print "===============\n";
print "*** End of Test\n";
print "===============\n";

quit;