1% Copyright (c) 2015, The Chancellor, Masters and Scholars of the University 2% of Oxford, and the Chebfun Developers. 3% Copyright (c) 2016 The Gonum Authors 4% All rights reserved. 5% 6% Redistribution and use in source and binary forms, with or without 7% modification, are permitted provided that the following conditions are met: 8% * Redistributions of source code must retain the above copyright 9% notice, this list of conditions and the following disclaimer. 10% * Redistributions in binary form must reproduce the above copyright 11% notice, this list of conditions and the following disclaimer in the 12% documentation and/or other materials provided with the distribution. 13% * Neither the name of the University of Oxford nor the names of its 14% contributors may be used to endorse or promote products derived from 15% this software without specific prior written permission. 16% 17% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 18% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 19% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 20% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR 21% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 22% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 23% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 24% ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 25% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 26% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 27 28 29function [x, w, v] = hermpts(n, varargin) 30%HERMPTS Hermite points and Gauss-Hermite Quadrature Weights. 31% HERMPTS(N) returns N Hermite points X in (-inf, inf). By default these are 32% roots of the 'physicist'-type Hermite polynomials, which are orthogonal with 33% respect to the weight exp(-x.^2). 34% 35% HERMPTS(N, 'PROB') normalises instead by the probablist's definition (with 36% weight exp(-x.^2/2)), which gives rise to monomials. 37% 38% [X, W] = HERMPTS(N) returns also a row vector W of weights for Gauss-Hermite 39% quadrature. [X,W,V] = HERMPTS(N) returns in addition a column vector V of 40% the barycentric weights corresponding to X. 41% 42% [X, W] = HERMPTS(N, METHOD) where METHOD is one of 'GW', 'REC', 'GLR', or 43% 'ASY' allows the user to select which method is used. 'GW' will use the 44% traditional Golub-Welsch eigenvalue method [1], best when n<=20. 'REC' 45% uses Newton's method with polynomial evaluation via the 3-term 46% recurrence for Hermite polynomials. 'GLR' uses Glaser-Liu-Rokhlin 47% fast algorithm which is much faster for large N [2]. 'ASY' uses Newton's 48% method with polynomial evaluation via asymptotic formula. 'ASY' is the 49% fastest for N>=200, 'GLR' is the most accurate for nodes close to 0. 50% By default HERMPTS uses 'GW' when N <= 20, 'REC' for 21<=N<200, and 51% 'ASY' when N>=200. 52% 53% References: 54% [1] G. H. Golub and J. A. Welsch, "Calculation of Gauss quadrature 55% rules", Math. Comp. 23:221-230, 1969, 56% [2] A. Glaser, X. Liu and V. Rokhlin, "A fast algorithm for the 57% calculation of the roots of special functions", SIAM Journal 58% on Scientific Computing", 29(4):1420-1438:, 2007. 59% [3] A. Townsend, T. Trogdon and S. Olver, Fast computation of Gauss 60% nodes and weights on the whole real line, submitted, 2014. 61% 62% See also CHEBPTS, LEGPTS, LAGPTS, and JACPTS. 63 64% Copyright 2015 by The University of Oxford and The Chebfun Developers. 65% See http://www.chebfun.org/ for Chebfun information. 66% 67% 'GW' by Nick Trefethen, March 2009 - algorithm adapted from [1]. 68% 'GLR' by Nick Hale, March 2010 - algorithm adapted from [2]. 69 70% Defaults: 71method = 'default'; 72type = 'phys'; 73 74if ( n < 0 ) 75 error('CHEBFUN:hermpts:n', 'First input should be a positive integer.'); 76end 77 78% Return empty vector if n = 0. 79if ( n == 0 ) 80 [x, w, v] = deal([]); 81 return 82end 83 84% Check the inputs 85while ( ~isempty(varargin) ) 86 s = varargin{1}; 87 varargin(1) = []; 88 if ( strcmpi(s, 'GW') ) 89 method = 'GW'; 90 elseif ( strcmpi(s,'glr') ) 91 method = 'GLR'; 92 elseif ( strcmpi(s,'rec') ) 93 method = 'REC'; 94 elseif ( strcmpi(s,'asy') ) 95 method = 'ASY'; 96 elseif ( strncmpi(s, 'phys', 3) ) 97 type = 'phys'; 98 elseif ( strncmpi(s, 'prob', 3) ) 99 type = 'prob'; 100 else 101 error('CHEBFUN:hermpts:input', 'Unrecognised input string; %s.', s); 102 end 103end 104 105% Three cases: 106% 107% N <= 20: Use GW 108% 21<=N<200: Use REC 109% N>=200: Use ASY 110if ( n == 1 ) 111 % n = 1 case is trivial 112 x = 0; 113 w = sqrt(pi); 114 v = 1; 115 116elseif ( (n < 21 && strcmpi(method,'default')) || strcmpi(method,'GW') ) 117 % GW, see [1] 118 119 beta = sqrt(.5*(1:n-1)); % 3-term recurrence coeffs 120 T = diag(beta, 1) + diag(beta, -1); % Jacobi matrix 121 [V, D] = eig(T); % Eigenvalue decomposition 122 [x, indx] = sort(diag(D)); % Hermite points 123 w = sqrt(pi)*V(1, indx).^2; % weights 124 v = abs(V(1, indx)).'; % Barycentric weights 125 v = v./max(v); % Normalize 126 v(2:2:n) = -v(2:2:n); 127 128 % Enforce symmetry: 129 ii = 1:floor(n/2); 130 x = x(ii); 131 w = w(ii); 132 vmid = v(floor(n/2)+1); 133 v = v(ii); 134 if ( mod(n, 2) ) 135 x = [x ; 0 ; -x(end:-1:1)]; 136 w = [w, sqrt(pi) - sum(2*w), w(end:-1:1)]; 137 v = [v ; vmid ; v(end:-1:1)]; 138 else 139 x = [x ; -x(end:-1:1)]; 140 w = [w, w(end:-1:1)]; 141 v = [v ; -v(end:-1:1)]; 142 end 143 144elseif ( strcmpi(method,'GLR') ) 145 % Fast, see [2] 146 147 [x, ders] = alg0_Herm(n); % Nodes and H_n'(x) 148 w = (2*exp(-x.^2)./ders.^2)'; % Quadrature weights 149 v = exp(-x.^2/2)./ders; % Barycentric weights 150 v = v./max(abs(v)); % Normalize 151 if ( ~mod(n, 2) ) 152 ii = (n/2+1):n; 153 v(ii) = -v(ii); 154 end 155 156elseif ( (n < 200 && strcmpi(method,'default')) || strcmpi(method,'REC') ) 157 158 [x, w, v] = hermpts_rec( n ); 159 160else 161 162 [x, w, v] = hermpts_asy( n ); 163 164end 165 166% Normalise so that sum(w) = sqrt(pi) 167w = (sqrt(pi)/sum(w))*w; 168 169if ( strcmpi(type, 'prob') ) 170 x = x*sqrt(2); 171 w = w*sqrt(2); 172end 173 174end 175 176%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 177%% %%%%%%%%%%%%%%%%%%%%%%% Routines for GLR algorithm %%%%%%%%%%%%%%%%%%%%%%%%% 178%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 179 180% Driver for 'GLR'. 181function [roots, ders] = alg0_Herm(n) 182% Compute coefficients of H_m(0), H_m'(0), m = 0,..,N. 183 184Hm2 = 0; 185Hm1 = pi^(-1/4); 186Hpm2 = 0; 187Hpm1 = 0; 188for k = 0:n-1 189 H = -sqrt(k/(k+1))*Hm2; 190 Hp = sqrt(2/(k+1))*Hm1-sqrt(k/(k+1))*Hpm2; 191 Hm2 = Hm1; 192 Hm1 = H; 193 Hpm2 = Hpm1; 194 Hpm1 = Hp; 195end 196 197% allocate storage 198roots = zeros(n, 1); 199ders = zeros(n, 1); 200if ( mod(n,2) ) 201 % zero is a root: 202 roots((n-1)/2) = 0; 203 ders((n+1)/2) = Hp; 204else 205 % find first root: 206 [roots(n/2+1), ders(n/2+1)] = alg2_Herm(H,n); 207end 208 209% compute roots and derivatives: 210[roots, ders] = alg1_Herm(roots, ders); 211 212end 213 214%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 215 216% Main algorithm for 'GLR' 217function [roots, ders] = alg1_Herm(roots, ders) 218 219n = length(roots); 220s = mod(n, 2); 221N = (n - s) / 2; 222 223% number of terms in Taylor expansion 224m = 30; 225 226% initialise 227hh1 = ones(m + 1, 1); 228u = zeros(1, m + 1); 229up = zeros(1, m + 1); 230 231for j = (N + 1):(n - 1) 232 233 % previous root 234 x = roots(j); 235 236 % initial approx 237 h = rk2_Herm(pi/2,-pi/2,x,n) - x; 238 239 % scaling 240 M = 1/h; 241 242 % recurrence relation for Hermite polynomials 243 c1 = -(2*n+1-x^2)/M^2; 244 c2 = 2*x./M^3; 245 c3 = 1./M^4; 246 u(1) = 0; 247 u(2) = ders(j)/M; 248 u(3) = .5*c1*u(1); 249 u(4) = (c1*u(2) + c2*u(1))/6; 250 up(1) = u(2); 251 up(2) = 2*u(3)*M; 252 up(3) = 3*u(4)*M; 253 up(m+1) = 0; 254 255 for k = 2:m-2 256 u(k+3) = (c1*u(k+1) + c2*u(k) + c3*u(k-1))/((k+1)*(k+2)); 257 up(k+2) = (k+2)*u(k+3)*M; 258 end 259 260 % flip for more accuracy in inner product calculation 261 u = u(m+1:-1:1); 262 up = up(m+1:-1:1); 263 264 % Newton iteration 265 hh = hh1; 266 hh(end) = M; 267 step = inf; 268 l = 0; 269 z = zeros(m, 1); 270 while ( (abs(step) > eps) && (l < 10) ) 271 l = l + 1; 272 step = (u*hh)/(up*hh); 273 h = h - step; 274 % powers of h (This is the fastest way!) 275 hh = [M ; cumprod(M*h + z)]; 276 % flip for more accuracy in inner product calculation 277 hh = hh(end:-1:1); 278 end 279 280 % update 281 roots(j+1) = x + h; 282 ders(j+1) = up*hh; 283end 284 285% nodes are symmetric 286roots(1:N+s) = -roots(n:-1:N+1); 287ders(1:N+s) = ders(n:-1:N+1); 288 289end 290 291%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 292 293% find the first root (note H_n'(0) = 0) 294function [x1, d1] = alg2_Herm(Hn0, n) 295 296% advance ODE via Runge-Kutta for initial approx 297x1 = rk2_Herm(0, -pi/2, 0, n); 298 299% number of terms in Taylor expansion 300m = 30; 301 302% scaling 303M = 1/x1; 304% c = log10(n); 305% M = 1./x1.^(1-1.25/(c)); 306 307% initialise 308u = zeros(1,m+1); 309up = zeros(1,m+1); 310 311% recurrence relation for Legendre polynomials 312u(1) = Hn0; 313u(3) = -.5*(2*n+1)*u(1)/M^2; 314up(1) = 0; 315up(2) = 2*u(3)*M; 316for k = 2:2:m-2 317 u(k+3) = (-(2*n+1)*u(k+1)/M^2 + u(k-1)/M^4)/((k+1)*(k+2)); 318 up(k+2) = (k+2)*u(k+3)*M; 319end 320 321% flip for more accuracy in inner product calculation 322u = u(m+1:-1:1); 323up = up(m+1:-1:1); 324 325z = zeros(m, 1); 326x1k = [M ; cumprod(M*x1 + z)]; 327step = inf; 328l = 0; 329% Newton iteration 330while ( (abs(step) > eps) && (l < 10) ) 331 l = l + 1; 332 step = (u*x1k)/(up*x1k); 333 x1 = x1 - step; 334 % powers of h (This is the fastest way!) 335 x1k = [1 ; cumprod(M*x1 + z)]; 336 x1k = x1k(end:-1:1); 337end 338 339% Update derivative 340d1 = up*x1k; 341 342end 343 344%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 345 346% Runge-Kutta for Hermite Equation 347function x = rk2_Herm(t, tn, x, n) 348m = 10; 349h = (tn-t)/m; 350for j = 1:m 351 k1 = -h/(sqrt(2*n+1-x^2) - .5*x*sin(2*t)/(2*n+1-x^2)); 352 t = t + h; 353 k2 = -h/(sqrt(2*n+1-(x+k1)^2) - .5*x*sin(2*t)/(2*n+1-(x+k1)^2)); 354 x = x + .5*(k1 + k2); 355end 356end 357 358%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 359%% %%%%%%%%%%%%%%%%%%%%%%% Routines for ASY algorithm %%%%%%%%%%%%%%%%%%%%%%%%% 360%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 361 362function [x, w, v] = hermpts_asy(n) 363% HERMPTS_ASY, fast algorithm for computing Gauss-Hermite nodes and weights 364% using Newton's method with polynomial evaluation via asymptotic expansions. 365% 366% x = Gauss-Hermite nodes, w = quad weights, v = bary weights. 367% 368% See [3]. 369 370[x, w, v] = hermpts_asy0( n ); 371 372if mod(n,2) == 1 % fold out 373 x = [-x(end:-1:1);x(2:end)]; 374 w = [w(end:-1:1) w(2:end)]; w = (sqrt(pi)/sum(w))*w; 375 v = [v(end:-1:1);v(2:end)]; v = v./max(abs(v)); 376else 377 x = [-x(end:-1:1);x]; 378 w = [w(end:-1:1) w]; w = (sqrt(pi)/sum(w))*w; 379 v = [v(end:-1:1);-v]; v = v./max(abs(v)); 380end 381 382% debug 383%tic, exact = hermpts(n); toc 384%semilogy(abs(exact-x)) 385end 386 387function [x, w, v] = hermpts_asy0(n) 388% Compute Hermite nodes and weights using asymptotic formula 389 390x0 = HermiteInitialGuesses(n); % get initial guesses 391t0 = x0./sqrt(2*n+1); 392theta0 = acos(t0); % convert to theta-variable 393for k = 1:20 394 [val, dval] = hermpoly_asy_airy(n, theta0); 395 dt = -val./(sqrt(2)*sqrt(2*n+1)*dval.*sin(theta0)); 396 theta0 = theta0 - dt; % Newton update 397 if norm(dt,inf) < sqrt(eps)/10, break; end 398end 399t0 = cos(theta0); 400x = sqrt(2*n+1)*t0; % back to x-variable 401 402ders = x.*val + sqrt(2)*dval; 403%ders = dval; 404w = (exp(-x.^2)./ders.^2)'; % quadrature weights 405 406v = exp(-x.^2/2)./ders; % Barycentric weights 407end 408 409function [val, dval] = hermpoly_asy_airy(n, theta) 410% HERMPOLY_ASY evaluation hermite poly using Airy asymptotic formula in 411% theta-space. 412 413musq = 2*n+1; 414cosT = cos(theta); sinT = sin(theta); 415sin2T = 2*cosT.*sinT; 416eta = .5*theta - .25*sin2T; 417chi = -(3*eta/2).^(2/3); 418phi = (-chi./sinT.^2).^(1/4); 419const = 2*sqrt(pi)*musq^(1/6)*phi; 420Airy0 = real(airy(musq.^(2/3)*chi)); 421Airy1 = real(airy(1,musq.^(2/3)*chi)); 422 423% Terms in (12.10.43): 424a0 = 1; b0 = 1; 425a1 = 15/144; b1 = -7/5*a1; 426a2 = 5*7*9*11/2/144^2; b2 = -13/11*a2; 427a3 = 7*9*11*13*15*17/6/144^3; 428b3 = -19/17*a3; 429 430% u polynomials in (12.10.9) 431u0 = 1; u1 = (cosT.^3-6*cosT)/24; 432u2 = (-9*cosT.^4 + 249*cosT.^2 + 145)/1152; 433u3 = (-4042*cosT.^9+18189*cosT.^7-28287*cosT.^5-151995*cosT.^3-259290*cosT)/414720; 434 435%first term 436A0 = 1; 437val = A0*Airy0; 438 439%second term 440B0 = -(a0*phi.^6.*u1+a1*u0)./chi.^2; 441val = val + B0.*Airy1./musq.^(4/3); 442 443%third term 444A1 = (b0*phi.^12.*u2 + b1*phi.^6.*u1 + b2*u0)./chi.^3; 445val = val + A1.*Airy0/musq.^2; 446 447%fourth term 448B1 = -(phi.^18.*u3 + a1*phi.^12.*u2 + a2*phi.^6.*u1 + a3*u0)./chi.^5; 449val = val + B1.*Airy1./musq.^(4/3+2); 450 451val = const.*val; 452 453%% Derivative 454 455eta = .5*theta - .25*sin2T; 456chi = -(3*eta/2).^(2/3); 457phi = (-chi./sinT.^2).^(1/4); 458const = sqrt(2*pi)*musq^(1/3)./phi; 459 460% v polynomials in (12.10.10) 461v0 = 1; v1 = (cosT.^3+6*cosT)/24; 462v2 = (15*cosT.^4-327*cosT.^2-143)/1152; 463v3 = (259290*cosT + 238425*cosT.^3 - 36387*cosT.^5 + 18189*cosT.^7 -... 464 4042*cosT.^9)/414720; 465 466%first term 467C0 = -(b0*phi.^6.*v1 + b1.*v0)./chi; 468dval = C0.*Airy0/musq.^(2/3); 469 470% %second term 471D0 = a0*v0; 472dval = dval + D0*Airy1; 473 474% %third term 475C1 = -(phi.^18.*v3 + b1*phi.^12.*v2 + b2*phi.^6.*v1 + b3*v0)./chi.^4; 476dval = dval + C1.*Airy0/musq.^(2/3+2); 477 478%fourth term 479D1 = (a0*phi.^12.*v2 + a1*phi.^6.*v1 + a2*v0)./chi.^3; 480dval = dval + D1.*Airy1/musq.^2; 481 482dval = const.*dval; 483 484end 485 486function x_init = HermiteInitialGuesses(n) 487%HERMITEINTITIALGUESSES(N), Initial guesses for Hermite zeros. 488% 489% [1] L. Gatteschi, Asymptotics and bounds for the zeros of Laguerre 490% polynomials: a survey, J. Comput. Appl. Math., 144 (2002), pp. 7-27. 491% 492% [2] F. G. Tricomi, Sugli zeri delle funzioni di cui si conosce una 493% rappresentazione asintotica, Ann. Mat. Pura Appl. 26 (1947), pp. 283-300. 494 495% Gatteschi formula involving airy roots [1]. 496% These initial guess are good near x = sqrt(n+1/2); 497if mod(n,2) == 1 498 m = (n-1)/2; bess = (1:m)'*pi; a = .5; 499else 500 m = n/2; bess = ((0:m-1)'+.5)*pi; a = -.5; 501end 502nu = 4*m + 2*a + 2; 503T = @(t) t.^(2/3).*(1+5/48*t.^(-2)-5/36*t.^(-4)+(77125/82944)*t.^(-6) -... 504 108056875/6967296*t.^(-8)+162375596875/334430208*t.^(-10)); 505airyrts = -T(3/8*pi*(4*(1:m)'-1)); 506 507airyrts_exact = [ -2.338107410459762 % Exact Airy roots. 508 -4.087949444130970 509 -5.520559828095555 510 -6.786708090071765 511 -7.944133587120863 512 -9.022650853340979 513 -10.040174341558084 514 -11.008524303733260 515 -11.936015563236262 516 -12.828776752865757]; 517airyrts(1:10) = airyrts_exact; % correct first 10. 518 519x_init = sqrt(nu + 2^(2/3)*airyrts*nu^(1/3) +... 520 1/5*2^(4/3)*airyrts.^2*nu^(-1/3) +... 521 (11/35-a^2-12/175*airyrts.^3)/nu +... 522 (16/1575*airyrts+92/7875*airyrts.^4)*2^(2/3)*nu^(-5/3) -... 523 (15152/3031875*airyrts.^5+1088/121275*airyrts.^2)*2^(1/3)*nu^(-7/3)); 524x_init_airy = real(x_init(end:-1:1)); 525 526% Tricomi initial guesses. Equation (2.1) in [1]. Originally in [2]. 527% These initial guesses are good near x = 0 . Note: zeros of besselj(+/-.5,x) 528% are integer and half-integer multiples of pi. 529% x_init_bess = bess/sqrt(nu).*sqrt((1+ (bess.^2+2*(a^2-1))/3/nu^2) ); 530Tnk0 = pi/2*ones(m,1); 531nu = (4*m+2*a+2); 532rhs = (4*m-4*(1:m)'+3)./nu*pi; 533 534for k = 1:7 535 val = Tnk0 - sin(Tnk0) - rhs; 536 dval = 1 - cos(Tnk0); 537 dTnk0 = val./dval; 538 Tnk0 = Tnk0 - dTnk0; 539end 540 541tnk = cos(Tnk0/2).^2; 542x_init_sin = sqrt(nu*tnk - (5./(4*(1-tnk).^2) - 1./(1-tnk)-1+3*a^2)/3/nu); 543 544% Patch together 545p = 0.4985+eps; 546x_init = [x_init_sin(1:floor(p*n));x_init_airy(ceil(p*n):end)]; 547 548 549if mod(n,2) == 1 550 x_init = [0;x_init]; 551 x_init = x_init(1:m+1); 552else 553 x_init = x_init(1:m); 554end 555 556% debug: 557%y = hermpts(n); 558%semilogy(abs(y - x_init)); 559%yhalf = -y(m:-1:1); 560%semilogy(abs(yhalf - x_init)); 561end 562 563%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 564%% %%%%%%%%%%%%%%%%%%%%%%% Routines for REC algorithm %%%%%%%%%%%%%%%%%%%%%%%%% 565%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 566 567function [x, w, v] = hermpts_rec(n) 568% Compute Hermite nodes and weights using recurrence relation. 569 570x0 = HermiteInitialGuesses(n); 571x0 = x0.*sqrt(2); 572 573for kk = 1:10 574 [val, dval] = hermpoly_rec(n, x0); 575 dx = val./dval; 576 dx(isnan(dx)) = 0; 577 x0 = x0 - dx; 578 if norm(dx, inf)<sqrt(eps), break; end 579end 580x = x0/sqrt(2); 581w = (exp(-x.^2)./dval.^2)'; % quadrature weights 582v = exp(-x.^2/2)./dval; % Barycentric weights 583 584if mod(n,2) == 1 % fold out 585 x = [-x(end:-1:1);x(2:end)]; 586 w = [w(end:-1:1) w(2:end)]; w = (sqrt(pi)/sum(w))*w; 587 v = [v(end:-1:1);v(2:end)]; v = v./max(abs(v)); 588else 589 x = [-x(end:-1:1);x]; 590 w = [w(end:-1:1) w]; w = (sqrt(pi)/sum(w))*w; 591 v = [v(end:-1:1);-v]; v = v./max(abs(v)); 592end 593 594 595end 596 597function [val, dval] = hermpoly_rec(n, x0) 598% HERMPOLY_rec evaluation of scaled Hermite poly using recurrence 599 600% evaluate: 601Hold = exp(-x0.^2/4); H = x0.*exp(-x0.^2/4); 602for k = 1:n-1 603 Hnew = (x0.*H./sqrt(k+1) - Hold./sqrt(1+1/k)); 604 Hold = H; H = Hnew; 605end 606% evaluate derivative: 607val = Hnew; 608dval = (-x0.*Hnew + n^(1/2)*Hold); 609end