Path integral Monte Carlo method for the quantum anharmonic oscillator

The Markov chain Monte Carlo (MCMC) method is used to evaluate the imaginary-time path integral of a quantum oscillator with a potential that includes both a quadratic term and a quartic term whose coupling is varied by several orders of magnitude. This path integral is discretized on a time lattice and calculations for the energy and probability density of the ground state and energies of the first few excited states are carried out on lattices with decreasing spacing to estimate these quantities in the continuum limit. The variation of the quartic coupling constant produces corresponding variations in the optimum simulation parameters for the MCMC method and in the statistical uncertainty for a fixed number of paths used for measurement. The energies and probability densities are in excellent agreement with those obtained from numerical solutions of Schr\"odinger's equation.


I. INTRODUCTION
The harmonic oscillator is one of the few exactly solvable quantum mechanical systems. Solutions for the energy eigenvalues and wave functions are obtained either by the analytic solution of the time-independent Schrödinger equation in terms of Hermite polynomials, or through a more abstract algebraic procedure based on raising and lowering operators 1 . The operator method makes a direct connection to the statistics of bosons, a conceptual building block of modern physics.
The bosonic character of the quantum harmonic oscillator is used in the quantization of smallamplitude vibrations in molecules and solids, the latter leading to phonons 2 , the quantum theory of radiation [3][4][5] , quantum field theory 6 , and as an illustration of the correspondence principle 7 . Such diverse applications attest to the pivotal role of the harmonic oscillator in the conceptual and computational development of quantum physics.
Sufficiently small fluctuations in any system around a stable equilibrium point may be described in terms of decoupled harmonic oscillators (normal modes), regardless of the shape of the confining potential. However, there are regimes where the harmonic oscillator paradigm breaks down. For example, during the thermal expansion of solids, the transformations between solid phases, and chemical reactions, the displacements of atoms from their equilibrium positions cannot be regarded as small. As harmonic interactions are derived by truncating the Taylor series expansion of the interatomic potential at second order 2 , we look to higher-order terms to augment the interaction potential. A basic one-dimensional model that incorporates the leading (quartic) correction to the harmonic potential has the Schrödinger equation, in which m is the mass of the particle, ω = (k/m) 1/2 is the natural frequency of the harmonic oscillator, where k is the stiffness of the potential, and λ is the coupling constant for the quartic term of the potential.
In the absence of an elementary analytic method for obtaining solutions of (1), attention turned to various approximate calculations. In Rayleigh-Schrödinger perturbation theory, the ground-state energy E 0 (λ) of (1) has the formal expansion where the expansion coefficients A n must be determined. Bender and Wu calculated 8 these coefficients to 75th order and found a rapid increase in |A n |, which suggests that (2) diverges. In fact, they showed that this series diverges for any λ > 0 because of the singularity at λ = 0 that separates the regions where (1) has an infinite sequence of bound states (λ > 0) from the region where there are only metastable states that can escape to infinity (λ < 0) [8][9][10] . This explanation is a particular case of an argument first advanced by Dyson 11 for perturbation expansions in quantum electrodynamics. The divergence of (2) has provided the impetus for alternative perturbation expansions 12 , resummation methods 13 , and other computational schemes 14,15 for quantum anharmonic oscillators.
In this paper, we adopt a somewhat different approach to solving (1) by evaluating the imaginarytime path integral for this system using the Markov chain Monte Carlo (MCMC) method 16 . Calculations are carried out on lattices with decreasing spacing to estimate the energy and probability density of the ground state and energies of low-lying excited states in the continuum limit. Comparisons with numerical integrations of Schrödinger's equation demonstrate the accuracy of our approach. These calculations are also pedagogical: the wide range of λ-values in (1) used in this study help to develop an intuition about how the potential affects parameters and convergence.
Although we are studying a specific case, the MCMC method can be applied to a quantum mechanical particle in any potential with minimal change to the basic procedure described here 16 .
Imaginary-time time path integrals can also be formulated for classical dynamics [17][18][19][20] , as well as providing a bridge to quantum field theory 21 .
The organization of this paper is as follows. We briefly outline the derivation of real-time and imaginary-time path integrals in Sec. B, referring to our earlier work 22,23 for a more comprehensive discussion. The MCMC method is summarized in Sec. C, and in Sec. D we discuss how varying the coupling constant of the quartic term in the potential energy affects the parameters in this method. Section E presents the correlation functions that can be calculated from the imaginary-time formalism to obtain the energy eigenvalues of (1). We have calculated the energy and probability density of the ground state and the energies of the first two excited states. We summarize our results and discuss other applications of imaginary-time path integrals in Sec. E.

