The statement of the problem (
) can now be refined to
incorporate the approximations of v and h over the mesh. The
problem of formulating a numerical solution to the continuity
equation of ice is to find 
, where
such that 
on 
and the following
equation is satisfied.
In addition 
is given by
where for all 
, 
must be zero on
The system of equations is arrived at by substituting equation
() and () into the variational statement
equation () and simplifying. Using the fact that 
is
an arbitrary trial function so that the values of the constants,
, are arbitrary, N algebraic equations are obtained.
The N equations are obtained by systematically
selecting values of 
such that only one constant is non-zero
for each of the N equations. These N equations
can be expressed in a matrix form as
Following traditional engineering nomenclature, the matrix 
is called the stiffness matrix, the matrix 
is called
the capacitance matrix, and the vector 
is called the load vector.
Components are given by the following expressions.
There are some special features of this matrix equation which
should be noted. First, both the stiffness matrix K and the
capacitance matrix C are sparse. This is due to the fact that
the global basis functions are nonzero only when nodes i and j
belong to the same element. Therefore, the entry 
(or 
)
will be nonzero only on elements containing both i and j.
Secondly, K and C are symmetric matrices. These facts are important
in that major economies are realized in both storage and solution speed for the
symmetric, sparse set of equations.
An important consequence of the form of the basis set chosen
is that the integrals in equations (), (),
and () can be calculated by adding contributions from
integrals over each individual element. This can be shown
by considering an element of the mesh, 
.
Substituting interpolants of v and h over the element of the form
where 
are local shape functions for 
and
is the number of nodes in 
, the linear system for
the element 
is obtained.
where
Here, 
, 
, and 
are the element
stiffness matrix, element capacitance matrix, and element load vector
which will be combined to form the global stiffness matrix, global
capacitance matrix, and global load vector. The element flux vector,
, represents the actual flux 
across
. The 
's are those pieces of the basis
functions, 
's, which are non-zero over the element 
.
Like the basis functions, the shape functions are 1 at exactly one node,
and zero at all the rest. Hence, for a four-node element, there are four
shape functions whose products must be integrated for equations
(), (), (), and ().
The components of the global system of equation () are
obtained by
summing equations (), (), (), and
() over all the elements of the mesh. Thus, the matrix and
vector entries of equation () are given by
where E is the number of elements.
The entries of the element stiffness matrix (equation
()), the element capacitance matrix (equation
()), and the element load vector (equation ())
are computed using Gaussian quadrature of order 3. The boundary
integral for the natural boundary conditions (equation ())
is computed in the same way and this value is added into the appropriate
element load vector.
In the steady state case equation() reduces to the
matrix equation
For the time dependent case, an unconditionally stable backward
difference scheme is used to approximate 
. That
is,
where T indicates the time step. Substituting equation ()
into equation () and simplifying yields
Both the steady state equation () and the time dependent
equation () are simple matrix equations which can be solved
by traditional matrix methods.
