# MATLAB Codes

There are numerous MATLAB resources for spectral and pseudospectral methods.
(1) Gautschi, W. Algorithm 726: ORTHPOL—A Package of Routines for Generating Orthogonal Polynomials and Gauss-Type Quadrature Rules, ACM Trans. Math. Software 20, 21-62 (1994). The source codes are here.

(2) MATLAB Differention Suite by Weideman and Reddy ACM Trans. Math. Sci. 26, 465-519 (2000). The source codes are here.

(3) Trefethen, L. N. Spectral Methods in Matlab, SIAM (2000).

(4) Gautschi, W. Orthogonal Polynomials (in MATLAB), J. Comput. Appl. Math 178, 215-234 (2005).

(5) CHEBFUN suite of Matlab codes from Oxford University at chebfun.org

MATLAB codes that accompany Spectral Methods in Chemistry and Physics.

Chapter 3   Chapter 4   Chapter 5   Chapter 6    (Click on a chapter to be redirected.)

Chapter 2 – Click on a file and save, changing the file extension from *.txt to *.m.

 ab_bimode Calculation of recurrence coefficients, βn, for full range bimodal polynomials written to bimodal.dat. Coefficients for ε = 0.02 and 0.005 are in bimodal_eps02 and bimodal_eps005 and shown graphically in beta_n_eps02 and beta_n_eps005. bimode_polyn_with_pts_eps02 Calculation of full range bimodal polynomial of order n for ε = 0.02 with recurrence coefficients in bimodal_eps02.dat. sample output bimodal_polyn_with_pts_eps02. Quadrature points are shown by the solid circles. Maple code bimode-beta_n-eps02-Maple bimode_polyn_with_pts_eps005 Calculation of full range bimodal polynomial of order n for ε = 0.005 with recurrence coefficients in bimodal_eps005.dat. sample output bimodal_polyn_with_pts_eps005. Quadrature points are shown by the solid circles. LnSqwtspts MultiExp quadrature by Gill and Chien (2003). Sample output for N = 10, LnSqwtspts10 – see also Gill web site ab_maxwell_p0 Recurrence coefficients, αn and βn, for Maxwell polynomials, p=0; output data abmaxp0 – more precise data from Gautschi abhrhermite ab_maxwell_p2 Recurrence coefficients, αn and βn, for Maxwell polynomials, p=2; output data abmaxp2 maxwell_p0_with_pts Maxwell polynomial for p = 0 and  N = 10 in the graph Maxwell_p0_with_pts compared with Hermite  (N = 20) in Fig. 2.4 for x positive. Rys_pw Recurrence coefficients, αn and βn, for Rys polynomials, xmin = -1 for full range, J_n(x), and xmin = 0 for half range, R_n(x). For N = 13, xmin = -1, recurrence coefficients Rys_fr_ab1 and quadrature points Rys_fr_pw1  are to be compared with Sagar and Smith. For N = 40, xmin = -1, recurrence coefficients Rys_fr_ab2 and quadrature points Rys_fr_pw2  are to be compared with Schneider and Nygaard. Rys_Jn_with_pts Full range, J_n(x), n = 20, c=2, polynomials as shown in Fig. 2.9  Jn_with_pts. Polynomials are not normalized to unity. Rys_Rn_with_pts Half range, R_n(x), n = 10, c = 2, polynomials as shown in Fig. 2.9 Rn_with_pts. Polynomials are not normalized to unity. ab_GK_rad_transfer Recurrence coefficients, αn and βn, for w(x)=exp(-3/2x) introduced by Gander and Karp for radiative transfer problems. Output in ab_GK is to be compared with results with Maple in Table 5 in Gander and Karp Legptswts Calculation of Legendre quadrature weights and points. Sample output for N = 10, Legptswts10 – see also the Keisan online calculator Hermptswts Calculation of Hermite quadrature weights and points. Sample output for  N = 10, Hermptswts10 – see also the Keisan online calculator Lagptswts Calculation of Laguerre quadrature weights and points. Sample output for  N = 10, Lagptswts10 – see also the Keisan online calculator lagpoly_with_pts Graph for the Associated Laguerre polynomial with quadrature points – click to download sample output hermite_with_pts Graph for the Hermite polynomials with quadrature points – click to download sample output Legendre_with_pts Graph for the Legendre polynomials with quadrature points – click to download sample output