A. Feynman path integral
The time-dependent Schrödinger equation, for a Hamiltonian operatorĤ =p The completeness relation, for the orthonormal eigenketsĤ|n = E n |n , when used in the trace in (B8), yields the partition function in the usual form

C. Correlation functions and propagators
The energies of the ground state and excited states of the quantum anharmonic oscillator are encoded in correlation functions that are expectation values of products of the position operator: The energy E 0 of the ground state can be obtained from the expectation of the Hamiltonian (1) which, in conjunction with the virial theorem 1 , is expressed in terms of correlation functions of x An alternative expression for the group-state energy, more closely related to those for excited states derived below, is The correlation functions x n for n = 2 and 4 are, from (13), By writingx which is the imaginary time counterpart of the Heisenberg picture for the time-dependence of operators, again using (11), and taking the limit β → ∞, we obtain which is the expectation value ofx n in the ground state.
The two-point correlation function for calculating the first excited state is where the subscript "c " denotes cumulant, which in diagrammatic analysis correspond to a connected diagram. Again using (17) and (11), yields Hence, the energy of the first excited state is obtained as The second excited state is obtained from the four-point connected correlation function By proceeding as above, we obtain The second excited state is, therefore, determined from The probability density of the ground-state wave function is obtained by following the development of Creutz and Freedman 16 . The probability P (x) of a particle being found between positions x − 1 2 ∆x and x + 1 2 ∆x at any time t in a real-time interval [0, t] is given by the time-average The numerator in the integrand counts the paths that begin at (x i , 0), end at (x f , t), and pass through (x , t ) for 0 ≤ t ≤ t. The propagator in the denominator counts all paths between (x i , 0) and (x f , t).
If ∆x is assumed to be small enough so that the integral over x can be evaluated by keeping terms only to first order in ∆x, we obtain The propagators are now written in terms of the eigenfunctions ofĤ by using the completeness relation (11) for the eigenkets ofĤ: By using this expression for each propagator in (26) and continuing to imaginary time τ = β, the long-imaginary-time/low-temperature limit yields the probability density of the ground state: Figure 1 provides a schematic illustration of how the wave function is calculated by assigning a path along a time lattice to spatial bins with width ∆x.

