BACKGROUND OF THE INVENTION Using a description language such as a netlist and device models an analog circuit can be first designed in terms of its predetermined inputs and expected outputs. The analog circuit design is then simulated before it is physically fabricated on a silicon chip.

One of the most difficult challenges in analog circuit simulation is the analysis of the circuits that operate on multiple time scales. Typical examples of this type of circuits are switched-capacitor filters and circuits used in RF (radio frequency) communications systems.

Applying standard transient analysis to a circuit of this type requires simulation of the detailed responses of the circuit over hundreds of thousands of clock cycles (millions of time points).

Many circuits of engineering interest are designed to operate near a time-varying, but quasi-periodic, operating point. Some of these circuits can be analyzed under the assumption that one of the circuit inputs produces a periodic response that can be directly calculated by steady-state algorithms, thus avoiding long transient simulation times. Under this assumption, all other time-varying circuit inputs are treated as small-signal by linearizing the circuit around the periodic operating point.

Existing algorithms are able to find periodic operating points and to perform periodic time-varying small-signal analysis. However, many circuits cannot be analyzed with the periodic-operating-point-plus-small-signal approach, because the above-described assumption may not apply. For example, predicting intermodulation distortion of a narrowband circuit, such as a mixer-plus-filter circuit, involves calculating the nonlinear response of the mixer circuit, driven by an LO (local oscillator), to two high-frequency inputs that are closely spaced in frequency. The steady-state response of such a circuit is quasi-periodic.

The analog circuit simulation is further complicated by the fact that many multi-timescale circuits have a response (again mixers and switched-capacitor filters are typical examples) that is highly nonlinear with respect to at least one of the exciting inputs, and so steady-state

approaches, such as the multi-frequency harmonic balance approach, do not perform well. To circumvent these difficulties, mixed frequency-time (MFT) algorithms have been proposed.

Specifically, the MFT algorithms exploit the fact that many circuits of engineering interest have a strongly nonlinear response to only one input, such as the clock in the case of a switched- capacitor circuit, or local oscillator in the case of a mixer, but respond only in a weakly nonlinear manner to other inputs.

Unfortunately, existing MFT algorithms suffer from several drawbacks that prevent their application to practical circuits, particularly large circuits. In existing MFT algorithms, poor sample point selection leads to ill-conditioned simulation environment, in which simulation values may be unsolveable with acceptable accuracy. In addition, existing MFT algorithms are based on a matrix-explicit linear solver (via Gaussian elimination) whose computational cost (or time) is proportional to an order of N3 for each Newton iteration, where N is the number of nodes of the circuit in simulation.

A new class of algorithms has been developed for simulating multi-timescale circuits by converting the circuit DAE (differential-algEbraic equation) into an equivalent multi-variable partial differential equations (M-PDE). However, the effectiveness of the M-PDE method to simulate large circuits has yet to be proven. In addition, there is evidence that, for some circuits, the M-PDE method generates inaccurate simulation results.

There is, therefore, a need in the art for a method and apparatus that utilizes the MFT method to accurately simulate large circuits.

There is another need in the art for a method and apparatus utilizing the MFT method to simulate large circuits with reduced computational cost and increased speed.

There is still another need in the art for a method and apparatus for generating an efficient linear problem solver structured such that the MFT method can accurately simulate large circuits with improved convergence.

The present invention provides a method and apparatus to meet these and other needs.

SUMMARY OF THE INVENTION To overcome the shortcomings of the available art, the present invention discloses a novel method and apparatus for simulating analog circuits by using a mixed frequency/time approach.

In broad terms, the present invention provides a method for simulating responses of a circuit, the circuit receiving a periodic sample signal and at least one information signal. The method comprises the steps of : selecting a set of distinct time points; defining a set of reference

time points, wherein each of the reference time points is associated with one of the distinct time points; establishing a first set of relationships between the values at the distinct time points and the values at the reference time points; establishing a second set of relationships between the values at the distinct time points and the values at the reference points; combining the first and second relationships to establish a system of equations in terms of the values at the distinct time points; and finding responses of the circuit at the distinct time points by solving the established system of equations.

The present invention also provides a corresponding apparatus for performing steps in the method described above.

