Modeling and Simulation of Fiber Networks
In order to get our heads round the simultaneous deformation of a
network of fibers we must, I feel, introduce notation that is a
bit difficult to express in html. Hence, I ask that you go back and
forth between these pdf notes and the
demonstration codes linked below.
Our first example considers the deformation of a simple
Tent
Our next example considers the deformation of a simple
Trailer
We next animate this scene by applying an incremental force
Moving Trailer Code
The resulting animation is
Moving Trailer Trailer
Gaussian Elimination
Looking back at our fiber code we see that all of the `math' is hidden
under the innocuous backslash - namely x=S\f. We strive in this
subsection to reveal the hidden work by examining systems of equations
of increasing complexity.
The easiest systems to solve are the uncoupled, or
diagonal systems, i.e., systems of the form
Sx = f where only the diagonal elements of
S are nonzero. If S is n-by-n we
simply write
for j=1:n
x(j) = f(j)/S(j,j);
end
or, even simpler,
x = f./diag(S)
Of course this procedure breaks down if S has a zero on
its diagonal. A diagonal matrix with one or more zeros on its diagonal
is said to be singular . In such a case,
Sx=f is not solvable unless f has a
corresponding zero. Although, even this case we note that
x is left undetermined. For example
[2 0][x(1)] = [4]
[0 0][x(2)] = [0]
determines x(1) but not x(2).
After diagonal matrices the next most easily handled type is the class
of triangular matrices. These are matrices all of whose
nonzeros lie either on and above (upper) or below
(lower). In the following upper triangular system
[4 2][x(1)] = [4]
[0 3][x(2)] = [6]
we start at the bottom and notice that x(2)=6/3=2.
The first equation then reads 4x(1)+2x(2)=4 or
4x(1)=4-2x(2) or
x(1) = (4-2x(2))/4 = 0. This procedure is commonly
known as backsubstition. With respect to coding, it
requires only slightly more effort than in the diagonal case.
x = zeros(n,1);
x(n) = f(n)/S(n,n);
for j=n-1:-1:1
tmp = 0;
for k=j+1:n
tmp = tmp + S(j,k)*x(k);
end
x(j) = (f(j) - tmp)/S(j,j);
end
or
x = zeros(n,1);
x(n) = f(n)/S(n,n);
for j=n-1:-1:1
x(j) = (f(j) - S(j,j+1:n)*x(j+1:n))/S(j,j);
end
As in the diagonal case we find ourselves dividing by each diagonal
element of S. Hence, a triangular matrix with one or
more zeros on its diagonal is said to be singular .
Gaussian elimination brutally seeks to reduce a given linear system to a
triangular linear system via repeated application of two
elementary row operations.
The first is Row Swapping, for example,
tmp = S(j,:)
S(j,:) = S(k,:)
S(k,:) = tmp
and the second is Row Mixing, for example
S(j,:) = S(j,:) + magicnumber*S(k,:)
In each column we swap in order to bubble up a big
pivot then we repeatedly mix until the
pivot has eliminated every element below the diagonal.
The systematic use of these two utilities is diagrammed in the following
pseudo code:
function x = gauss(S,f)
n = length(f);
S = [S | f] Augment S with f, i.e., copy f onto column n+1 of S
Why? Because f must go along for the ride during
elimination.
for k=1:n-1 k is column in which elimination is carried out
r = row number, larger than or equal to k,
with largest value (in magnitude) in column k
if this largest value is really small then warn the user
swap row r and row k
for j=k+1:n
mix row k into row j in order to eliminate S(j,k)
end
end
if S(n,n) is really small then warn the user
strip off the changed f, i.e., copy column n+1 of S onto f
x = trisolve(S,f)
return
Here is a small example where S and f are:
2 4 2 2
1 2 4 and 3
4 2 1 1
Here is S augmented with f:
2 4 2 2
1 2 4 3
4 2 1 1
rows 1 and 3 are swapped:
4 2 1 1
1 2 4 3
2 4 2 2
our first pivot is 4 and elimination occurs in column 1:
4 2 1 1
0 3/2 15/4 11/4
0 3 3/2 3/2
rows 2 and 3 are swapped:
4 2 1 1
0 3 3/2 3/2
0 3/2 15/4 11/4
our second pivot is 3 and elimination occurs in column 2:
4 2 1 1
0 3 3/2 3/2
0 0 3 2
our third and final pivot is also 3, our
matrix is now triangular, and so we may trisolve:
so x(3) = 2/3
so 3*x(2) + (3/2)*(2/3) = 3/2 so x(2) = 1/6
so 4*x(1) + 2*(1/6) + 1*(2/3) = 1 so x(1) = 0
In the case that you encounter a small (smaller than eps) pivot
your warning message might look something like that found in my
singular diary.
Optimal Design of Fiber Nets
Biological fiber nets have (likely) evolved to best fit their niche.
How one measures `fitness' and how one identifies nature's `design variables'
are each matters of considerable debate. We will
focus on the simplest such measure, namely compliance, or work done by
the load, and suppose that nature may tune the radii of her fibers.
Recalling that work is merely the product of force and distance
we express compliance as
(1) C = xTf
where x is the solution to
(2) ATK(a)Ax = f
where we have stressed the dependence of K, the elemental
stiffness matrix, upon a, the vector of fiber cross sectional
areas. We will presume throughout that each fiber has the same
Young's modulus. Our interest then is in minimizing the C of
(1) subject to x obeying the equilibrium equation (2). A moment's
reflection will confirm that the compliance can be reduced to zero by
simply choosing infinitely fat fibers. We preclude this possibility by
constraining the volume of the net, via
(3) LTa = V
where L is the column vector of fiber lengths.
Let us solve a few small design problems by hand before opening the
big toolbox. In particular, the displacement of the loaded symmetric
(&theta=&pi/4) tent is
[x(1)] s [ a(1)+a(2) a(2)-a(1) ][f(1)]
[ ] = ------- [ ][ ]
[x(2)] a(1)a(2) [ a(2)-a(1) a(1)+a(2) ][f(2)]
where s2=1/2 as usual. And so
s(f(1)-f(2))2 s(f(1)+f(2))2
C = ------------- + --------------
a(2) a(1)
If we now invoke the volume constraint
a(1)+a(2) = sV
we find that
s(f(1)-f(2))2 s(f(1)+f(2))2
C = ------------- + --------------
sV-a(1) a(1)
and so we have reduced ourselves to a one-dimensional minimization problem.
Let us consider a pair of special loadings.
Design against falling debris: In this case f=[0 -1]' and so
s s s2V
C = --- + ------ = ------------
a(1) sV-a(1) a(1)(sV-a(1))
which attains its (positive) minimum at
a(1) = sV/2 and so a(2) = sV/2 as well.
This tells us that the two tent legs should be equally fat to best
sustain a vertical load.
Design against bears: In this case f=[1 1]' and so
4s
C = ---
a(1)
which approaches its minimum as a(1) approaches infinity.
This of course leads to the absurdity that a(2) approaches negative
infinity. If we impose the additional constraints that a(1) and
a(2) each be nonnegative then we arrive the optimal design
a(1) = sV a(2) = 0.
This corresponds to the fact that, as the load is perpendicular to the
tent's right leg, the right leg does no work. In reality, in order to
account for variations in the loading, one typically sets a nonzero
lower bound on each a(j), and for cost or construction contingencies
one likewise also imposes a finite upper bound on each a(j).
In other words, we ask that
(4) alo &le a(j) &le ahi for each j
We now arrive at the General Case
min xTf subject to ATK(a)Ax = f
a LTa = V
alo &le a(j) &le ahi for each j.
Matlab has a large toolbox for solving such optimization
problems, type help optim.
The routine most pertinent to our problem is fmincon.
Its application to our simple tent is coded here.
Running this code with f=[.5;-1] reveals the figure above.
