A multi-interval chebyshev collocation method for efficient high-accuracy rf circuit simula

A multi-interval Chebyshev collocation method for efficienthigh-accuracy RF circuit simulationBaolin Y angCadence Design Systems,Inc.,SanJose,CA95134byang@Joel PhillipsCadence Berkeley Laboratories,SanJose,CA95134jrp@ABSTRACTMost RF circuit analysis tools use either shooting-Newton or har-monic balance methods.Neither can efficiently achieve high ac-curacy on strongly nonlinear circuits possessing waveforms with rapid transitions.We present a multi-interval-Chebyshev(MIC) method that discretizes the circuit equations by dividing the simu-lation domain into a set of intervals whose size is adaptively chosen and using Chebyshev polynomials to represent the solution in each interval.The MIC method has excellent stability properties,is as effective at solving nonlinear problems as shooting techniques,can achieve high resolution on a wide variety of circuits,and in con-junction with an appropriate preconditioner can be combined with matrix-implicit Krylov-subspace solvers to analyze large circuits with moderate computational cost.1.INTRODUCTIONThe exploding demand for high performance wireless products has increased the need for efficient and accurate simulation techniques for RF integrated circuits.RF circuit simulation is difficult because RF circuits typically contain signals with multiple-timescale prop-erties,as usually the data and carrier signals in a system are sepa-rated in frequency by several orders of magnitude.Special-purpose RF simulators exploit the sparsity of the spectrum in order to make the computations tractable[1].For example,instead of performing a long transient analysis of a circuit driven by a periodic source, we may seek tofind the steady-state directly.In this paper we consider the central problem of periodic-steady-state computation. The periodic-steady-state problem is of interest in its own right, but is also important to study because it forms the basis of more advanced analyses,such as cyclostationary noise analysis,multi-frequency distortion analysis,or envelope simulation.The method presented in this paper can be advantageously applied to all the typical RF simulation problems.Two numerical methods commonly used in steady-state computa-tion are the shooting-Newton method,based on low-orderfinite To appear at37th Design Automation Conference,June2000.Not to be redistributed.difference discretizations such as the second-order Gear method, and the harmonic balance method,based on high-order spectral discretizations[1,2].The advantage of the low-order polynomial-based methods is that they can select time-points based on localized error estimates and as a result can easily handle sharp transitions in circuit waveforms.The harmonic balance method,on the other hand,has the advantage of attaining spectral accuracy for smooth waveforms.Recent developments in matrix-free Krylov-subspace algorithms[5,6]have made these methods even more popular as they can be used to analyze circuits with thousands of devices. High precision computations are often necessary in RF circuit sim-ulation.For example,accurately computing the noisefigure of a highly nonlinear RF circuit often requires accurate determination of the periodic operating point up to many multiples of the funda-mental frequency,as noise may be translated from very high fre-quencies by the mixing action of the time-varying circuit elements to appear at the output.High accuracy is easy to achieve for smooth waveforms using spec-tral approximation in the harmonic balance method.However,non-linear circuits often produce waveforms that have sharp-transition regions.The waveforms may only be in the case of a pulse, which means that the waveforms are continuous but the deriva-tives of the waveforms are not continuous at some points.In this case the spectral accuracy of the harmonic balance method is lost and acceptable accuracy is obtained only at the cost of including very many harmonics and/or timepoints in the simulation.To help handle sharp transitions,Nastov and White introduced the time-mapped harmonic balance(TMHB)method[4].In the TMHB method,a time-map function is used to map the solution in the original time space to a space where the solution is more smooth. This method is an improvement over traditional harmonic balance, but it is still restricted,since,theoretically speaking,a point in the original time space is still a point in the mapped time space,and so one still cannot achieve spectral accuracy.Basically it can only put more points around the points to localize the oscillations inevitably associated with using a spectral method to approximate functions.And,as with all methods that utilize global basis functions,for reasons of numerical stability there are also strong restrictions on how rapidly the timestep can be changed locally.The low-order multi-step methods,such as the Gear methods,typ-ically used in the shooting method offer excellentflexibility in lo-cally adapting the timesteps.For low to moderate accuracy,these methods are preferred for problems that are highly nonlinear and/or contain sharp transitions.However,low-order multi-step methodsPermission to make digital/hardcopy of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage, the copyright notice, the title of the publication and its date appear, and notice is given that copying is by permission of ACM, Inc. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee.DAC 2000, Los Angeles, California(c) 2000 ACM 1-58113-188-7/00/0006..$5.00are not efficient at high precision because veryfine discretizations are needed,resulting in loss of speed and increased memory re-quirements.In addition,the multistep discretizations used in theshooting method must be causal,and so have stability problems as-sociated with using one-side approximations,i.e.,performing diff-entiation using backward differences.It is generally accepted thatmultistep methods of order greater then three or four are not suf-ficiently stable,and in practice even the third order methods haverelatively stringent stability restrictions on how rapidly timestepsmay be varied.As a result,multi-step method of order higher than are not widely used in circuit simulation.Often,because of the backward-looking nature of the approximation,it is necessary todecrease the order after sharp-transition points or points,often tofirst order,resulting in veryfine grids in those regions.We seek a method that has none of the above disadvantages.Be-cause of inevitable discontinuities in the solution waveforms,resid-uals,or device equations,we perform the discretization by break-ing the simulation interval,of length,into timestep intervals of length,.Within each interval,Chebyshev poly-nomials are used to represent the solution waveform.Chebyshev polynomials have excellent approximation properties,in particular showing spectral convergence as the order is increased.We will show that this method has several desirable properties that lead to an excellent tradeoff between efficiency and accuracy,and further-more can be easily combined with matrix-implicit Krylov-subspace techniques to allow analysis of large circuits.In a sense we are able to exploit the best characteristics of the shooting-Newton and harmonic-balance methods,with the drawbacks of neither.2.CIRCUIT EQUATION DISCRETIZATION Circuit behavior is usually described by a set of nonlinear dif-ferential algebraic equations(DAEs)that can be written as(1)where represents contribution of the dynamic elements such as inductors and capacitors,is the vector representing static elements such as resistors,is the vector of inputs,and the vector contains state variables such as node voltages.The periodic steady-state solution is the solution of Eq.(1)that also satisfies the two-point boundary condition for some period.We now contrast three ways of discretizing the time-derivative op-erator:linear multistep methods,harmonic balance methods,and Chebyshev methods.As representative of linear multistep meth-ods we consider the Gear methods.Given the values,,at the time points:,to perform the discretization at ,we seek an interpolating polynomial of degree based on.The time derivative is approx-imated by evaluating the derivative of the interpolating polynomial at.Thus the time-derivative of the approximate solution is required to match the differential equation at the point.We call such a point a collocation point.The Gear methods are implicit in only“one-point-at-a-time”.They are called one-stage methods.At each timepoint,for a method of order,the solution vector at only points is involved,and if previous timepoints are known,a nonlinear system of size must be solved tofind the solution at the next timepoint.Each step of a shooting method involves the factorization of the circuit Jacobian matrix.Because this matrix is usually very sparse,it can be factored with fewfill-ins,in nearly time.The order of polynomial approximation is usuallyfixed at a low value,and increased accuracy is obtained by locally changing the timesteps based on estimates of the truncation errors induced by the order-polynomial approximation.Note that higher-order Gear methods can not be used at points after or sharp-transition points because approximation of the differential operator using an-order Gear method involves preceding points.Higher-order approximation across points or sharp-transition points leads to oscillation of the solution and makes the solution even less accurate than afirst-order solution.A similar problem occurs if the timestep ratios are changed too rapidly.For contrast,consider a particular(time-domain)construction of the harmonic balance method,known as Fourier-collocation.We seek to form the solution by an interpolating trigonometric polynomial of degree:(2)where denotes an-point interpolation operator.By comput-ing the derivative of the trigonometric polynomial(3)at the Fourier collocation points,, and enforcing the differential equation at these points we have spec-ified the numerical method.We may solve for the Fourier coeffi-cients directly,or express the Fourier coefficients in terms of the values at the Fourier collocation points as(4)The differentiation procedure may be performed efficiently to high order by noting that in the frequency domain,the differentiation op-erator becomes a multiplication operator,so by using the the FFT to convert the to,performing the differentiation,then in-verse transforming,we may accomplish the time-derivative oper-ation in operations.Alternatively,we can use this prescription to obtain a time-domain differentiation matrix such that approximates the value of at the collocation points, where vector is a vector of the values of at the collocation points.Then the derivatives can be obtained directly by a matrix-vector multiplication,requiring operations.Note that, in practice,using the FFT does not provide much advantage over using the matrix explicitly unless the order of approximationis fairly large.Accuracy in the harmonic balance method is usually achieved by increased the order of approximation.For a given set of colloca-tion points,the maximal possible order of approximation is used. With waveforms,this has the very desirable effect of achieving spectral convergence,meaning that as the spacing between the col-location points decreases,then asymptotically the error decreases faster than any polynomial approximation(in particular,faster than the Gear methods).However,note that though the expansion in Fourier series conveniently means automatically enforcing the pe-riodic boundary condition,it also means that the basis functions, sines and cosines,must cover the entire domain.Unlike the Gearmethods,which couple only local points,every point in the so-lution interval is coupled to every other through each basis function.The spectral convergence of the Fourier series is very desirable for smooth waveforms and nearly linear circuits,but un-fortunately,a single discontinuous point can destroy the spectral convergence properties.Because of the global nature of the basis functions,spectral convergence is very difficult to recover once in-terrupted.Most semiconductor device models possess some degree of discontinuity in the derivatives,and so true spectral convergence is rarely observed in practice.An alternative class of basis functions that also possess good ap-proximation properties,including the possibility of spectral conver-gence property,are the Chebyshev polynomials.However,Cheby-shev polynomials have the advantage that they can be applied when the simulation domain has been decomposed into multiple intervals in order to resolve complex or singular behavior,as has been done for electromagnetic wave problems in complex geometries[9]. Here we give the basics of the Chebyshev collocation method.Given the values,,at the Chebyshev Gauss-Lobatto collocation points,,we seek an interpolating polynomial of degree:(5) and compute the derivative of the polynomial as(6)The coefficients,,are computed from the coefficients by us-ing the backward recurrence relation[10](7)We can also obtain a derivative matrix by noting that the inter-polating polynomial can be expressed in terms of a Lagrange inter-polation polynomial.The entries of the Chebyshev differen-tiation matrix are computed by taking the analytical derivative of and evaluating it at the collocation points for,i.e.,.The entries of the matrix are(8)where when and otherwise.It is pointed out in[8]that calculating the matrix directly in this form results in unnecessarily large numerical error,and a formulation that yields more accuracy is given.The key to the multi-interval Chebyshev(MIC)approach is to break the simulation domain into multiple intervals,.In each interval the Cheby-shev collocation method is used to discretize the differential equa-tion.Continuity of the solution is enforced at the boundaries of any two intervals.For thefirst and last interval,this has the effect of enforcing the periodic boundary condition.The reason this de-composition is useful is that waveforms in circuit simulation typ-ically have regions of different properties.There are smooth re-gions where high-order methods can be used and there are sharp-transition regions where interval boundaries must be placed.Our particular combination of high-order Chebyshev discretizations in multiple time intervals provides several important advantages. First,as in harmonic balance,high accuracy is achievable because in each interval all the desirable approximation properties of Cheby-shev polynomials can be exploited.Chebyshev polynomials are usually close to the optimal polynomial approximant of a given order,and spectral convergence can be achieved if order-based(-type)refinement is used in each interval.Because of the rapid con-vergence,once a proper set of intervals is chosen high accuracy can be obtained at little marginal cost.Conversely,there will be few intervals in regions where the solution is smooth,because the intervals there can be quite large.Second,as in the lower-orderfinite difference schemes,the timestep intervals can be adaptively chosen to resolve nonlinear or abrupt transition behavior.The impact of nonlinear behavior can some-times be difficult to diagnose.Solution waveforms,such as volt-ages,may be quite smooth,but residuals,such as currents,will be much more difficult to represent with a Fourier series if the circuitequations are very nonlinear.In the MIC method,if the intervals are chosen such that discontinuities occur only at the interval bound-aries,then the discontinuities can be handled efficiently.Around a transition or strongly nonlinear region,more intervals can be placed in a denser timestep spacing,and possibly the order of the methodlowered.Third,the method has excellent numerical stability properties.It can be shown that our scheme is a member of a particular class of implicit Runge-Kutta(IRK)methods,the collocation-IRK meth-ods[3].Higher-order implicit Runge-Kutta methods are multi-stagemethods,that is,they are implicit in more than one point at a time. Because of the multiple implicit states,and unlike multi-step meth-ods such as the Gear methods,IRK methods can be constructed that are A-stable even at high order.By A-stable,we mean a scheme that is numerically stable for systems with poles located anywhere in the left-half-plane,including the entire imaginary axis.The par-ticular schemes we use in this paper are not A-stable,but we have been able to prove that the schemes of all orders used here are what is called-stable,meaning that if is a pole of the system be-ing integrated,then the numerical scheme is stable with timestep for all in the left-half-plane,except possibly within an an-gle of the imaginary axis.For our methods is usually less thanone degree,and the stability region only excludes a small semi-circle near the imaginary axis,which can be avoided by appropriate choice of timesteps.In addition,because we choose the Gauss-Lobatto-Chebyshev points for collocation,implying we always place a collocation point at the forward interval boundary,our particular construction gives methods that are L-or stiffly stable,just as the Gear methods are. Therefore the method is suitable for solving differential-algebraic equations as initial value problems.This property allows us to use shooting techniques to solve boundary value problems or construct preconditioners.In fact,in a certain sense that is important in practical applications,the MIC approach is more stable than the Gear methods.Since unlike the Gear methods the MIC scheme retains excellent stability properties even when the timesteps arerapidly varied,higher-order schemes can be used near sharp tran-sition points.Finally,the discretization scheme is robust.Unlike harmonic balance schemes,where sufficiently sharp transitions can induce oscillations in the global basis that are next to impossible to eliminate,in the worst possible case the discretization in a single MIC interval degenerates to the backward-Euler scheme which is robust and well-tested.In practice we have never found a need to drop below the second order scheme.3.EFFICIENT IMPLEMENTATIONFor simplicity of presentation,let us discuss Chebyshev-collocation methods for the case of two intervals:and .The number of collocation points in is and the numberof collocation points in is.The solution in is represented by and the solution in is represented by.The continuity condition is en-forced at the internal interval boundary,and the periodicity con-dition relates the two interval endpoints.Assume thatis the differentiation matrix in and is the differentiation matrix in ,Eq.(1)is discretized into the following equation:(9) where is the identity matrix.Newton’s method is used to solve this nonlinear equation.By defining the matrices(10) at each collocation point we can write the the Jacobian matrix for the Newton iteration asNote that and.At each Newton iteration step,we must solve a linear system(11) where is the residual of the previous solution.This step is prob-lematic because in each interval the discretization of the differentia-tion operator is global.That is,the approximation of at every collocation point involves at all the collocation points. Thus the differentiation matrix,like the Fourier collocation differ-entiation matrix,is a full matrix(on each interval).The resulting larger,less-sparse equations systems lead to difficulties in large-scale circuit simulation and is the reason IRK methods are notwidely used to date.Ironically,it is the sheer scale of RF simu-lation problems that makes IRK methods more suitable in the RFcontext.Matrix-implicit Krylov subspace methods,such as GM-RES or QMR,are now commonly used to simultaneously solvethe linear system of equations(Eq.11)resulting from discretiza-tion of the circuit equations over a fundamental period.We nowshow that it is possible to easily construct a preconditioner that ismore effective than those typically used in harmonic balance andcomparable to the preconditioners used for shooting methods,thusassuring rapid convergence of the iterative solvers.For example,ona similar test(the DC-DC converter),with a comparable approach to preconditioner construction,our method requires two orders ofmagnitude fewer GMRES iterations than reported in[4].In each interval we use a“local”preconditioning approach thatis motivated by techniques used for large-scale harmonic balanceanalysis.But because we can decompose into multiple in-tervals,we are able to devise a preconditioner that is much moreefficient than the corresponding harmonic balance preconditioner.In this way,our“global”preconditioning approach is similar tothat used in[5]for shooting methods,where a preconditioner isobtained forfinite difference equations by dropping the top-rightcorner entries in the left-hand-side matrix of the linear system.In harmonic balance,it is common to construct a preconditioner by assuming the capacitance and conductance matrices respec-tively can be well represented by(piecewise)-constant approxima-tions.We adopt both these techniques,for example using a single as an approximation of all the capacitance matrices in interval and as an approximation of all the capacitance matrices in interval,the conductance matrices being treated similarly. When this is done,note that the two large blocks of the precondi-tioning matrix corresponding to the two intervals are only coupled by the subdiagonal block.Hence it suffices to discuss how to in-vert the top diagonal block.To accomplish this,consider forming the Schur decomposition of the relevant portion of the Chebyshev differentiation matrix such that(12) where is a Schur matrix and is a unitary matrix.has2-by-2blocks on the diagonal and all other subdiagonal entries are zeros.By applying this decomposition to the top block of,we can reduce it to an almost upper block-triangular structure with2N-by-2N blocks on the diagonal.These2N-by-2N blocks on the diagonal can be further transformed into a N-by-N complex matrix that has the same sparsity structure as a single-timepoint circuit Jacobian as generated in a transient analysis.It can be factored efficiently,with fewfill-ins,just as the Jacobian matrix in a low-order method. The advantage of this preconditioner over harmonic balance pre-conditioners is that we can approximate the capacitance and con-ductance matrices with averages separately in each interval at lit-tle computational cost,in contrast to using global averages which make the preconditioner very inefficient in simulations of nonlinear circuits.Because the intervals are chosen to resolve the nonlinear behavior,we can expect the cost of applying the preconditioner to be relatively independent of the degree of nonlinearity of the cir-cuit.If the circuit is very nonlinear,there will be very many in-tervals,yet precisely for this reason the preconditioner will still be a good approximation to the original matrix,and in fact will be nearly as cheap to apply as the shooting-method preconditioner.In a practical implementation it is not advantageous to attempt to solve the MIC equations directly.Instead,we use a low-order shooting-Newton method to obtain an initial approximation to the solution waveform.For steady-state problems,the shooting-Newton method has a good global convergence rate in PSS simulations and can rapidly reach the neighborhood of the steady state.From the result of the low-order solution method,we decompose the time interval into a number of subintervals.The objec-tive of thisfirst decomposition is to avoid sharp transitions and points in each interval.Next we determine the order of accuracy desired in each interval.At this point,if a interval contains many time points of the low-order shooting method,the solution in this interval must be relatively smooth and we will choose a relatively high order in this interval.If a interval contains only a few time points,it means that we are in a sharp-transition region and we can choose low-order IRK methods in the interval.Once we have dis-cretized the circuit equations in all the intervals and the boundaries are connected with continuity conditions and the periodicity condi-tions then we can proceed with the Newton iteration.Once an initial solution is obtained,we use an adaptive refinement scheme to split the intervals in any additional sharp-transition re-gions that may appear and increase the order of approximation in other regions where accuracy is not satisfactory.In this stage, we use both h-type(interval decomposition)and p-type(order in-crease)refinement to achieve the required accuracy in the most ef-ficient way.Error control can be performed by examining the trailing Cheby-shev coefficients as indicators of potential numerical truncation er-ror.Our initial experiments have been conducted using the follow-ing simple error criteria,(13)where is the order of the Chebyshev method used in interval, and are the th and th Chebyshev expansion coefficients,normalized by the local maximum in all intervals,is the length of the interval,is the time period,and is the desired level of accuracy.This takes into account both h-type refinement, ,and the p-type refinement,,which is proven to be effective in numerical experiments.If the error criteria is not satisfied in a given interval,then we check the absolute value of and.If,we increase the order of the method in the interval,as this indicates that the order of Chebyshev coefficients is decreasing,so it is reasonable to seek spectral convergence.Otherwise,the solution is likely non-smooth, and we decompose the interval into two new intervals at the middle point.4.NUMERICAL RESULTSIn this Section we show the spectral accuracy of the multidomain Chebyshev collocation method in periodic steady state computa-tions.Efficiency of using h-type(domain decomposition)and p-type(polynomial order increase)refinement in resolving waveforms with sharp transitions is shown.Our method’s ability to resolve sharp transitions with high accuracy is verified.We also demon-strate the efficiency of the preconditioning technique in the solu-tion of the linear system using the matrix-free Krylov-subspace ap-proach.In thefirst example,we show the spectral accuracy of the multido-Figure1:Spectral accuracy of multidomain Chebyshev collo-Figure2:Comparison of h-type,p-type,and hp-type refine-main Chebyshev collocation methods by computing the periodicsteady-state solution of a crystalfilter.It can be seen in Fig.1that our method can achieve machine precision using only time points(we use intervals and16th order in each interval).Spectralconvergence is evident in the exponential decay of the error.In the second example,we show the efficiency of using both h-typeand p-type refinement in resolving solutions with sharp transitions.We introduce nonlinearity into the previousfilter example by plac-ing mixer before it.We compare three refinement strategies:h-type refinement only,p-type refinement only,and hp-type(h-typeand p-type)refinement.In Fig.2,one can notice that both p-typerefinement and hp-type refinement can achieve spectral accuracy.In the hp-type refinement,we use Chebyshev polynomials with or-der varying from to in the intervals and we achieve the same level of accuracy that is achieved using p-type refinement wherethe order as high as in used in two intervals.Note that for thesame total number of time steps,the hp-type refinement is moreefficient.We can clearly see from Fig.2that both hp-type refine-ment and p-type refinement is superior to h-type refinement of a second order method.The advantages of using high-order methods are clear from this example.In the last example,we calculate the periodic steady-state of a verynonlinear DC-DC converter similar to a circuit analyzed in[4].Thefinal mesh of the proposed approach is obtained through three steps of hp-refinement.By using both h-type and p-type refinement,the。