BRIEF DESCRIPTION OF THE DRAWING The above mentioned advantages of the present invention as well as additional advantages will be more clearly understood as a result of a detailed description of the preferred embodiments of the invention when taken in conjunction with the following drawing in which: Figure 1 shows circuit 102, which will be used to illustrate a simulation process in accordance with the present invention; Figure 2 shows an exemplary envelope in response to two inputs shown in Figure 1 ; Figure 3 shows a scheme of selecting K (K = 2) sample points along the sample envelope shown in Figure 2 to obtain five distinct time points, in accordance with the present invention ; Figure 4 shows the effectiveness of the preconditioning process in compressing the eigenvalues for an RF receiver circuit, in accordance with the present invention; Figure 5 shows the effectiveness of the preconditioning process in reducing the number of iterations needed to solve each MFT Newton update equation, in accordance with the present invention; Figure 6 shows a simulation result for a filter circuit in accordance with the present invention; Figure 7 shows the intermodulation distortion of a high-performance receiver; Figure 8 shows a flowchart illustrating a process of simulation the responses of a circuit, in accordance with the present invention; Figure 9 shows a flowchart illustrating an exemplary process of solving the system of equations established in Figure 9, in accordance with the present invention ; and Figure 10 shows block diagram illustrating an exemplary computer system, which can be used as a hardware platform for executing the program that performs the processes shown in Figures 8 and 9, in accordance with the present invention.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS Figure 1 shows circuit block 102 including N circuit nodes, which will be used to illustrate a simulation process in accordance with the present invention. Circuit 102 receives S information signals 11, 12,..., Is from inputs 104.1,104.2,..., 104. s, respectively. Circuit 102 also receives a clock (or reference) signal C from input 106. In response to receipt of the input signals, circuit 102 generates a signal at output 108.

1. The improved MFT algorithm The behavior of a circuit (such as circuit 102 in Figure 1) can be described by a set of nonlinear differential-algebraic equations (DAEs) that can be written as d Q(v(t))+I(v(t))+u(t) = 0, Where Q (v (t)) E N 1S typically the vector of sums of capacitor charges at each node, I (v (t)) E stN is the vector of sums of resistive currents at each node, u (t) E flN is the vector of inputs, v (t) e SN iS the vector of node voltages, and N is the number of circuit nodes.

The present invention is particularly advantageous in situations where the input signal u (t) is quasiperiodic. A signal is L-quasiperiodic if it can be written as a Fourier series with L fundamental frequencies. RF circuits are generally influenced by one periodic timing signal, often referred to as the LO (local oscillator) or the clock, and one or more information signals. If denotes the clock signal frequency (such as the signal at input 106 shown in Figure 1), and f.. f are the S information signals (shown in Figure 1), then the (S + l)-quasiperiodic input can be written as: In a preferred embodiment, the present invention utilizes two conditions to improve the simulation of quasi-periodic circuit operating conditions. The first condition is that the circuit of interest possesses a quasiperiodic steady-state response. That is, v (t) is an S + 1 quasiperiodic signal with fundamentals fi,..'f, f.. The second condition is that all physical circuits have a finite bandwidth. Using these two conditions, the present invention selects only a finite number

