1function [wts,LebesgueConstant] = get_lmmwts(xpts, npts, nqloc,rlen,dt,delta,winData,computeLebesgueCst)
2% GET_LMMWTS gives the weights for the LMM-style treatment of [0,dt]
3% integrals or [0,L] integrals where L<dt
4%
5% wts = get_lmmwts(xpts, npts, nqloc,rlen,phit,blendingTol,dt)
6% performs barycentric interpolation using 'xpts' (a vector of size 'npts')
7% as interpolation nodes, then applies GL quadrature with 'nqloc'
8% quadrature nodes to evaluate the integral over
9% [xpts(npts-1),xpts(npts-1) + rlen].
10% The weights are modulated by the blending function needed for the windowing
11% in the analytic split of the solution. Here 'phi_wts' are cheb weights
12% to compute the blending function.
13% 'chebApproxInfo': a struct with fields domain [ax,bx] over which cheb
14% weights are computed; ninters, number of subintervals, and nord,
15% order of polynomial interpolation or number of cheb nodes.
16% 'computeLebesgueConstant' is a button to compute Lebesgue constant if
17% needed (set to one).
18%
19% Note: function rewritten from Fortran function provided by Leslie
20% Greengard
22if nargin==0, test_getlmmwts; return; end
24rintmatloc = zeros(nqloc, npts); % Allocate memory for integration matrix
26% Obtain Gauss-Legendre nodes and weights -> on [-1,1] interval
27gl = glwt_prep(nqloc);
28glnodest = gl.x0;
29glweightst = gl.w0;
31% precompute the barycentric weights
32X=repmat(xpts,1,npts);
33w = 1./prod(-X+X.'+eye(npts),1);
35% Loop over intervals (xpts(1),..,xpts(npts))
36aa = xpts(npts-1); % here the interval of integration [aa,bb] = [xpts(npts-1),xpt(npts)]
38% GL nodes adjusted to interval
39xt = aa + rlen*(glnodest + 1.0)/2.0;
40phi = generalwindow((dt - xt)/delta,winData);
41phimult = (1 - phi);
43glwloc = rlen*glweightst.*phimult;
45tol = 1.0e-20;
46for iquad = 1:nqloc % loop over each GL node
47 rr = prod((xt(iquad) - xpts));
48 rintmatloc(iquad, :) = rr*w./(xt(iquad) - xpts');
49 rintmatloc(iquad,(abs(xt(iquad) - xpts) < tol)) = 1;
50end
52% Compute the Lebesgue constant
53if(computeLebesgueCst == 1)
54 LebesgueConstant = norm(rintmatloc,inf);
55else
56 LebesgueConstant = NaN;
57end
59% Compute weights: GL quadrature applied to each lj
60wts = glwloc'*rintmatloc;
62end % end of get_lmmwts function
64function test_getlmmwts
65clf;
66ax = 0; bx = 1;
67h = 1e-4;
68xpts = ax:h:bx;
69rlen = h;
70tol = 1e-12; % error tolerance
71gam = .5; % pad gam*Nyquist for window
72theta = log(1/tol);
73W = ceil(2*theta/(pi*gam));
74[phit1,phitt1] = window(theta,W);
75chebApproxInfo = struct('ninters',12,'nord',5,'domain',[-36,100]);
76[phi_wts,~,~]= prepWindowChebWts(phit1,phitt1,chebApproxInfo,tol);
78% npts = 7;
79% nqloc = npts;
80% [wts,LebesgueConstant] = get_lmmwts(xpts, npts, nqloc,rlen,phi_wts,h,chebApproxInfo);
81% bar(wts);
82% title(sprintf('Order %2d LMM weights for grid points\n with separation %1.2e',npts,h));
84p = 8;
85for i = 1:p
86 npts = 2*i;
87 nqloc = npts;
88 [wts,LebesgueConstant] = get_lmmwts(xpts, npts, nqloc,rlen,phi_wts,h,chebApproxInfo);
90 subplot((p/2),2,i)
91 bar(wts);
92 str = '\Lambda';
93 title(sprintf('order $%2d$, $%s_{%d} = %1.1f$',npts,str,npts,LebesgueConstant),'Interpreter','latex');
94end
96sgtitle(sprintf('LMM weights for grid points\n with separation %1.2e',h))
97saveas(gcf,'/Users/nalhassanieh/Desktop/stabilityResults/LMM_wts2','epsc');
98end