Chapter 3: – Click on a file and save, changing the file extension from *.txt to *.m.

 legendre_norm_linear Quadrature evaluation of the normalization of two Maxwell polynomials, n=4 and n=6. Convergence with the number of Legendre quadrature points, N, and a linear map to [-1,1] as shown in Fig. 3.3. legendre_norm_expmap Quadrature evaluation of the normalization of two Maxwell polynomials, n=4 and n=6. Convergence with the number of Legendre quadrature points, N, and an exponential map to [-1,1] as shown in Fig. 3.3. legendre_norm_Becke Quadrature evaluation of the normalization of two Maxwell polynomials, n=4 and n=6. Convergence with the number of Legendre quadrature points, N, and the Boyd-Becke map to [-1,1] as shown in Fig. 3.3. legendre_norm_MK Quadrature evaluation of the normalization of two Maxwell polynomials, n=4 and n=6. Convergence with the number of Legendre quadrature points, N, and the Mura-Knowles map to [-1,1] as shown in Fig. 3.3. radint_maxp2_vsN_a Convergence of the radial integral, Eq. (3.10), for 2 Gaussians, Eq. (3.14), versus N for several scale factors s as shown in Fig. 3.6(A). The linear fit parameters A(s) and B(s)  are printed to the screen and are shown in Table 3.2. radint_maxp2_vsN_b Convergence of the radial integral, Eq. (3.10), for 3 Gaussians, Eq. (3.15), versus N for several scale factors s as shown in Fig. 3.6(B). The linear fit parameters A(s) and B(s)  are printed to the screen and are shown in Table 3.2. lindh_integ_max2 Integration of a Gaussian,  Eq. (3.11), with Maxwell  (p = 2) quadratures; with results in Fig.3.4. Data file with recurrence coefficients for Maxwell quadratures is in abmaxp2 lindh_integ_max2B Integration of a Gaussian, Eq. (3.11),  with Maxwell (p = 2) quadratures versus the number of quadrature points, N for different values. sigma_O_O2 Energy dependence of the reactive cross section, Eq. (3.25), for collisional dissociation of O-O2 reactive collisions in comparison with the line-of-centers cross section, Eq. (3.23), as shown in Fig. 3.7. integrand_O_O2 Variation of the integrand for the calculation of the equilibrium O-O2 dissociative rate coefficient as shown in Fig. 3.7 for three temperatures. sr_rate_O_O2 Convergence of the Simpson’s rule integration of the O-O2 dissociative rate coefficient as shown in Table 3.4. lag_rate_O_O2 Convergence of the equilibrium O-O2 dissociative rate coefficient with Laguerre quadratures as shown in Table 3.4. lag_scaled_rate_O_O2 Convergence of the equilibrium O-O2 dissociative rate coefficient with Laguerre quadratures as shown in Table 3.4 versus the scaling parameter, s. nurx_tdep_speed This code evaluates the thermal average, Eq. (3.22), of the nuclear cross section, Eq. (3.33). The convergence of the nuclear fusion reaction rate evaluated with a Maxwell quadrature based on the weight function w(x) = xexp(-x*x). Inp ut parameters  are;  tk = temperature in 10**9 K, sc = quadrature scale factor and irx identifies a reaction from a choice of four. The results are in Table 3.7

Chapter 4: – Click on a file and save, changing the file extension from *.txt to *.m.

 sin_expand_hermite Expansion of sin(x) in Hermite polynomials; Fig. 4.4. The MATLAB code is run three times to produce this figure. Some changes by hand are required between each run to format the graph as noted. plot_kappa Plot_kappa plots the Kappa distribution for 5 kappa values in comparison with the Maxwellian expand_kappa_lag Expansion of the Kappa distribution in Laguerre polynomials demonstrating the nonconvergence; Fig. 4.8 for kappa = 30 with s = 1 and nmax = 25. For kappa = 50, the plotting format has to be changed. wave_packet_cosines Plot the graph for the sum of cosines,  Eq. (4. 105) as shown in Fig. 4.12 wave_packet A model for a quantum mechanical wavepacket as shown in Fig. 4.13 fid A model free induction decay curve and the Fourier transform are shown in Fig.4.14. The frequency spectrum of the signal is recovered with the FFT. fs_step Fourier sine series of a step function shown in Fig. 4.15 as an introduction to Gibbs phenomenon runge_uniform Demonstration of the Runge phenomenon with a uniform interpolation grid for the function, Eq. (4.151)  as shown in Fig. 4.20 (upper graphs). runge_cheb Interpolation of the function, Eq. (4.151) with Chebyshev quadrature points as shown in Fig. 4.20 (lower graph)

Chapter 5: – Click on a file and save, changing the file extension from *.txt to *.m.

 noneq_rx Nonequilibrium effects for an elementary chemical reaction as shown in Fig. 5.4(A) noneq_spectral Spectral convergence of the solution of the Boltzmann equation for an elementary chemical reaction as summarized in Fig. 5.4(C) linbe_gr_eig Eigenvalues of the linearized Boltzmann collision operator, Eq. (5.65). Spectral convergence of the “spectral gap”, λ2. Sonine-Laguerre basis functions are used. linbe_eig Eigenvalues of the linearized Boltzmann collision operator with the pseudospectral method based on the Maxwell polynomial. Maxwell polynomial quadratures are used, Eqs. (5.99) and (5.100). This requires the recurrence coefficients up to N = 200 in ab_maxp2_200 as a *.dat file.(This code could be further vectorized.) linbe_eig_gamma1 Eigenvalues of the linear Boltzmann collision operator, Eqs. (5.108) and (5.110), for mass ratio =1. This requires the recurrence coefficients up to n = 200 in ab_maxp2_200 as a *.dat file. (This code could be further vectorized.) rt_chandra Eigenvalues of the radiative transfer equation, Eq. (5.19), and the extrapolation length. Eq. (5.25)

Chapter 6: – Click on a file and save, changing the file extension from *.txt to *.m.

 Harm_Osc_Fourier Convergence of the eigenvalues for the quantum harmonic oscillator with Fourier basis functions as shown in Fig. 6.14(A). Harm_Osc_Herm1 Eigenvalues of the quantum harmonic oscillator with Hermite basis functions with the representation as in Eq. (6.175) as in Baye and Heenen (1986). The results obtained with this Matlab code shows “ghost” levels as in Table 6.6. Harm_Osc_Herm2 Exact calculation of the eigenvalues of the quantum harmonic oscillator with Hermite basis functions with the representation as in Eq. (6.178). The derivative operator is constructed as in Eq. (3.138) without explicit reference to the harmonic potential. The eigenvalues are exact as are the eigenfunctions  shown in Fig.  (6.15) calculated with this code. SE_electron_potl The convergence of the eigenvalues for the Schroedinger equation with the potential, Eq. (6.180),  based on the  nonclassical polynomials orthogonal with respect to the weight function, Eq. (6.181). The recurrence coefficients for this weight function are in the file ab_se_electron. This Schroedinger equation is isospectral with the Fokker-Planck equation for electron relaxation for a hard sphere cross section, namely Eq. (6.63), with sigma(x) = 1