Least squares method
The Collocation Method is quite simple, but it does not restrict the magnitude of £, £v, £n> and Ер outside the collocation points. A natural extension of this method is to equate to zero the sum of the residuals within sub-domains. However, this approach may conceal quite large residuals and lead to computational difficulties (Linz, 1979). The Least Squares Method overcomes that inconvenience by minimizing the sum of squared errors.
Consider first how to find an approximation for H(r, t) if some of its values are known at few points. In this case, the number and the type of base functions hj(r, t) of the approximant h(r, t) are previously elected (according to the expected physical behavior to be described and to the desired accuracy). Then, the sum of the squared
Figure 3.4 Interpolated hydraulic head map from a I00x |00m equally spaced grid nodes of calculated head values.
residuals e2 only depends on arbitrary variations of its numerical coefficients aj. It is sufficient to find which set of coefficients aj minimizes this sum within the flow domain and its boundaries. This set is the solution of the simultaneous equations below' (see Section 3.6.1).
In the last constraint, the coefficients denote the direction cosines of the outward-normal n to the boundary element dS. It must be emphasized that as the analytic expression for H(r, t) is unknown, the integrals /H(r, t) must be found numerically from the known values of H(r, t) at monitoring points.
The next example clarifies this method for a simple case.
The following graph (see Figure 3.5) shows the time-variation of the zinc content, from September 30, 2004 to March 31, 2005, found in groundwater samples collected each week (but not at the same weekday) in a point located 259 m downstream from an industrial process.
Vales of dissolved Zn content in groundwater samples, as a function of the elapsed time t, can be approximated by a cubic polynomial fitted to monitored data by the least square method, but considering the irregular time interval between successive samples (varying from 1 to 7 days). Denoting by Pj its unknown numerical coefficients, a general polynomial expression for j varying from 0 to Np can be written as (is this example Np=4):
Minimizing the sum of the squared residuals e2 over the total timespan, considering a cubic approximant, yields four simultaneous equations:
Figure 3.5 Zinc content time-variation monitored in a control point. Faulty sampling, defective analysis, or rain infiltration dilution may explain the disproportionate data spread.
In the above equations, T denotes the time elapsed since sampling began. The following expression in these equations corresponds to trivial numerical integrations multiplied by the weighting factors t°, t1, t2, and t3, but considering the irregular time interval 5t between consecutive Zn values.
From a structural point of view, these four simultaneous equations have a precise meaning. The first equation expresses the equivalence of the approximant areas and by the known solution (in this case, defined by discrete values, despite its lousy accuracy). The second equation forces the centers of gravity of these areas to be equally distant from the ordinate axis (equivalence or first moments). The third one requires equivalence of inertia moments (second moments) and the fourth, of third moments.
Solving these simultaneous equations gives the unknown numerical coefficients of the cubic approximant:
Figure 3.6 shows the resulting approximant.
Figure 3.6 Cubic polynomial fitted to the scattered data, taking into account the irregularvariationofthetimeintervalStbetweensuccessiveZn contents values.
The preceding simultaneous equations can be written using the symbolic continuous inner product notation, and they may also be combined and pondered by weighting coefficients wv and ws:
These minimum constraints within the flow domain or the minimum constraints at the flow domain boundaries are exactly satisfied by
the above equations. How well they approximate each location’s solution depends on the relative magnitude of the weighting coefficients wv and ws.
The Least Squares Method also holds when the differential operators Dy and DN are applied to the residuals [h(r, t)-H(r, t)]. In this case, the unknown coefficients aj after minimizing the sum of the squared residuals £v2 and eN2 corresponds to the solution of the following simultaneous equations:
Both requirements can also be combined (occasionally using weighting coefficients). Then, using the convenient inner product notation, these conditions can be written as:
The next example clarifies this method.
Parallel linear drains buried under an embankment must be designed to keep the rising water table caused by the infiltration of 75% of a continuous and extremely exceptional 100mm rain for four days, below a predetermined maximum level of 0.5 m, corresponding to an average infiltration rate of 18.75 mm/day.
Being T=4 days, a quadratic polynomial H(x, t, T) may approximate the rise of the water as a product Fi(x) F2(t). In this case, x is the distance measured from the drains centerline and t is the elapsed time since the start of the continuous rain:
In the above formula, HL is the vertical distance from the embankment impervious base to the drain base level. Obviously, the units of the numerical coefficients must be compatible, that is, ao: [L], ap [-], a?. [L_1], HL: [L] and x: [L],
Flow symmetry implies zero horizontal discharge at the drain’s centerlines, that is, qn=0, at x = 0, corresponding to the following Neumann condition at x = 0:
Substituting H for the selected approximation:
Taking the derivate:
From which a( = 0.
The draining capacity is designed to assure zero hydraulic pressure at the drain centers. Then, the following Dirichlet condition at x = L is:
From which a 2 = -ao/L2.
By replacing a and a2 for their values, the approximant reduces to:
Now, the Least Squares Method can be employed to find a0. The unconfined Dupuit’s continuity equation in (x, t) and the Neumann condition at (L, t) for a 1 m wide strip are, respectively:
To effectively drain qn must balance all water accretion coL, that is, qn = coL.
Now, replacing the last expression H(x, t, T) in the Dupuit and Neumann conditions give:
Then, minimizing the squared residuals within the flow domain and the drain boundaries, making hj = 1 and combining these two requisites with equal weighting coefficients:
Evaluating these two integrals gives:
Solving for ao:
Finally, substituting ao in the last expression for H(x, t, T): Figure 3.7 shows results for progressive time-periods.