[Date Index] [Thread Index] [Author Index]

tensegrity and FindRoot for quadratic systems

• To: mathgroup at smc.vnet.net
• Subject: [mg100249] tensegrity and FindRoot for quadratic systems
• From: alessio <alessio at kth.se>
• Date: Fri, 29 May 2009 20:58:00 -0400 (EDT)

I have implemented the following code for simulating a 3x3 lattice of vertices, connected by 16 rigid bars to form 8 triangles, fixed at two points, and under the action of the gravity (the configuration can be seen in the animation that the code creates at the end). My problem is that for this particular geometry, FindRoot is not able to find a solution within the prescribed tolerance. Does this mean that I am in a local minimum? What can I do to avoid this situation?



len[p1_, p2_] := 
 
 (p2[[1]] - p1[[1]])*(p2[[1]] - p1[[1]]) +
  (p2[[2]] - p1[[2]])*(p2[[2]] - p1[[2]]) +
  (p2[[3]] - p1[[3]])*(p2[[3]] - p1[[3]])

strain[q_] := {
  
  len[q[[1]], q[[2]]] - 25,
  len[q[[2]], q[[4]]] - 25,
  len[q[[3]], q[[4]]] - 25,
  len[q[[1]], q[[3]]] - 25,
  len[q[[1]], q[[4]]] - 50,
  len[q[[2]], q[[5]]] - 25,
  len[q[[5]], q[[6]]] - 25,
  len[q[[4]], q[[6]]] - 25,
  len[q[[2]], q[[6]]] - 50,
  len[q[[6]], q[[7]]] - 25,
  len[q[[4]], q[[7]]] - 50,
  len[q[[7]], q[[8]]] - 25,
  len[q[[4]], q[[8]]] - 25,
  len[q[[3]], q[[8]]] - 50,
  len[q[[3]], q[[9]]] - 25,
  len[q[[8]], q[[9]]] - 25,
  q[[1, 1]] - 0,
  q[[1, 2]] - 0,
  q[[1, 3]] - 0,
  q[[9, 1]] - 0,
  q[[9, 2]] - 10,
  q[[9, 3]] - 0
  }

undeformed = {{0, 0, 0}, {5, 0, 0}, {0, 5, 0}, {5, 5, 0}, {10, 0, 
    0}, {10, 5, 0}, {10, 10, 0}, {5, 10, 0}, {0, 10, 0}};

force = {
   {0, 0, 98},
   {0, 0, 98},
   {0, 0, 98},
   {0, 0, 98},
   {0, 0, 98},
   {0, 0, 98},
   {0, 0, 98},
   {0, 0, 98},
   {0, 0, 98}
   };

EL[u_, p_, q_, \[Lambda]_] := Flatten[{
   
   {(q - p)/dt - (p - u)/dt} - \[Lambda] .D[
      strain[q], {q}] - {force},
   {strain[q]}
   }]

deformed = {
   {x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3},
   {x4, y4, z4}, {x5, y5, z5}, {x6, y6, z6},
   {x7, y7, z7}, {x8, y8, z8}, {x9, y9, z9}
   };
dt = 1/1000;

multipliers = {
   {l1, l2, l3, l4, l5, l6, l7, l8, l9, l10, l11, l12, l13, l14, l15, 
    l16, l17, l18, l19, l20, l21, l22}
   };

hess[u_, p_, sym_, lag_] := D[
  EL[u, p, sym, lag],
  {Flatten[{sym, lag}]}
  ]

hess[undeformed, undeformed, deformed, multipliers];

