% LEAST_SQUARES_SOB Polynomial least squares approximation in Sobolev space. % % This routine generates the (n+1)x(N*(s+1)) array phat of % Sobolev least squares approximants of degrees <=n, and % their derivatives of orders <=s, evaluated at the N % abscissae t_k of the discrete Sobolev inner product, and % the (n+1)-vector c of Sobolev_Fourier coefficients. Here, % s is the highest-order derivative appearing in the % Sobolev inner product and is determined by the routine % automatically from the size of the input array xw. The N % values of the derivative of order sig of the (nu)th- % degree approximant, nu<=n, are output in positions % (nu+1,sig+1:s+1:N*(s+1)) of the array phat. The Nx(s+1) % array of the input array f contains the N values of the % given function and its first s derivatives at the N % points t_k. The abscissae t_k and the weights w_k^{(sig)} % of the Sobolev inner product are input via the Nx(s+2) % array xw containing the abscissae in the first column % and the successive weights for the derivatives in the % subsequent s+1 columns. The user also has to provide % the NxN upper triangular array B of the recurrence % coefficients for the Sobolev orthogonal polynomials, % which for s=1 can be generated by the routine % chebyshev_sob.m, and for arbitrary s by the routine % stieltjes_sob.m. % function [phat,c]=least_squares_sob(n,f,xw,B) N=size(xw,1); s=size(xw,2)-2; Ns=N*(s+1); p=zeros(n+1,Ns); p2=zeros(1,n+1); % % Generate the matrix of Sobolev orthogonal polynomials % and their derivatives, along with the array of their norms % sN=1:s+1:Ns; p(1,sN)=1; p2(1)=sum(xw(:,2)); t=xw(:,1)'; for k=1:n for sig=1:s+1 sigN=sig:s+1:Ns; bsum=zeros(1,N); for j=1:k bsum=bsum+B(j,k)*p(k-j+1,sigN); end if sig==1 p(k+1,sigN)=t.*p(k,sigN)-bsum; else p(k+1,sigN)=t.*p(k,sigN)+(sig-1).*p(k,sigN-1)-bsum; end p2(k+1)=p2(k+1)+sum(xw(:,1+sig)'.*(p(k+1,sigN).^2)); end end % % Compute the matrix of least squares approximants and their % derivatives, along with the array of Fourier coefficients % c=zeros(n+1,1); phat=zeros(n+1,Ns); e=f; for k=1:n+1 for sig=1:s+1 sigN=sig:s+1:Ns; c(k)=c(k)+sum(xw(:,1+sig)'.*e(:,sig)'.*p(k,sigN)); end c(k)=c(k)/p2(k); if k==1 for sig=1:s+1 sigN=sig:s+1:Ns; phat(1,sigN)=c(1)*p(1,sigN); end else for sig=1:s+1 sigN=sig:s+1:Ns; phat(k,sigN)=phat(k-1,sigN)+c(k)*p(k,sigN); end end if k==n+1, return, end for sig=1:s+1 sigN=sig:s+1:Ns; e(:,sig)=e(:,sig)-c(k)*p(k,sigN)'; end end