Stress example

> with(linalg):

The framework is a rectangle with an interior vertex joined to the other four.

> v1:=[1,1]; v2:=[-2,1]; v3:=[-2,-2]; v4:=[1,-2]; v5:=[0,0];

Compute the stress matrix.

> A:=matrix([[op(v2-v1),op(v1-v2),0,0,0,0,0,0],[0,0,op(v3-v2),op(v2-v3),0,0,0,0], [0,0,0,0,op(v4-v3),op(v3-v4),0,0], [op(v4-v1),0,0,0,0,op(v1-v4),0,0], [op(v5-v1),0,0,0,0,0,0,op(v1-v5)], [0,0,op(v5-v2),0,0,0,0,op(v2-v5)], [0,0,0,0,op(v5-v3),0,0,op(v3-v5)], [0,0,0,0,0,0,op(v5-v4),op(v4-v5)]]);

> rank(A);

Since the matrix has 10 columns and rank 7, the dimension of the nullspace is 3, which equals (3 choose 2). So the framework is infinitesimally rigid.

> B:=nullspace(transpose(A));

The nullspace is one-dimensional, so there is a non-trivial stress.

> multiply(B[1],A);

Construct the stress polynomial and verify its property.

> l12:=B[1][1]; l23:=B[1][2]; l34:=B[1][3]; l14:=B[1][4]; l15:=B[1][5]; l25:=B[1][6]; l35:=B[1][7]; l45:=B[1][8];

> l11:=-l12-l14-l15; l22:=-l12-l23-l25; l33:=-l23-l34-l35; l44:=-l14-l34-l45; l55:=-l15-l25-l35-l45;

> b:=l12*x1*x2+l23*x2*x3+l34*x3*x4+l14*x1*x4+l15*x1*x5+l25*x2*x5+l35*x3*x5+l45*x4*x5+ l11/2*x1^2+l22/2*x2^2+l33/2*x3^2+l44/2*x4^2+l55/2*x5^2;

> factor(b);

> M:=transpose(matrix([[op(v1),1], [op(v2),1], [op(v3),1], [op(v4),1], [op(v5),1]]));

> g:=grad(b, vector([x1,x2,x3,x4,x5]));

> multiply(M,g);

>