step[u_, p_] := FindRoot[
  
  EL[u, p , {
    {x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3},
    {x4, y4, z4}, {x5, y5, z5}, {x6, y6, z6},
    {x7, y7, z7}, {x8, y8, z8}, {x9, y9, z9}
    }, {
    {l1, l2, l3, l4, l5, l6, l7, l8, l9, l10, l11, l12, l13, l14, l15,
      l16, l17, l18, l19, l20, l21, l22}
    }],
  
  {{x1, p[[1, 1]]}, {y1, p[[1, 2]]}, {z1, p[[1, 3]]},
   {x2, p[[2, 1]]}, {y2, p[[2, 2]]}, {z2, p[[2, 3]]},
   {x3, p[[3, 1]]}, {y3, p[[3, 2]]}, {z3, p[[3, 3]]},
   {x4, p[[4, 1]]}, {y4, p[[4, 2]]}, {z4, p[[4, 3]]},
   {x5, p[[5, 1]]}, {y5, p[[5, 2]]}, {z5, p[[5, 3]]},
   {x6, p[[6, 1]]}, {y6, p[[6, 2]]}, {z6, p[[6, 3]]},
   {x7, p[[7, 1]]}, {y7, p[[7, 2]]}, {z7, p[[7, 3]]},
   {x8, p[[8, 1]]}, {y8, p[[8, 2]]}, {z8, p[[8, 3]]},
   {x9, p[[9, 1]]}, {y9, p[[9, 2]]}, {z9, p[[9, 3]]},
   {l1, 0}, {l2, 0}, {l3, 0}, {l4, 0},
   {l5, 0}, {l6, 0}, {l7, 0}, {l8, 0},
   {l9, 0}, {l10, 0}, {l11, 0}, {l12, 0},
   {l13, 0}, {l14, 0}, {l15, 0},
   {l16, 0}, {l17, 0}, {l18, 0},
   {l19, 0}, {l20, 0}, {l21, 0}, {l22, 0}}
  
  , Jacobian -> hess[u, p, {
     {x1, y1, z1},
     {x2, y2, z2},
     {x3, y3, z3},
     {x4, y4, z4},
     {x5, y5, z5},
     {x6, y6, z6},
     {x7, y7, z7},
     {x8, y8, z8},
     {x9, y9, z9}
     }, {
     {l1, l2, l3, l4, l5, l6, l7, l8, l9, l10, l11, l12, l13, l14, 
      l15, l16, l17, l18, l19, l20, l21, l22}
     }]
]

tri = {
   {0, 0, 0}, {x2, y2, z2}, {x4, y4, z4},
   {x3, y3, z3}, {0, 0, 0}, {x4, y4, z4},
   {x6, y6, z6}, {x2, y2, z2}, {x5, y5, z5},
   {x6, y6, z6}, {x7, y7, z7}, {x4, y4, z4},
   {x8, y8, z8}, {x3, y3, z3}, {x9, y9, z9},
   {x8, y8, z8}, {x7, y7, z7}
   };

edgeT = {tri} /. {
    x2 -> 5, y2 -> 0, z2 -> 0,
    x3 -> 0, y3 -> 5, z3 -> 0,
    x4 -> 5, y4 -> 5, z4 -> 0,
    x6 -> 10, y6 -> 5, z6 -> 0,
    x5 -> 10, y5 -> 0, z5 -> 0,
    x7 -> 10, y7 -> 10, z7 -> 0,
    x8 -> 5, y8 -> 10, z8 -> 0,
    x9 -> 0, y9 -> 10, z9 -> 0
    };

vtx = {
   {x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3},
   {x4, y4, z4}, {x5, y5, z5}, {x6, y6, z6},
   {x7, y7, z7}, {x8, y8, z8}, {x9, y9, z9}
   };

vold = undeformed
vnew = undeformed;

For[i = 0, i < 10, i++,
 
 sol = step[vold, vnew];
 vold = vnew;
 vnew = vtx /. sol;
 edge = tri /. sol;
 
 Print[i];
 AppendTo[edgeT, edge];
 ]

Animate[
 Graphics3D[Line[edgeT[[i]]],
  PlotRange -> {{-2, 12}, {-2, 12}, {-1, 10}}],
 {i, 1, 10},
 AnimationRunning -> False
 ]