function [po,ERROR] = ortpol ( X , F , W , NTERMS )
%*****************************************************************************80
%
%MATLAB version of SUBROUTINE ORTPOL on page 263f of conte/deboor
%************** I N P U T ************
% (X(i), F(i)), i=1:length(X) are the abscissae and ordinates of
% the given data points to be fitted
% W is the vector of weights to be used in the weighted scalar product
% NTERMS-1 is the upper bound on the degree of the polynomial approximant p
%************** O U T P U T ************
% po is a structure containing the coefficients D of the polynomial p
% with respect to an orthonormal polynomial basis, as well as the
% arrays B and C for the three-term recurrence that generates that
% orthonormal basis.
% The call ORTVAL(po,x) provides the value of p at x.
% ERROR(i) is the difference between Y(i) and p(X(i)).
%************** M E T H O D ************
% The sequence P0, P1, ..., PNTERMSm1 of orthogonal polynomials with
% respect to the discrete inner product
% (P,Q) := SUM ( P(X(I))*Q(X(I))*W(I), I=1:length(X) )
% is generated in terms of their three-term recurrence
% PJp1(x) = (x - B(J+1))*PJ(x) - C(J+1)*PJm1(x)
% and the coefficient D(J) of the weighted least squares approximation
% to the given data is obtained correspondingly as
% D(J+1) = (F,PJ)/(PJ,PJ), J=0:NTERMS-1.
% Actually, in order to reduce cancellation, (F,PJ) is computed as
% (ERROR,PJ), with ERROR = F initially and, for each J , ERROR is reduced
% by D(J+1)*P(J) as soon as D(J+1) becomes available.
%
% Reference:
%
% Samuel Conte, Carl de Boor,
% Elementary Numerical Analysis,
% Third Edition,
% SIAM, 2017,
% ISBN: 978-1-611975-19-2.
%
W = W(:); X = X(:).'; F = F(:).'; B = X*W; C = 0; D = F*W; S = sum(W);
D = D/S; ERROR = F - D;
if (NTERMS==1), po = struct('B',B,'C',C,'D',D); return, end
B = [B/S,zeros(1,NTERMS-1)]; PJM1 = ones(size(X)); PJ = X - B(1);
C = [C,zeros(1,NTERMS-1)];
D = [D,zeros(1,NTERMS-1)];
S = [S,zeros(1,NTERMS-1)];
for J=2:NTERMS
P = PJ.*(W.'); D(J) = ERROR*(P.');
P = P.*PJ; B(J) = X*(P.'); S(J) = sum(P);
D(J) = D(J)/S(J);
ERROR = ERROR - D(J)*PJ;
if J==NTERMS, continue, end
B(J) = B(J)/S(J); C(J) = S(J)/S(J-1);
P = PJ; PJ = (X-B(J)).*PJ - C(J)*PJM1; PJM1 = P;
end
po = struct('B',B,'C',C,'D',D);
return
end