Least Squares Regression

Input := 

A = 68.7672;
c = 3.0022;
L = 299.2239;
fc = 625.0925;
f = fc - A(Cosh[c x/L]- 1);

The parabola which best models the Gateway Arch will have symmetry with respect to the y-axis. Thus, it will be a quadratic function with linear coefficient zero. A general quadratic of this type can be expressed as follows, where p1 and p2 are parameters.

Input := 

g = p1 x^2 + p2;

One way of finding the parabolic curve which best fits the catenary is to measure the difference between the curves; i.e., consider f(x) - g(x). However, if the curves intersect, then on some intervals this difference will be positive and on some it will be negative. Even though the differences may be large in magnitude, the opposite signs allow for the possibility that the differences will in effect "cancel".

One way to force the difference to be positive is to express the difference as f(x) - g(x) on subintervals where f(x) > g(x) and as g(x) - f(x) on subintervals where g(x)> f(x). This requires finding the point(s) of intersection of the curves (which is not necessarily an easy task.) Another way to eliminate the difficulty of the negative differences is to use the absolute value of the difference f(x) - g(x), but this causes complications in the subsequent integration calculations if using current computer algebra systems.

These problems are eliminated by squaring the difference. The expression
( f(x) - g(x) )^2
will give a (positive) measure of the square difference between the two curves at any point. The total difference is given by integrating over the length of the curves. Due to the symmetry in our problem, it is necessary only to integrate from 0 to the right end point. The parabola will best fit the catenary when total square error is minimized. Since we are performing what is know as a least squares regression, we refer to the total square error as SE (for Square Error), and we remark that this quantity depends on the parameters p1 and p2 for the candidate fitting curve.

Input := 

SE = Integrate[(f - g)^2, {x, 0, L}]

Output =

          7             9
7.31452 10  - 5.58561 10  p1 + 
 
            11   2
  4.79746 10   p1  - 
 
  277612. p2 + 
 
            7                   2
  1.78607 10  p1 p2 + 299.224 p2

To minimize the square error, we find its critical points by finding the zeroes of its derivatives.

Input := 

derSEp1 = D[SE, p1];
derSEp2 = D[SE, p2];

By setting each partial derivative equal to zero and solving the resulting system of equations simultaneously, we find the values of p1 and p2 that minimize the sum of the square differences.

Input := 

sol = Solve[{derSEp1 == 0,derSEp2 == 0}, {p1, p2}]

Output =

{{p1 -> -0.00633079, 
 
   p2 -> 652.829}}

How does the least squares parabola compare to the actual arch? We will retrieve the values for p1 and p2 from the above solution and then plot each curve to see. The true arch is the solid curve and the approximating curve is dashed.

Input := 

p1star = p1 /. sol[[1,1]];
p2star = p2 /. sol[[1,2]];
bestg = g /. {p1 -> p1star, p2 -> p2star}

Output =

                      2
652.829 - 0.00633079 x

Input := 

Plot[{f, bestg}, {x, -L, L},
 AspectRatio -> Automatic,
 PlotStyle -> {{AbsoluteThickness[2]},
               {AbsoluteThickness[3],
                AbsoluteDashing[{3, 6}]}},
 AxesLabel -> {"[ft]", "height [ft]"}];

Its important to note that the fitting curve is only best fitting over the length of the true arch. It's a bit more complicated to get a best fitting curve all the way to the ground in the least squares sense. We offer a solution to that case in ISSUES IN SOLUTIONS.