of Fourier series terms to approximate v (t) while maintaining the necessary accuracy. Thus : where V (kl,..., k,, kd E=-e. (An interesting property of the MFT algorithm is that it is not necessary to truncate to a finite number of harmonics of fc.) Assume that v (t) is sampled at a discrete set of points t', z = to + nTc, where Tc = lut is the clock period, to c [0, Tc) and n runs over the integers, to obtain a discrete signal ' (t). Since the"envelope"v (t') is S-quasiperiodic and can be represented as a Fourier series in only the "information"fundamentals. The clock fundamental has disappeared, thereby forming an envelope.

Figure 2 shows an exemplary envelop in response to the signals at input 104.1 (assuming that the inputs at 104.2,..., 104. S are zero inputs) and input 106. In Figure 2, the hollow, circular dots represent samples (or distinct points), while the continuous waveform (dashed line) is the waveform having V as its Fourier coefficients, or equivalently, obtained by Fourier interpolation of the sample points.

In principle, because there are only K (2K3 + 1) Fourier coefficients to represent v, once the value of K distinct points tl,..., tk along the sample envelope are known, then the full envelope can be recovered. The envelope corresponding to the quasiperiodic operating point is obtained by obtaining K sample points that lie on the solution to the DAE given by equation (1).

Figure 3 shows a particular scheme of selecting K (K=2) sample points along the sample envelope shown in Figure 3 to obtain five distinct points (or distinct time points) [5 = (K x 2) + 1], in accordance with the present invention. In Figure 3, the distinct points are denoted by circles, and the reference points (reference time points) are denoted by squares.

Let us define the state transition function # (v0, tk, tf) = v (tf) = V (tf) : v (t) that satisfies equation (1) for t # [tk,tf] and v (tk) = vo In particular, define the vector #0 = [#T(t1), ..., #T(tK)]T = [vT(t1), ..., vT(tK)]T, (6 where superscript T denotes matrix transpose, to contain v at the K sample points t7,, = t0 + nkTc, k = 1... K, nk E Z. The value of the K points that follow by one cycle can be obtained from the transition function, 10#TcT = [vT(t1+Tc), ..., vT(tK + Tc)]T = [#(v(t1), t1, t1 + Tc)T, ..., #(v(tK), tK +<BR> <BR> <BR> <BR> <BR> (7)

which may be written more compactly by introducing the multi-cycle transition function that is the collection of the K transition functions from tk to tk as VT, Tc(#0). (8) Note that for each mode n, the vector of signals on that node, at the sample time plus one clock cycle, #Tn, is a delayed version of the signals at the sample points (this will be discussed in greater detail below). There exists, therefore, a linear operator D that maps v0 to VT <BR> <BR> <BR> <BR> #Tcn = DTc#0n.<BR> <P> (9) Note that Dy is a real matrix and independent of node n. Hence equation (9) holds for each n = 1,..., N, and represents a boundary condition on a solution to (1).

Combining equations (8) and (9) gives (DTc # IN)#0 - #Tc(#0) = 0, (10)

where 0 is the Kronecker product and IN is the N by N identity matrix. Equation (10) is a system of KN nonlinear equations and KN unknowns that can be solved for the envelope sample points. From these sample points and the transition functions, the circuit's quasiperiodic operating point (in particular, the spectrum of v) can be recovered.

2. Sample point selection To construct the matrix Dito, referred to as the delay matrix, consider the Fourier series of vo and Tc. Referring to equation (4), equation (11) holds where Tc(k1, ..., ks) = e-j2#k1f1Tc ... e-j2#ksfsTc Thus if is the matrix mapping sample points on the envelope to Fourier coefficients, then the delay matrix may be constructed as DTc = #-1#Tc#. ( In particular r may be constructed as the Kronecker product of one-dimensional (2Ks+l)-point Fourier-transform matrices (s) = ej2#mf stn/(2 Ks+1) mn as r = r'reS) (15) From the properties of Kronecker products, r-1 is likewise a Kronecker product of the inverses of the rocs'. In the existing MFT algorithms, no particular consideration was given to the choice of the sample points tk, so that the r6S)'s there are ill-conditioned matrices corresponding to an "almost-periodic"Fourier transform. By contrast, the improved MFT algorithm of the present invention performs a process of choosing well-conditioned sample points.

Assume the K sample points can be arranged into an S-dimensional array i (kl,..., ks),- such that for a given dimension s, there exists an integer p, and <BR> <BR> <BR> <BR> <BR> #(###,ks + 1,###) - #(###,ks,###) =<BR> <BR> 2A +1 holds. In this case, the entries of the r (S) matrices are: <BR> <BR> <BR> <BR> <BR> ) = ej2#mn/(2Ks+1)<BR> inn That is, they are the DFT matrices, and the matrix r :C2K1+1 # ### # C2Ks+1 #C2K1+1 # ### # C2Ks+1 represents an S-dimensional DFT. Thus r has a condition number of one; it is perfectly well-conditioned.

3. Matrix-implicit solution procedure The Newton's method can now be employed to solve (DTc # IN)#0 - #Tc(#0) = 0, (18) At iteration i, the Jacobian matrix is given by Recall from (13) D = r-lQT, which is fixed through all Newton iterations. Let J = ####0 = #0i be obtained from the multicycle transition function by e

Note that Jis block-diagonal. Defining b = -(DTc#IN)#0i - #(#0i), the Newton iteration is performed by solving the equation ((DTc # IN) - J)#(#0i) = b, using an iterative Generalized Minimal Residual (GMRES) solver, and setting #0i+1 = #0i + ##0i. (22) Each iteration of GMRES requires a matrix-vector multiplication. For a vector q E f N, the term (DTe 0 Iv) q is calculated by first applying a K dimensional DFT N times, then scaling each row with C2T., and finally applying a K dimensional inverse DFT N times.

Let q be partitioned into q = [ql-..., qKT]T,qk # RN, for 1 # k # K. Then The calculation of each 8 V-+§t) qk can be carried through matrix-vector multiplication and backsolving without explicitly forming the matrix.

4. Preconditioning For many problems, the GMRES algorithm is not efficient for solving equation (21) without an effective preconditioner. To analyze the reason, consider the case where the state transition function of the circuit, over one clock cycle, is approximately linear, that is # (X, t, t + e) = HX (t). Linear circuits are an obvious example of a case where this is true, and while nonlinear circuits will have nonlinear state-transition functions, if the method performs poorly for linear circuits it surely will not work well for nonlinear circuits either. However, many nonlinear circuits have a state-transition function that is nearly linear, a fact which is exploited below to construct an effective preconditioner. The convergence of the GMRES method will depend on the location of the eigenvalues of the Jacobian matrix, DT-J. IfXHisan

eigenvalue of the matrix H, then ei#tc-#H, where # = 2#(k1f1 + k2f2+###ksfs) will be an eigenvalue ofDy-J, for every kl, k2, ### in the MFT analysis. Thus unless all the secondary input frequencies are nearly commensurate with the clock frequency, the eigenvalues of DRY-J will be"fanned out"by delay matrix. This will cause severe convergence problems for the GMRES solver. Roughly speaking, the GMRES algorithm in the MFT algorithm with K total harmonics will take K times as many iterations to coverage than the GMRES iteration for the steady-state problem with only the clock excitation applied. This follows because the eigenvalues of H are typically inside the unit circle of the complex plane. The delay matrix replicates the eigenvalue structure K times, each shift being a complex number of order unity, and generally causing the convex hull of the eigenvalues of-to enclose the origin.

The following lemmas about the properties of Kronecker products are needed to perform the formal analysis.

Lemma 5.1 If A1, A2,...,Ap # Fmxm, B1,B2,...,Bp # Fnxn then A1A2...ap)#(B1B2...Bp) = (A, # B1)(A2 # B2)...(Ap # Bp).

Lemma 5.2 If A # Fmxm, B # Fnxn then (a) (A # In)(Im#B) = (Im # B) (A @ In).

(b) (A#B)-1 = (A-1#B-1).

Theorem 5.3 If xi is an eigenvalue of a vO, then ei25 «'kTe is an eigenvalue of the MFT Jacobian matrix.

The proof is as follows. For linear circuits, the diagonal blocks of,)) are the same, i. e., <BR> <BR> <BR> <BR> <BR> 0<BR> .-, y Denote a diagonal block as H, then the Jacobian matrix is equal to<BR> <BR> vol<BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> (#-1#Tc#)#IN-(IK#H) (24)<BR> <BR> <BR> <BR> <BR> = (#-1#IN)(#Tc # IN)(##IN)-(IK#H) (25)<BR> <BR> <BR> <BR> <BR> = (#-1#IN){(#Tc#IN)-(#-1#IN)-1(IK#H)(##IN)-1}(##IN) (26)<BR> <BR> <BR> <BR> <BR> = (#-1#IN){(#Tc#IN)-(##IN)(IK#H)(##IN)-1}(##IN) (27)<BR> <BR> <BR> <BR> <BR> <BR> = (#-1#IN){(#Tc#IN)-(IK#H)(##IN)(##IN)-1}(##IN) (28) = (r-' sIN) {(#Tc#IN)-(IK#H)}{##IN). (29)

Equation (24) to equation (25) holds because of IN=INININ and Lemma 1. Equation (26) to equation (27) holds due to Lemma 2 (b), and equation (27) to equation (28) holds due to Lemma 5.2 (a). Since (#-1 # IN) is unitary and its inverse is (r # IN)-', the right hand side of equation (29) has the same spectrum as (QTe 0 IN)- (IK 0 H). It is easy to verify that (QTC #IN)-(IK#H) is block diagonal, hence its eigenvalues are the union of eigenvalues of all the blocks, e 2Zr§9kT IN-H, for k = 1,---, K.

The preceding analysis suggests a good way of preconditioning for solving the Newton equation (21). Solving equation (21) is equivalent to solving {QTC- #IN)J#(#-1#IN))}α = rb, (3 0) wherey=rA-vi. A good choice of preconditioner is P=(#Tc#IN)-(IK#H) where Hcan be chosen as the Jacobian matrix from the steady-state analysis in the initial guess stage discussed in Section 5, or any of the diagonal blocks, ######, for i = 1, ..., K, of ###. In particular, if the single-cycle state-transition function is linear and time invariant, then the Newton equation can be solved in a single GMRES iteration. Note that the preconditioner presented here is effective if the Jacobian of the state-transition function is nearly constant over multiple cycles. The circuit behavior inside each clock cycle is hidden from the preconditioner. This is not the case in, for example, the time-or frequency-averaged preconditioners typically used in modern harmonic balance codes. For this reason the preconditioner presented here may perform well under much weaker assumptions about the circuit behavior, in particular at higher power levels.

For each GMRES iteration, a system Pu = v has to be solved. Since P is block diagonal, it needs to solve a sequence of K systems (ej2##kTc IN-H)uk=vk, for k = 1,..., K, where uT=[u1T,...,uKT]T and vT=[v1T,...,vKT]T. The preconditioner can be applied very efficiently by incorporating a Krylov subspace reuse algorithm, as the linear systems to be solved are the same as arise in the small-signal analysis for periodic time-varying systems. The basic idea of the algorithm is that the Krylov subspaces associated with the matrices eiT-H are very similar for different C9k. Essentially, the Krylov-subspace re-use algorithm allows the preconditioner for the matrix (DTC'5i 0 IN)-Jto be applied with only slightly more cost than an iterative solve with the matrix H.

Figure 4 shows the effectiveness of the preconditioning process in compressing the eigenvalues for an RF receiver circuit. Specifically, Figure 4 shows eigenvalue distribution before and after the preconditioning process. hi figure 4, the eigenvalues are very tightly

clustered around unity, indicating excellent performance of the preconditioner and very rapid GMRES convergence.

Figure 5 shows the effectiveness of the preconditioning process in reducing the number of GMRES iterations needed to solve each MT Newton update equation, for the same RF circuit mentioned above. Only three iterations are needed to reduce the residual by a factor of 10 2, whereas without the preconditioning process, over 400 iterations are needed to achieve any reduction in the residual at all. Since the MFT circuit equations are not solved exactly, on average there is a performance advantage in the MFT method to using approximate solves of the Newton update equation, and therefore GMRES is converged to a relatively loose tolerance.

5. Improving Newton convergence Rapid convergence of Newton's method can only be assured with a good initial estimation. To achieve a good initial estimation, the present invention first calculates the periodic steady state response of the circuit with the clock signal applied, while suppressing other non-DC signals. Using the steady state solution as an operating point, a small-signal analysis is performed by treating non-clock fundamentals as small signals. As a result of the small signal analysis, amplitudes atfs + ksfc, for -Ks#Ks#Ks, 1#s#S, are generated. These amplitudes are transformed into time domain initial conditions via inverse multidimensional discrete Fourier transform (DFT). At higher input power levels, using a Newton continuation method, with the amplitude of the non-clock signals as the continuation parameter, is generally effective in securing convergence.

6. Spectrum calculation After the solution is converged, the values = [v (tl) T, v (t2), V (tK) T] T and the integration solution in # = [v(t1)T, v(t2)T,###,v(tK)T]T are available. From these pieces of information, the spectrum v (t) can be obtained. Let Define v (z) = [v (tl f T), V (t2 + T) T,..., v (tK + z) T) T. Tlien

Then for each KN-vector V (', kC) where-K < kc < KC which is collection of all N-vectors V (kl,..., ks, k^), where -K1#k1#K1,...,-Ks#ks#Ks (the actual order is determined by the Fourier transform), Forming {(#(#)-1#)#IN}#(#) requires the values for v (tel-),..., v (tact), or synchronized time steps between cycles. The total computational costs is one KN-vector integration and MFourier transforms, where M is the number of synchronized time points.

The synchronized time step requirement may not be easily met in practice. One alternative is to use interpolation schemes. However, these schemes potentially lose accuracy.

Another alternative is to utilize integration instead of multidimensional discrete Fourier transforms. Specifically, it is easy to verify that where Ep is a KNxN block matrix whosepth Nx N block is IN and other blocks zero, and p is determined by (kl,..., ks) from the Fourier transform. Calculating equation (34) does not require synchronized time points. The total cost of calculating V (., kc) is K KN-vector integrations plus one final Fourier transform. However, it might be more expensive since integrations normally cost more than Fourier transforms.

7. Simulation Utilizing A Preferred Embodiment Of MFT Method The first example is a low-pass switched-capacitor filter of 4kHz bandwidth and having 238 nodes, resulting in 337 equations. To analyze this circuit, the improved MFT of the present invention analysis was performed with an 8-phase 100kHz clock and a IF sinusoidal input at lOOHz.

The 1000 to 1 clock to signal ratio makes this circuit difficult for traditional circuit simulators to analyze. In the improved MFT method, three harmonics were used to model the input signal. The eight-phase clock resulted in the need to use about 1250 timepoints in each transient integration. This brings the total number of variables solved by the analysis to slightly less than three million (337 x (2 x 3 + 1) x 1250 = 2,948,750). The simulation took a little less than 20 minutes CPU time to finish, on a Sun UltraSparcl workstation with 128 Megabyte memory and a 167MHz CPU clock. Figure 6 shows the output spectrum of the filter.

The second example is a high-performance image rejection receiver. It consists of a low- noise amplifier, a splitting network, two double-balanced mixers, and two broadband Hilbert transform output filters combined with a summing network that is used to suppress the undesired side-band. A limiter in the LO path is used for controlling the amplitude of the LO. It is a rather large RF circuit that contains 167 bipolar transistors and uses 378 nodes. This circuit generates 987 equations in the simulator.

To determine the intermodulation distortion characteristics, the circuit was driven by a 780MHz LO and two 50mV closely placed RF inputs, at 840MHz and 840MHz+10KHz, respectively. Three harmonics wee used to model each of the RF signals. 200 time points were used in each transient clock-cycle integration, considered to be conservative in terms of accuracy for this circuit. As a result, nearly ten million unknowns (987 x (2 x 3 + 1) 2 x 200 = 9,672,600) were generated. It took 55 CPU minutes to finish on a Sun UltraSparelO workstation with 128 Megabytes of physical memory and a 300MHz CPU clock. Figure 7 shows 3rd and 5th order distortion products.

To understand the efficiency of the improved MFT method of the present invention, consider that traditional transient analysis would need at least 80,000 cycles of the LO to compute the distortion, a simulation time of over two days. In contrast, the MFT method of the present invention is able to resolve very small signal levels, such as the 5th order distortion products show in Figure 7.

Solving the MFT equations by direct factorization methods is also impractical, as the storage needed for the factored rank-50,000 (987 x (2 x 3 + 1) 2 = 48,363) MFT Jacobian of Equation 19 is several gigabytes. Forming the Jacobian matrix by direct methods would also require computation time proportional to the cost of 50,000 transient integration cycles, again a number on the order of days.

Figure 8 shows a flowchart illustrating a process of simulating the responses of a circuit such as the circuit shown in Figure 1, in accordance with the a preferred embodiment of the present invention.

Step 804 selects a set of evenly spaced distinct time points shown as the circle dots in Figure 3. The details of step 804 can be found in equation (16) and related descriptions.

Step 806 defines a set of reference time points shown as the square dots in Figure 3. As shown in Figure 3, each of the reference time points is associated with a respective distinct time point, and each of the reference time points is one signal period (or clock cycle) away from its respective distinct time point.

Step 808 establishes a first set of relationships between the values at the distinct time points and the values at the reference time points. The details of step 808 can be found in equation (8) and related descriptions.

Step 810 establishes a second set of relationships between the values at the distinct time points and the values at the reference time points. The details of step 810 can be found in equation (9) and related descriptions.

Step 812 combines the first and second sets of relationships to establish a system and equations that contain the values at the distinct time points only. The details of step 808 can be found in equations (10) and (18) and related descriptions.

Step 814 finds (or generates) the simulated responses of the circuit at the distinct time points by solving the established system of equations. If a circuit includes N internal circuit nodes and M outputs, step 814 can End (or generate) the simulated responses for all of the N internal circuit nodes and M outputs. The details of step 814 can be found in equations (18)- (22) and (30) and related descriptions.

Figure 9 shows a flowchart illustrating an exemplary process of solving the established system of equations in step 814 of Figure 8, in accordance with the present invention.

Step 904 selects a set of estimated values to reflect estimated circuit responses at the distinct time points. The details of step 904 can be found in Section 5.

Step 906 establishes a system of linear equations at the estimated values. The details of step 906 can be found in equation (21) and related descriptions.

Step 908 preconditions the system of linear equations to improve the convergence of solution to the system of linear equations. The details of step 908 can be found Section 4.

Step 910 solves the system of linear equations to generate the correction values to adjust the estimated circuit responses at the distinct time points. The details of step 910 can be found in equations (21)- (22) and related descriptions.

Step 912 adjusts the estimated values as newly estimated values to reflect the estimated circuit responses at the distinct time points. The details of step 912 can be found in equations (21)- (22) and related descriptions.

Step 914 determines whether the adjusted estimated values have an acceptable accuracy to represent the circuit responses. If the determination is negative, the process is led to step 906.

If the determination is positive, the process is led to step 916. The estimated values and adjusted estimated values are in time domain.

Step 916 converts the estimated values from time domain to frequency domain. The details of step 916 can be found in Section 6, equations (31)- (34).

8. Hardware platform Figure 10 shows a block diagram illustrating an exemplary computer system 1000, which can be used as a hardware platform for executing the program that performs the processes shown in Figures 8 and 9.

As shown in Figure 10, computer system 1000 includes system bus 1001, processing unit 1002, memory device 1004, disk drive interface 1006, hard disk 1008, display interface 1010, display monitor 1012, bus interface 1014, mouse 1016, keyboard 1018, and network communication interface 1020.

The hard disk 1008 is coupled to disk drive interface 1006; monitor display 1012 is coupled to display interface 1010; and mouse 1016 and keyboard 1018 are coupled to bus interface 1014. Coupled to system bus 1001 are processing unit 1002, memory device 1004, disk drive interface 1006, display interface 1010, and network communication interface 1020.

The memory device 1004 stores data and programs. Operating together with disk drive interface 1006, hard disk 1008 also stores data and programs. However, memory device 1004 has faster access speed than hard disk 1008, while hard disk 1008 has higher capacity than memory device 1004.

Operating together with the display interface 1010, display monitor 1012 provides visual interfaces between the programs being executed and users, and displays the outputs generated by the programs.

Operating together with bus interface 1014, mouse 1016 and keyboard 1018 provide inputs to computer system 1000.

The network communication interface 1020 provides an interface between computer system 1000 and network 104 in accordance with predetermined networking protocols.

The processing unit 1002, which may include more than one processor, controls the operations of computer system 1000 by executing the programs stored in memory device 1004 and hard disk 1008. The processing unit also controls the transmissions of data and programs between memory device 1004 and hard disk 1008.

In the present invention, the program for performing the steps shown in Figures 8 and 9 can be stored in memory device 1014 or hard disk 1018. The program can be executed by processing unit 1002.

9. Summary The present invention improves the existing MFT method. The MFT method of the present invention is an efficient approach to analyzing multi-frequency nonlinear effects such as intermodulation distortion. Making the MFT method computationally efficient on problems of engineering interest required careful construction of the delay matrix, matrix-implicit Krylov subspace iterative linear solvers, and a preconditioner tailored to the MFT method and the circuits it typically analyzes. As a result, nonlinear systems comprising tens of millions of unknowns can be solved in less than an hour with computational resources commonly available to engineering designers.

One salient advantage of the MFT method in the present invention is in computing the functions 0 and the product of the Jacobian of 4) with some vectors. Both computations are essentially the solution of an initial value problem. Each application of the operator Du-jour calculation of the Newton residual, involves solving K such initial value problems, that is, integrating K sets of DAEs forward in time over one clock period. Each of the K problems, however, is essentially decoupled. Parallel implementations of the MFT will therefore enjoy very efficient processor utilization. This decoupling also assists the implementation of out-of- core solvers. In fact, it has been observed it is possible to implement the MFT algorithm as an out-of-core algorithm with over 80% average CPU utilization.

While the invention has been illustrated and described in detail in the drawing and foregoing description, it should be understood that the invention may be implemented through alternative embodiments within the spirit of the present invention. Thus, the scope of the present invention is not intended to be limited to the illustration in this specification, but is to be defined by the appended claims.