kill(all);

  y▲
   │
n10├─────┐n9           ┬
   │   ▲ │             │ a
   │   │ │             │ 
   ├───┼─┤             ┼
   │   │ │             │ a
   │  G└─┼────►        │
   ├─────┼─────┐n4     ┼
   │     │     │       │ a
   │     │     │       │
  O└─────┴─────┴───►x  ┴
   n1          n3  
      a     a
   ├─────┼─────┤

assume(a>0,E>0);

integration of a generic function across the area

aint(f):=
     integrate(integrate( f(x,y) ,x,0*a,2*a),y,0*a,1*a)
    +integrate(integrate( f(x,y) ,x,0*a,1*a),y,1*a,3*a)$

integrate the [x,y] -> 1 function

area: aint( lambda([x,y],1) );

integrate the 
[x,y] -> x-xG and
[x,y] -> y-yG functions
to obtain the first order moments of area

SxG: aint( lambda([x,y],x-xG) );

SyG: aint( lambda([x,y],y-yG) );

find the axial center of elasticity coordinates

linsolve([SxG=0,SyG=0],[xG,yG]),globalsolve=true;

integrate the
[x,y] -> (y-yG)*(y-yG) ,
[x,y] -> (x-xG)*(y-yG) , and
[x,y] -> (x-xG)*(x-xG) functions
to obtain the second order moments of area

Jxx : aint( lambda([x,y],(y-yG)*(y-yG) ) );

Jxy : aint( lambda([x,y],(y-yG)*(x-xG) ) );

Jyy : aint( lambda([x,y],(x-xG)*(x-xG) ) );

beam bending stiffness matrix

EJmat : matrix([E*Jyy,E*Jxy],[E*Jxy,E*Jxx])/(E*Jxx*E*Jyy-(E*Jxy)^2);

its eigendecomposition to obtain principal moments, and the associated directions

[[vals,vmult],vecs]:eigenvectors(EJmat);

linsolve(E*J11=1/vals[1],J11);

linsolve(E*J22=1/vals[2],J22);

vecs[1][1]$
[%[1],%[2]]/sqrt(%[1]^2+%[2]^2)$
[v1x,v1y]:%;

vecs[2][1]$
[%[1],%[2]]/sqrt(%[1]^2+%[2]^2)$
[v2x,v2y]:%;

right angle moving from vxx to v11

alpha_rad:atan2(v1y,v1x),numer;
alpha_grad: alpha_rad/%pi*180,numer;

Mixed moment of inertia

J12 : aint( 
    lambda(
        [x,y],
         (v1x*(x-xG)+v1y*(y-yG)) /* x1-x1G */
        *(v2x*(x-xG)+v2y*(y-yG)) /* x2-x2G */
    ) 
);

[J11,J22] : aint( 
    lambda(
        [x,y],
        [
            (v2x*(x-xG)+v2y*(y-yG))^2,/* (x2-x2G)^2 */
            (v1x*(x-xG)+v1y*(y-yG))^2 /* (x1-x1G)^2 */
        ]
    ) 
);

nodes defining the convex envelope of the cross section, as in figure

ni:[1,3,4,9,10]$
nx:[0,2,2,1, 0]*a-xG$
ny:[0,0,1,3, 3]*a-yG$
N:length(ni)$

a few test cases, for FE comparison

[Mx,My]:[1,0]*a*(E*a^2)$
[cx,cy]:[Mx*E*Jyy+My*E*Jxy,Mx*E*Jxy+My*E*Jxx]/(E*Jxx*E*Jyy-(E*Jxy)^2);
ne:-cy*nx+cx*ny;
[lmax(ne),lmin(ne)]$
%,numer;

[Mx,My]:[0,1]*a*(E*a^2)$
[cx,cy]:[Mx*E*Jyy+My*E*Jxy,Mx*E*Jxy+My*E*Jxx]/(E*Jxx*E*Jyy-(E*Jxy)^2);
ne:-cy*nx+cx*ny;
[lmax(ne),lmin(ne)]$
%,numer;

[Mx,My]:[v1x,v1y]*a*(E*a^2)$
[cx,cy]:[Mx*E*Jyy+My*E*Jxy,Mx*E*Jxy+My*E*Jxx]/(E*Jxx*E*Jyy-(E*Jxy)^2);
ne:-cy*nx+cx*ny;
[lmax(ne),lmin(ne)]$
%,numer;

[Mx,My]:[v2x,v2y]*a*(E*a^2)$
[cx,cy]:[Mx*E*Jyy+My*E*Jxy,Mx*E*Jxy+My*E*Jxx]/(E*Jxx*E*Jyy-(E*Jxy)^2);
ne:-cy*nx+cx*ny;
[lmax(ne),lmin(ne)]$
%,numer;