III. MARKOV CHAIN MONTE CARLO METHOD
Monte Carlo simulations are carried out on a time lattice with N τ time increments δτ with lattice points x n = nδτ for n = 0, 1, 2, . . . , N τ . The nth time increment is x n+1 − x n . Periodic boundary conditions are imposed on the lattice, whereby increment N τ + 1 equals increment 1. The imaginary-time path integral in (8), and (9) is the continuum limit of the discretized expression where Appendix A. Probability Density 40 The cross terms are O(1/ E 10 ), where E 10 = E 1 E 0 . By arguments similar to This is the Feynman-Kac formula 24,25,30 , which establishes a connection between the propagator (the Green's function for the Hamiltonian) and a path integral.
We work with a dimensionless action. With each variable expressed in terms of a suitable power of the lattice spacing δτ , we work in units where = 1 = c, and introduce the dimensionless By combining (1), (30), and (31), the action for the anharmonic oscillator becomes in whichλ is also dimensionless.
The MCMC method is based on determining the statistics of observables from paths (x 1 , . . . , x Nτ ) that are representative of the distribution in (29) and (30). This requires generating reliable sequences of (pseudo) random numbers. We have used the Mersenne twister 33 .
We begin with an initial path P (0) , which may be an array of random numbers ("hot" start) or zeros ("cold" start). This path is updated by applying the Metropolis-Hastings algorithm 31,32 to each element x i of the path in random order, called a "sweep." There are two steps in the updating process: 1. Generate a random number u from a uniform distribution in the interval [−h, h], where h, called the hit size, must be chosen judiciously. If h is too large, few changes will be accepted; too small and the exploration of phase space will be slow. For the sets of simulations here, hit sizes were chosen to obtain an acceptance rate of 50-60%. One sweep produces the next path, e.g. P (1) from P (0) . Each path is determined only by the immediately preceding path, so the complete sequence of paths forms a Markov chain, but the paths are correlated. The accuracy of the MCMC method relies on sampling from the distribution in (29) and (30) which, in turn, relies on the paths being stationary and independent. The initial path "thermalizes", that is, attained equilibrium after N therm sweeps. To counteract the inherent autocorrelation in a Markov chain, a number N sep of paths between successive paths used for

IV. PARAMETERS FOR MCMC SIMULATIONS
For all calculations reported here,m =ω = δτ and N τ δτ = 250 in (32). Calculations have been performed forλ = 0, 1, 50 and 10 3 , which range from the harmonic oscillator to the strong quartic limit (Fig. 2). Such a large variation ofλ affects not just the quantum mechanical behavior of the oscillator, but also several parameters used in the MCMC method: the hit size h needed to achieve an acceptance rate of 50-60%, and the number N sep of paths that must be discarded between successive paths used for calculations. The coupling constant also affects the convergence to the continuum limit of the probability density of the ground state and the energy levels. Figure 3 compares the acceptance rate versus hit size for the anharmonic oscillator with a weak (λ = 1) and a strong (λ = 10 3 ) quartic coupling constant for several discretizations. These data were obtained as follows: (i) begin with a cold start and h = 0.2, (ii) ignore several initial sweeps, (iii) calculate the acceptance rate, as given in the pseudocode in Ref. 22, (iv) calculate the arithmetic mean of the acceptance rate for every 100 sweeps. The hit size is then increased by 0.2 and the process (i)-(iv) is repeated.
Most apparent in Fig. 3 is that, for each discretization, the target hit size is much smaller for the larger quartic coupling constant, which leads to larger changes in the action with increasing  coupling and, thus, suppresses the acceptance rate. In other words, the exploration of phase space is slower, which is a natural consequence of the more localized potential associated with a stronger coupling constant (Fig. 2).
A key element of the MCMC method is the selection of paths used for calculations. These paths must be representative of the equilibrium distribution, as defined by the partition function, so the paths must first equilibrate from some initial configuration. The number N therm of sweeps required to attain equilibrium is determined when the measured quantity fluctuates about a steady state, which depends on the observable, but generally increases as either δτ orλ decreases. Typically, the equilibration of several observables is plotted and the maximum number of sweeps to equilibrium is used for all simulations. Figure 4 shows the equilibration for x 2 for the anharmonic oscillator with coupling constants λ = 1 (a,b) andλ = 10 3 (c,d) for lattice spacings δτ = 0.2 (a,c) and δτ = 1 (b,d). Note the differences in scales of x 2 in each panel, which shows that the equilibrium value is much smaller for the system with the larger coupling constant. Also immediately apparent is how many fewer sweeps from the initial path are required to attain equilibrium with increasing coupling constant.
The trends in Figs. 3 and 4 can be explained by the sharpening of the potential energy profile with increasingλ (Fig. 2). All of our simulations began with cold starts. Thus, with increasing λ, the accessible configuration space decreases, so the attainment of equilibrium is correspondingly quicker. Similarly, the acceptance rate for a fixed hit size h is lower for largerλ because of the increasing confinement near the origin, at the expense of the classically forbidden region. This also explains why, for an observable such as x 2 (Fig. 4), sweeps over Markov chains show smaller fluctuations around their mean for largerλ. Finally, the increasing equilibrium value of x 2 withλ is a direct result of the increasingly localized probability density of the ground-state wave function (Sec. V A) caused by the sharpening potential energy profile (Fig. 2).
The number N sep of sweeps discarded between successive measurements, and the hit size h for the discretizations δτ used for the calculations described in the following section are compiled in Table I  The properties of our Markov chains enable an indirect evaluation of the long-imaginary-time limit (32) of the expression (29) for the probability density of the ground state. In particular, the Markov chains in Sec. C, which are aperiodic, irreducible, and positively recurrent, guarantee not only that any initial chain approaches equilibrium, as the discussion accompanying Fig. 4 demonstrates, but that long-time limits may be replaced by ensemble averages (ergodicity) 35,36 . the enhanced localization of the wave function forλ = 10 3 compared with that for the harmonic oscillator (λ = 0). This is to be expected from the narrowing of the potential with increasingλ ( Fig. 1). Moreover, the rate of convergence to the continuum limit is considerably slower for the large quartic term. The errors for δτ = 1 are small for the harmonic oscillator, but are substantial for the large quartic potential, with small discrepancies remaining in the tail of the distribution even for δτ = 0.1.

B. Ground-state energy
The energies of anharmonic oscillators were expressed as correlation functions in Sec. II C. The MCMC method necessitates evaluating discrete approximations to these quantities, which are then used to estimate the energies at decreasing lattice spacing δτ . The results with the simulation parameters in Table I are shown in Fig. 6 forλ = 0, 1, 50, and 10 3 . Also shown are the exact Next, we will explore the plots of the ground state probability density (see Section 3.3 for the procedure). We may make a comparison between the results from the PI approach with those obtained from MATLAB code based on SE 1  Again, the plots presented are not THE best results. They were generated for my choice of free parameters given in the starting of this chapter. In all the histograms I have chosen 50 bars and the bar width is automatically decided by MATLAB. Alternatively, bar width could be fixed rather than the no. of bars, for example 0.1, as in Ref. [3,12]. Next, we will explore the plots of the ground state probability density (see Section 3.3 for the procedure). We may make a comparison between the results from the PI approach with those obtained from MATLAB code based on SE 1  Again, the plots presented are not THE best results. They were generated for my choice of free parameters given in the starting of this chapter. In all the histograms I have chosen 50 bars and the bar width is automatically decided by MATLAB. Alternatively, bar width could be fixed rather than the no. of bars, for example 0.1, as in Ref. [3,12]. 1 This may be assumed to be the continuum limit to a very approximation.

First Excited State
The procedure to find the higher energy eigenvalue is outlined in Section 4.

First Excited State
The procedure to find the higher energy eigenvalue is outlined in Section 4.5. Using values from Tables 5.3  results 16,20 for the harmonic oscillator ( Fig. 6(a)) and cubic spline fits to the data points for the anharmonic systems ( Fig. 6(b,c,d)). The spline fits were used to estimate the ground-state energies in the continuum limit on linear axes; the logarithmic axis for δτ is used in Fig. 6 for presentation purposes only.
The ground-state energies for the systems shown in Fig. 6 have been calculated by several methods, with the results compiled in Table II The ground-state energy obtained from the extrapolation of the spline are within a few tenths Table II. Ground-state energy of the anharmonic oscillator for the indicated values ofλ. The column labelled MCMC is the data point in Fig. 6 with the finest discretization, Spline is the value obtained by extrapolating the spline curve, SE is the energy obtained by the numerical integration of Schrödinger's equation, and the last column contains the energies calculated by the method in Ref. 14

C. First excited state
The first excited of the anharmonic oscillator is obtained by evaluating the logarithmic derivative of the long-time limit of G(τ ) in (20), which is the zero-temperature limit of the correlation function in (19). Equation (21) then establishes the difference E 1 − E 0 between the ground state and first excited state as the negative of the slope of G(τ ) in the long-τ limit. On a time lattice, (21) is approximated as The approximate solution of this equation, is independent of τ . As ∆τ = nδτ for some non-negative integer n, we can determine E 1 − E 0 by plotting log[G 2,∞ (n)] versus n. An example is shown in Fig. 7 forλ = 0, 1, 50, and 10 3 with δτ = 0.2. The linear behavior, evident for small n, enables estimates to be made for E 1 , given the values of E 0 in Table II.
The calculation of the energy differences E 1 − E 0 are shown in Fig. 8 and the resulting energies of the first excited states are shown in Table III. The differences between estimates based on the extrapolation of the cubic spline and the exact values are of the order of 1% or less. However, the error bars for these calculations, particularly forλ = 0 andλ = 1, are much larger than the corresponding calculations of the ground-state energies (Fig. 6). Note, in particular, that despite

D. Second excited state
The determination of the second excited state from (24) proceeds in a manner similar to that in the preceding section. The approximate solution analogous to (34) is where G 4,∞ is the approximate solution of and is again independent of τ . The energy difference E 2 − E 0 is obtained by plotting log[G 4,∞ (n)] versus n.
The results of such an analysis are shown in Fig. 9. The larger error bars for all of the oscillators are clearly evident, as is the more limited range of discretizations that are practical. Nevertheless, when estimates for E 2 − E 0 obtained from the extrapolation of the spline fits in Fig. 9 are combined with the corresponding values of E 0 in Fig. 6 are compared with exact calculations, our estimates are found to agree to within a few per cent (Table IV).
We have applied the MCMC method to the evaluation of the imaginary-time path integral for the anharmonic oscillator with a coupling strengths ranging from the harmonic limit to strongly quartic. Quantities calculated include the probability density of the ground state and the energies of the ground state and the first two excited states. The ground-state probability density becomes more localized as the potential becomes more localized due to the increase in the quartic coupling constant, and the energy levels increase accordingly, as expected from elementary quantum mechanics. Information about higher lying excited states becomes increasingly difficult to extract because they are exponentially suppressed in the imaginary-time path integral. Indeed, a comparison of Figs. 6, 8, and 9 illustrates how the statistical fluctuations become more pronounced away from the ground state for the same number of paths used for measurements. We will present calculations of higher excited states in a separate publication.
Nevertheless, the imaginary-time path integral formulation has many attractive features for applications in condensed matter and other fields in physics. The method we have used here, which is described in detail in Ref. 16 and 22, can be applied to a particle in any potential. The inclusion of interactions between several particles is also within the scope of the method and is of interest even in one-dimensional systems. For example, the behavior of electrons in quantum wires and carbon nanotubes, where electronic motion is allowed in one dimension, but strongly restricted in the two lateral dimensions, is a central concern in theoretical and experimental studies of these systems. The harmonic potential can be used as a model of electrons (and holes) in semiconductor quantum dots.

Appendix A: The Partition Function
The partition function is Inserting the completeness relation for the orthonormal eigenkets ofĤ, 1 = ∞ n=0 |n n|, into the partition function yields which is the more familiar form of the partition function.

Appendix B: Correlation Functions
Consider the 2n-point correlation function: The first term is defined as Using the completeness relation for the orthonormal eigenkets ofĤ twice in the right-hand side, We also have Again using the completeness relation twice in the right-hand side, Since the expression for x n (τ ) does not depend on τ (which is expected) which is the same as (B5). Thus, by combining (B1), (B3), (B5) and (B6), we obtain To take the low-temperature limit of G 2n , we first have that 39 Then, for the terms in (B7), the limit β → ∞ yields Similarly, Hence, For n = 1, For n = 2, Consider the matrix element 0|x 2 (0)|q . For the harmonic oscillator, the ground state has even parity, the first excited state has odd parity, and parity alternates between even and odd for all higher lying states. We expect the same pattern for the anharmonic oscillator. Therefore, as the operatorx 2 and the ground-state wave function are both even functions of x, the matrix element Simulations that follow this approach are said to use the Metropolis-Hastings algorithm.

Basic Concepts and Notation
A Markov chain is a sequence of events in which a given position in the chain depends only on the immediately preceding position, rather than on the history of the chain. A Markov chain is defined by three attributes: 1. A state space which defines the values that can be taken by the chain.

2.
A transition operator that defines the probability of one state progressing to another state.
3. An initial condition that specifies the state from which the progression of the chain begins.
Our goal is to calculate the expectation value Ô of some operator O over a finite time interval, where is a configuration of the system, which also represents a path, and P eq (x k ) is the probability distribution in the canonical ensemble. The MCMC solution to this problem to the calculation of (C1) is to construct a Markov chain on the state space X that has P eq (x) as the stationary distribution. In other words, given transition probabilities W (x, x ) for chains x and x in state space, we havê The discretized path integral representation of the canonical partition function in (B8), identifies the equilibrium probability of the configuration x = {x 1 , x 2 , . . . , x Nτ } as and averages of an operator O(x) are calculated as in which [Dx(τ )] signifies that the integral is over all paths in the state space X.
The Monte Carlo method introduced by Metropolis is based on the idea of "importance sampling", whereby the phase-space points x in (C5) are not selected completely at random, which would be impractical for higher-dimensional spaces, but are chosen to be more densely distributed where n is the number of states generated by the Markov chain. Of course, in actual calculations, only a finite number of states are used to estimate averages, as determined by the operator and the system. The next section provides the mathematical foundation for (C6).

Properties of Markov Chains
If the state space is finite or countable, then the elements of the Markov chain can be represented as vectors and the transition probability as a matrix that operates on these vectors. However, most applications of the MCMC method have uncountable state spaces in which the initial state is an unconditional probability distribution and the transitions are expressed in terms of a conditional probability distribution. Finite state spaces are simpler to present so, we will use this setting in what follows.
For a Markov chain that consists of n states, the transition operator is an n × n matrix P whose entries p ij signify the probability of moving from state s i to state s j in a single step. Similarly, the probability to move from state s i to s j in 2 steps is n k=1 p ik p kj , which is the (i, j)th element of P 2 . The generalization to m steps is the (i, j)th element of P m .
If a transition operator does not change across transitions, the Markov chain is called (time) homogeneous. An important consequence of homogeneity is that, as t → ∞, the Markov chain will reach an equilibrium called the stationary distribution of the chain. The stationary distribution of a Markov chain is important for sampling from probability distributions, which is at the heart of MCMC methods. In more formal terms, suppose that the current state of the Markov chain is represented by the probability vector u which means that u i is the probability of being in state s i . Now we want to know the probability v j of finding the chain in state s j after m-steps. In matrix form, this is given by If there is a vector such that w T = w T P, which also implies w T = w T P m , then the Markov chain is said to be in equilibrium and w is called a stationary vector. Computationally, to find such a vector, one calculates the left eigenstate of the matrix P corresponding to unit eigenvalue.
The continuous state analogue of p ij is the transition density p(x, y). The continuous quantity corresponding to the m-step transition probability P m is denoted as p (m) (x, y), which is defined by the Chapman-Kolmogorov recursion relation, for n ≥ 2.
The most important theorem behind the MCMC method is: 35,36 Fundamental limit theorem. An irreducible, aperiodic Markov chain with transition matrix P has a stationary distribution w satisfying w j > 0, j w j = 1, and w T = w T P, if and only if all its states are positively recurrent, and w is unique and identical to the limiting distribution with the following definitions: • An irreducible Markov chain is the one in which the probability to go from any state s i to any other state s j in a finite number of steps is non-zero. Thus, there exists a finite n such that • The states of a Markov chain are called positive recurrent if the probability to return to the same state is unity after a finite number of evolution steps.
• The states are periodic if the transition s i → s i is only possible in steps which are integral multiples of the period d(i). The period is defined as the highest common factor for all m for which (P m ) ii > 0. For an aperiodic state s j , d(j) = 1.
Note that the stationary distribution obtained in (C9) is independent of initial state s i . This theorem gives us the freedom to start in any ergodic (aperiodic, irreducible and positive recurrent) Markov chain and we are guaranteed to end up in an equilibrium state. The process of evolution of the chain into a stationary state is called thermalization. Although stated in terms of discrete states, with suitable modifications, this theorem is also valid for continuous states. We refer the reader to Ref. 36 for details.

Appendix D: Splines
There are several methods for fitting curves to a set of given data points that enable the prediction of values between these points (interpolation) and the estimate of values outside the range of the points (extrapolation). In spline interpolation, piecewise functions pass through a set of data points, often referred to as knots, such that the function is smooth at these data points and satisfies some specified conditions at or near the first and last points. Thus, if there are n points (x i , y i ), for i = 1, 2, . . . , n, a total of n − 1 functions is used, one for each interval between the ith and (i + 1)st data points, for i = 1, 2, . . . , n − 1. Cubic polynomials are most commonly employed to strike a balance between smoothness and avoiding Runge's phenomenon, which is the appearance of oscillations at the edges of an interval that occurs when using high-order polynomials for a set of equi-distant interpolation points. The cubic polynomials are written as Each polynomial has four unknown coefficients , so the complete spline has 4(n − 1) unknowns. These unknowns must be chosen according to conditions that the spline must satisfy: 1. P i (x i ) = y i and P n−1 (x n ) = y n for i = 1, 2, . . . , n − 1. These conditions guarantee that the spline interpolates the data points.
2. P i (x i ) = P i+1 (x i ) for i = 1, 2, . . . , n − 1. These conditions require that the values of adjacent polynomials are the same at the points where they meet, which ensure that the interpolating spline is continuous.
These conditions supplement the continuity of the spline by requiring that the slopes of adjacent polynomials are the same at the point where they meet. Thus, the spline is differentiable on (x 1 , x n ).

4.
for i = 2, . . . , n − 1 or, in abbreviated notation, This condition supplements the continuity and differentiability of the spline by requiring that the second derivatives of adjacent polynomials are the same at the point where they meet.
Thus, the spline has a second derivative at every point on (x 1 , x n ).
The four conditions on the spline provide 1. The natural spline boundary conditions are that is, the second derivatives vanishes at the end points. These boundary conditions seldom used since they does not provide a sufficiently accurate approximation near the end points of the data set.
2. The clamped spline boundary conditions set the first derivatives at the end points to a particular value, which may be known, or specified at the user's discretion: The equations determining the coefficients for cubic splines for either natural or clamped boundary conditions can be expressed as tridiagonal matrices. At each point (i.e. knot), the spline changes from the cubic polynomial in the preceding interval changes smoothly to the polynomial in the next interval. The basic idea of the not-a-knot boundary conditions is to not change the cubic polynomials across the second and penultimate points, which are first two interior points. Thus, these two points are not knots, which results in the first two intervals having the same spline function and the last two intervals having the same two spline functions. The mathematical expression of these boundary conditions are The default in-built MATLAB function performs a cubic not-a-knot spline fit.