JP2002318990A - Simulation method and simulator - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a simulation method and simulator capable of enhancing the stability and precision of the tracking of a solution to reduce the number of steps of processing. SOLUTION: This simulation method comprises determining a solution by the tracking of a solution curve using BDF method. Processing is executed from process S16 in a known state up to the j-th step. The prediction and correction element in the j+1-th step are calculated to judge whether the estimated value of error is smaller than an allowable value of error or not (processes S16-S20). When the determination in process S20 is 'Yes', the following step width or solution is calculated, and the processing is advanced by one step (processes S26-S31). When the determination in process S20 is 'No', the step width is contracted, or the processing is returned back by one step (processes S21-S25). Thereafter, the processing is continued from process S16.

Description

DETAILED DESCRIPTION OF THE INVENTION

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a simulation method, and more particularly, to a DC operating point analysis (hereinafter referred to as D
The present invention relates to a simulation method and a simulator suitable for a circuit simulator for performing C operating point analysis).

[0002]

2. Description of the Related Art In LSI design, operation verification at a transistor level is performed, and a SPICE-based circuit analysis simulator is used. In a simulator that performs DC operating point analysis using the conventional Newton-Raphson method (hereinafter referred to as the NR method), if the number of transistors in the circuit exceeds several hundreds, the convergence of the NR method is poor, and a solution must be obtained. Can not. For this reason, it is necessary for the circuit designer to set an appropriate initial node potential by trial and error, which is a major problem in automating circuit design and operation verification.

Further, in the NR method, since global convergence is not guaranteed, it is necessary to give an initial value sufficiently close to a true solution in order to obtain a solution without fail. In practice,
Since it is difficult to provide this initial value in a large-scale circuit, various methods using a homotopy method that guarantees global convergence have been proposed.

Here, the principle of the homotopy method will be briefly described. In the homotopy method, the circuit equation of the circuit to be analyzed is f (x) = 0... (1), y = (x ¹ ,..., X ⁿ , t), and the solution curve of the homotopy equation is , H (x, t) = f (x) − (1-t) f (a) = 0 (2) g (x, t) = 0 (3)

Here, the variable x is x = (x ¹ ,...,
x ⁿ ), which are the voltage and current at each node of the circuit to be analyzed. t is a homotopy variable, and the value between 0 and 1 is sequentially increased, and the equation of the solution curve is sequentially changed to finally obtain the equation to be analyzed. First, an appropriate initial value is given, and by the homotopy method, Equation 3
Is calculated on a sphere having a predetermined radius and analyzed by the NR method based on the predictor to obtain a new point (corrector) on the sphere and the solution curve. Next, a spherical surface centered on the point is drawn, and the next intersection is obtained again by the NR method. The radius of each sphere is the step width. If the sphere is sufficiently small, the NR method can be expected to converge. If the NR method does not converge, the sphere is redrawn to a smaller size, that is, the step width is changed, and the NR method is applied again. By repeating this operation, the solution curve can be tracked.

A flowchart shown in FIG.
A conventional simulation method using the homotopy method according to the technique described in Japanese Patent No. 47267 will be described. First, initialization such as setting of a variable and giving an initial value in a program is performed (step S61). An initial value and a predictor for obtaining a first solution in the homotopy method are obtained (step S62), and a calculation is performed by the NR method based on the predictor to obtain a corrector (step S63).
If the NR method does not converge (step S64), an error message is immediately displayed and the process ends.

When the NR method converges in step S64 and the first corrector is obtained, the process proceeds to a step for obtaining a second solution, and from the initial value and the first corrector,
A second predictor is determined, and a modifier is determined based on the second predictor by the NR method. In this way, the process of sequentially tracking the solution curve is executed (steps S65 to S7).
4). In these processes, the step width (the radius of the sphere) in the homotopy method is selected each time according to the degree of change of the solution curve. The state of the solution curve and the step width is shown in FIG. 3, for example.

[0008] In the simulation method described in the above publication, when the solution curve is traced, an angle φ formed by the predictor and the corrector is calculated, and the angle φ is compared with a predetermined allowable value. If the angle φ is small, it is determined that the processing in that step is correct, and the process proceeds to the next processing step. On the other hand, if the angle φ between the predictor and the corrector is large, the current processing result is discarded because the step width in that step is too large. The step width is reduced again, and a predictor and a corrector are obtained.

When the homotopy variable t starts from 0 and finally becomes 1, the corrector at that time is obtained as a final solution (steps S76 to S78).

[0010]

In the conventional circuit simulation method, there is a problem that in the processing relating to the step of obtaining the first solution, a calculation method effective for the set initial value has not been established. In this case, if it is not an appropriate initial value, an error occurs immediately and the process is terminated. Therefore, it is necessary to set an appropriate initial value by trial and error.

There is also a problem that unless a predictor is properly selected, it is impossible to track a correct solution curve. FIG.
Shows a state in which the solution curve is tracked by the simulation method of FIG.

When the correct solution curve 2 and the erroneous solution curve 5 have the positional relationship shown in FIG. 1, a predictor b and a corrector c may be obtained on the erroneous solution curve 5 in some cases. In such a case, even if the tracking moves from the correct solution curve 2 to the erroneous solution curve 5, the current processing step is not discarded, the next processing step is executed, and the erroneous solution curve 5 is further moved. Will be tracked.

Further, when a solution of the simultaneous equations is obtained by the NR method, an error in the matrix calculation may increase depending on the coefficient, and the operation of the simulator operating by the NR method may become unstable.

SUMMARY OF THE INVENTION The present invention has been made to solve the above-mentioned problems of the prior art, and provides a simulation in which a solution curve can be traced with high stability and accuracy and the number of processing steps can be reduced. It is an object to provide a method and a simulator.

[0015]

In order to achieve the above object, a simulation method according to a first aspect of the present invention uses a variable step BDF method to track a solution curve of a homotopy equation, thereby achieving simultaneous equations. In the simulation method for solving the above, the step of obtaining the (j + 1) th solution on the solution curve by tracing the solution curve in the BDF method includes calculating a predictor from the homotopy equation based on the solutions up to the jth, and A first step of calculating a corrector that forms an approximate solution based on the second step; and a second step of calculating an estimated value of an error with respect to a true solution of the corrector that forms the approximate solution; Greater than the allowed value,
And a third step of determining whether the predictor is small. If it is determined that the predictor is large, a small step width is calculated as compared with the step width when the predictor is obtained. Based on the width, the first and second steps are repeated to obtain a new corrector. If it is determined that the corrector is smaller in the third step, the corrector at that time is set to the (j + 1) th approximate solution. , J + 2th solution is obtained.

According to the simulation method of the present invention, the optimum step width is obtained when the solution curve is traced by calculating the step width based on the comparison result between the estimated value of the error and the allowable value of the error. , The tracking of the solution curve in the BDF method is stable and the accuracy in the tracking is improved.

In a preferred simulation method according to the present invention, the error tolerance E is a vector formed by the j-th solution and the predictor of the (j + 1) -th solution, and the corrector of the j-th solution and the (j + 1) -th solution. Θ is the angle formed by the vector
If the constant ξ is 0 <ξ <1 and the constant ε is 0 <ε, then E
= (1 + ξcos ² θ) ε, E = (1 + ξ | cosθ |) ε,
Alternatively, it is calculated by E = (1 + ξcos θ) ε. In this case, the step width in the direction perpendicular to the solution curve is made smaller than the traveling direction in which the solution curve is tracked, so that the number of processing steps is reduced without impairing the tracking accuracy, and the solution curve tracking is efficiently performed. I can do it.

Further, a preferred simulation of the present invention
In the method, the small step width h _smallIs small
H, the constant α is 0 <α <1, and the constant η is
η <1, the estimated value of the error is T_{j + 1}Next to the BDF method
If the number is k and the function max ｛｝ is among all the elements in parentheses
H as a function to select the maximum value from_small= Max ｛αh,
ηh (E / | T_{j + 1}｜)^{1 / (1 + k)}This place calculated by｝
In this case, while keeping the step width from becoming extremely small,
Stable tracking becomes possible.

Further, in the preferred simulation method of the present invention, after it is determined in the third step of obtaining the j + 1-th solution that the size is small, the step width h _big adopted in the step of obtaining the j + 2-th solution is j + 1 The step width is h, the constant β is 0 <β <1, the constant γ is γ> 1, the estimated value of the error is T _{j + 1} , the order of the BDF method is k, and the maximum step width is h _max. And the function min ｛｝ is a function that selects the minimum value from all the elements in parentheses, h
_big = min ｛βh (E / | T _{j + 1} |) ^{1 / (1 + k)} , γh, h
It is calculated by _max ｝. In this case, the next step width is made as large as possible within a range where the estimated value of the error is smaller than the allowable value of the error, so that the number of processing steps can be reduced.

Further, in the preferred simulation method of the present invention, the lower limit h _min of the step width is h ₋₁ , the constant μ is μ> 2, and the constant ν is 1
If 1>ν> 0, it is preferable to calculate h _min = max {με, νh ₋₁ }, and estimate the error of the solution obtained one step before using αh ₋₁ as the step width. If the value is smaller than the allowable value of the error, the step width h _big of the next step is h _big = min ｛βh (E / | T _{j + 1}
|) ^{1 / (1 + k)} , {h, h _max }. As described above, a step size as large as possible is adopted as long as no error occurs in tracking, thereby preventing a reduction in tracking efficiency.

Further, in a preferred simulation method according to the present invention, in the first step, a tangent of a solution curve at a point of t = 0 and an intersection of a plane of t = 0 and the solution curve are set, where t is a homotopy variable. The intersection at which t is positive among the intersections with a sphere having a predetermined radius centered on is set as a predictor. In this case, the predictor for obtaining the first solution is optimally set, and the corrector is easily obtained based on this.

A simulation method according to a second aspect of the present invention is a method for solving a homotopy equation using the NR method. For each of the simultaneous equations constituting the homotopy equation, a coefficient of a variable of each of the simultaneous equations is used. The maximum coefficient is obtained from among the above, and the coefficients of all the variables in each simultaneous equation are divided by the maximum coefficient. In this case, since the absolute difference between the coefficients in each equation is small, N
It can be expected that the error when solving the simultaneous equations by the R method becomes small. By reducing the calculation error, stable tracking of the solution curve becomes possible.

Further, a simulation method according to a third aspect of the present invention includes a step of creating a system of simultaneous equations for a circuit based on netlist information representing the characteristics of each electronic component constituting the electronic circuit and the connection relation thereof; A solution step of solving the system of simultaneous equations, and a simulation method for analyzing the operation of the electronic circuit comprising the step of solving the system of simultaneous equations, wherein the homotopy equation is constructed from the system of simultaneous equations using a homotopy method, and a variable step is performed. A step of finding a solution of the system of simultaneous equations by tracing a solution curve of the homotopy equation using a BDF method, and a step of finding a j + 1-th solution on the solution curve by tracing the solution curve in the BDF method, Based on the solutions up to j-th, a predictor is calculated from the homotopy equation, and a corrector that forms an approximate solution based on the predictor is calculated. A second step of calculating an error estimate for the true solution of the modifier forming the approximate solution, and a second step of calculating whether the error estimate is larger or smaller than an allowable value of the error. And if it is determined that the predictor is large, calculate a step width smaller than the step width when the predictor was obtained, based on the small step width The first and second steps are repeated to obtain a new corrector. If it is determined that the corrector is smaller in the third step, the corrector at that time is determined as the (j + 1) th approximate solution. Characterized by shifting to the step of finding the solution of In this case, the analysis of the operation of the electronic circuit using the homotopy method can find the solution faster and more reliably than in the past.

In the simulation method according to the first aspect of the present invention, instead of the above configuration, the allowable value of the error is set to a predetermined range. Alternatively, it may be determined whether the value is smaller than the minimum value of the range. In this case, if the estimated value of the error is within a predetermined range, the previous step width is automatically adopted.

[0025]

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS Hereinafter, a simulation method according to the present invention will be described based on an embodiment of the present invention with reference to the drawings. FIG. 1 shows the operation of a circuit simulator in which a simulation method according to an embodiment of the present invention is incorporated. The circuit simulator is loaded with a SPICE-based circuit analysis program, performs a DC analysis first, and is programmed to perform an AC analysis and a transient analysis based on the result of the DC analysis.

First, circuit data of a circuit to be analyzed is read from a circuit data file (step S1), and a circuit equation for DC operating point analysis is created based on the circuit data (step S2).

The circuit data includes, for example, net list information indicating the characteristics of each electronic component constituting the electronic circuit and the connection relationship thereof. The circuit equation is a simultaneous equation using the contact potential and the contact current of the analysis target circuit as variables, and has the same number of equations as the sum of the total number of node potentials of the analysis target circuit and the number of voltage sources used.

First, a first predictor is obtained based on an appropriate initial value and a gradient of a solution curve, and based on this predictor, a solution of a circuit equation is obtained using a normal NR method ( Step S3). When the NR method converges and a solution is obtained (step S4), the process proceeds to step S6 to execute the next process.

If no solution is found, a solution of the circuit equation is found by the homotopy method (step S5). In solving the circuit equation by the homotopy method, the simulation method according to the embodiment of the present invention is employed. The result obtained by solving the circuit equation is written in a result file, and the process ends (step S6).

FIG. 2 is a flowchart showing the details of step S5 in FIG. The simulation method of the present embodiment is executed based on the homotopy method.

First, an initial step width and parameters are set and initialized (step S11). The predictor in the first step is calculated by the processing described below (step S1).
2).

FIG. 7 schematically shows the first step of the homotopy method. Assuming that the homotopy variable is t and the intersection between the correct solution curve and the t = 0 plane is a, an intersection b between a tangent of the solution curve drawn to the intersection a and a sphere having a radius h ₀ centered on the intersection a is obtained. This is used as the first step predictor.

Specifically, it is calculated as follows. Jf (a) · x = −f (a) (4) where Jf is the Jacobian of f. Assuming that the initial step width is h, the predictor yp to be obtained is expressed by the following equation using the solution x of the equation (4). Where y ^{n + 1} = t = 1
It is. When i = 1, ..., n

(Equation 1) When i = n + 1

(Equation 2)

Using the predictor obtained in step S12 as an initial value, a homotopy is calculated using a solution curve tracing algorithm (hereinafter referred to as a BDF method) in which a predictor modifier method is applied to a first-order backward difference formula by the NR method. The equation is solved to determine a first-step corrector (step S13).

It is checked whether the NR method converges (step S14).
If the predictor in the first step can be calculated, the process proceeds to step S16. If the predictor cannot be calculated, the step width h is reduced (step S15), and the process is executed from step S12.

When the predictor in the first step can be calculated, 2
A predictor after the step is obtained (step S16). In the second step, a primary predictor is calculated, and in the third and subsequent steps, a secondary predictor is calculated.

[0037] j-th step y, h and y _j, is expressed as h _j, j + 1 th predictor yp _{j + 1} are shown as follows. The primary predictor in the second step is yp _{j + 1} = y _j + ω _{j + 1} (y _j −y _j−1 ) (7) The secondary predictors in the third and subsequent steps are: yp _{j + 1} = y _j + ω _{j + 1} (y _j −y _j−1 ) + ω _{j + 1} ω _j (1 + ω _{j + 1} ) / (1 + ω _j ) × ｛(y _j −y _j−1 ) −ω _j (y _j-1 −y _j-2 )｝ (8) where ω _j = h _j / h _j-1 .

Using the predictor obtained in step S16 as an initial value, the homotopy equation is solved by the NR method using the BDF method, and a corrector for the second and subsequent steps is obtained (step S1).
7). In the second step, the first-order BDF method is used to calculate
After the step, calculation is performed by the second-order BDF method.

The modifier yc _{j + 1} at the j + 1th step is j +
It can be obtained by solving the following equation with the predictor yp _{j + 1} of the first step as an initial value. h (yc _{j + 1} ) = 0 (9)

(Equation 3) The second step of the first-order BDF _{^{method, dy i j + 1 = yc}} i j + 1 -y i j-1 ···· (11) 3 steps subsequent second order BDF method, dy ⁱ _{j +1 = (1 + 2ω j +} 1) / (1 + ω j + 1) × yc i j + 1 - (1 + ω j + 1) y i j + ω 2 j + 1 / (1 + ω j + 1) × y i j-1 .... (12)

It is checked whether the NR method converges (step S18).
If the second and subsequent modifiers cannot be calculated, step S2
When the processing is executed from step 1 and calculation is possible, an estimated value (hereinafter, referred to as an error estimated value) _{Tj + 1 of} the local discretization error and an allowable value E of the error are calculated based on the predictor and the corrector. (Step S19).

If j is the number of steps, the estimated value of the error T
_{j + 1} is for the first step (j = 1), T _{j + 1} ｜| y−yc | ≒ | yp−yc | (13) For the second and subsequent steps (j> 1) , T _{j + 1} } | y−yc |｜ ω _{j + 1} (1 + ω _{j + 1} ) / ｛ω _{j + 1} (1 + ω _{j + 1} ) + (1 + 2ω _{j + 1} ) (ω _j−1 + ω _j ⁻¹ +1)｝ × | yp−yc | (14) Expression 13 or Expression 14 shows the estimated value of the error between the correct solution at the j + 1th step and the corrector.

The allowable value E of the error is defined as θ, where an angle formed between a vector composed of the predictors of the solutions of the j-th and j + 1-th steps and a vector composed of the correctors of the solutions of the j-th and j + 1-th steps is θ, and the constant ξ is 0. <Ξ <1, and the constant ε is 0 <ε.
The allowable value E of the error is calculated by one of the following equations. E = (1 + ξcos ² |) εε (15) E = (1 + ξ | cos θ |) ・ (16) E = (1 + ξcos θ)） (17)

[0043] investigated whether the estimated value T _{j + 1} Do error tolerance E is smaller than the error (step S20), the estimated value of the error T _{j + 1}
Is larger, the lower limit h of the step width is as shown below.
_min is calculated, and it is determined whether the step width h _{j + 1} of the ( _{j + 1) th} step is larger than the lower limit h _{min of the} step width (step S2).
1).

The lower limit h _min of the step width is defined as h, the step width of the current step is h, the step width one step before the current step is h ₋₁ , the constant μ is μ> 2, and the constant ν is 1> ν
If> 0, h _min = max {με, νh ₋₁ } (18)

The current step width h is the lower limit h of the step width.
_{If it is} larger than _min , the second and subsequent modifiers cannot be calculated, and if the determination in step S18 is "No", the step width is reduced to a new step width αh instead of the current step width h.

If the estimated value of the error is larger than the allowable value of the error and the determination in step S20 is "No", the step width is reduced as described below (step S22). Thereafter, the processing is performed from step S16, and the calculation of the (j + 1) th step is performed. However, the constant α is α <1.

The small step width h _small is a step width before being reduced, the constant α is 0 <α <1, the constant η is η <1, the estimated value of the error is T _{j + 1} , the BDF method If the order of k is k and the function max ｛｝ selects the maximum value from all elements, Is calculated by

The small step width h _small is h _small <h as shown in Equation 19. The meaning of each term used in Equation 19 is shown below. Is in the range of 0 <α <1, and 0.6 to 0.75 is used for the actual value.
When the constant α is smaller than 1, the first element αh is made smaller than h, and the small step width h _small is changed to the previous step width h.
Smaller than

The estimated error value T _{j + 1} is predicted to be smaller than the allowable error value E. Equation 19 shows an approximation of the largest step width among them, and guarantees that it is smaller than h.

The constant η is in the range of 0 <η <1 for safety, and 0.7 to 0.95 is used as the actual value. The function max ｛｝ suppresses the step width from becoming extremely small and the change in the step width from becoming large, so that the estimated value of the error is accurate and the solution curve can be stably tracked.

The current step width h is the lower limit h of the step width.
_{If it is} smaller than _min , the current processing is checked (step S23),
If the process returns one step, the calculation cannot be continued, so an error message is displayed and the process ends with an error. If the process is not one step back, the current process is returned one step from the jth step to the j-1th step (step S24).

Next, the step width is reduced by replacing the current step width h with the new step width αh (step S2).
5) After that, the processing is executed from step S16 to calculate the j-th step.

In the calculation of the (j + 1) th step, the following equation (20) is used.
, The next step width is calculated (step S26).

The large step width h _big calculated by the (j + 1) th step is such that the large step width is h _big , the step width before the increase is h, the constant β is 0 <β <1, and the constant γ
, Γ> 1, the maximum step width is h _max , the estimated value of the error is T _{j + 1} , the order of the BDF method is k, and the function min ｛｝
Selects the minimum value among all the elements, h _big = min ｛βh (E / | T _{j + 1} |) ^{1 / (1 + k)} , γh, h _max ｝ (20) ).

If it is determined that the estimated value of the error is smaller than the allowable value of the error, an approximate value of the solution by the modifier is set as the solution of the (j + 1) th step, and the large step width h _big is set to j +
Set as the step width of the second step.

The large step width h _big by the calculation of the (j + 1) th step can be obtained by using the expression 20 in a range where the estimated value of the error is smaller than the allowable value of the error if the estimated value of the error is smaller than the allowable value of the error. , And make the next step width as large as possible. The meaning of each term used in Equation 20 is shown below.

The estimated error value T _{j + 1} is predicted to be smaller than the allowable error value E. Equation 20 shows an approximation of the largest step width among them and guarantees that it is larger than h.

The constant β of the first element is set to 0 for safety.
<Β <1, and the actual value is 0.9 to 0.95.
Is used.

The constant γ of the second element is in the range of 1 <γ, and a value of 1.8 to 2.0 is used as an actual value. The second element γh keeps the ratio between the step width h before the enlargement and the large step width h _big within a certain range, thereby preventing the error estimation value from being inaccurate and the solution curve from becoming unstable. I do.

The third element h _max is a value that determines the absolute upper limit of the step width so that the step width does not become extremely large, and an actual value of about 10.0 is used.

The function min ｛｝ selects the minimum value from the first to third elements, thereby suppressing the step width from becoming extremely large and smoothing the change in the step width.

[0062] large step width by j-th step of computing h _big is the lower limit h _min is smaller than the step width h is the step width of the current _{step, h big = min {βh (} E / | T j + 1 |) 1 ^{/ (1 + k), ξh} , is calculated by h _max} ···· (21). ξ is a constant satisfying ξ <γ.

It is determined whether the homotopy variable has become 1.
Specifically, with respect to homotopy variable ^{_{t j + 1 = y n +}} 1 j + 1, t j <1 <t j + 1 or _{t j>1> t j +} 1 is either determined (step S27 ), If the homotopy variable is not 1,
The process is executed from step S30.

If the homotopy variable becomes 1, x _j and x
_{j + 1} is linearly interpolated to obtain an approximate value of the solution x at the time t = 1 (step S28), the equation 1 is solved by the NR method with the approximate value of the solution x as an initial value, and the result is output (step S29). ).

It is checked whether the number of iterations of the homotopy method indicating the number of loop iterations in steps S16 to S31 is larger than a preset value _Nmax indicating the upper limit of the loop (step S3).
0). The set value _Nmax is set to 6 to 17.

If the homotopy method repetition number is large, the process is terminated. If the homotopy method repetition number is small, y _{j + 1} = yc
_The current processing is advanced by one step as _{j + 1} (step S3
1). Thereafter, the processing is executed from step S16 to j + 2
The calculation of the step is performed.

In the above algorithm, the second-order BDF method has been described, but it is easy to use a third-order or higher-order BDF method.

FIG. 3 shows control for obtaining an optimum step width by the simulation method of FIG. Correct solution curve 2
The solid line indicates a tracked portion and the broken line indicates an untracked portion.

The simulation method uses the correct solution curve 2
When the curvature is large, the step width is reduced, and when the curvature is small, the step width is increased. In this case, the movement of the error to the solution curve is avoided, the error is suppressed, and the number of steps is reduced, so that the solution curve is efficiently tracked.

FIG. 4 shows control for keeping the corrector within an allowable error range by the simulation method of FIG. Since the simulation method evaluates the estimated value of the error and tracks the solution curve, the approximate curve 3 represented by the corrector is out of the tolerance 4 of the error indicating a certain range from the correct solution curve 2. Do not move to. In this case, the solution tracking accuracy is high.

FIG. 5 shows how the simulation method of FIG. 2 effectively tracks a solution. In the correct solution curve 2, a solid line indicates a tracked portion and a broken line indicates an untracked portion.

The simulation method is based on Equations 13 to 15.
The error tolerance is calculated based on the following formula, so that the error tolerance increases in the same traveling direction as the direction in which the solution curve is tracked, and the error tolerance in the vertical direction perpendicular to the solution curve tracking direction. The value decreases. In this case, since the step width in the vertical direction is made smaller than the traveling direction of tracking the solution curve, the number of processing steps is reduced, and the solution curve can be efficiently tracked.

FIG. 6 shows how the simulation method of FIG. 2 tracks again from the erroneous solution curve to the correct solution curve. The figure shows a state where steps S20, S21 and S22 are executed. The correct solution curve 2 and the erroneous solution curve 5 approach and have a predetermined positional relationship. 1A shows a sphere having a discarded step width, and 1B shows a sphere having a new step width.

As shown in step S19, the estimated value of the error is proportional to the difference bc between the predictor c and the corrector b. Process S
At 20, the difference b ₁ c ₁ between the spheres 1A of the step width is large, and the estimated value of the error is larger than the allowable value of the error. Step S22
The ball 1A having the step width is discarded, and the ball 1 having the step width is discarded.
B is newly adopted. The sphere 1B having the step width has a difference b ₂
c ₂ is small, the estimated value of the error is smaller than the tolerance of errors.

The simulation method uses an error solution curve 5
, The estimated value of the error becomes larger than the allowable value of the error, and a small step width is calculated based on Equation 19, and the correct solution curve 2 is tracked again.

In the simulation method, when the estimated value of the error is smaller than the allowable value of the error, the process is advanced by one step, and the step width is calculated based on Expression 20. Equation 21 is the first
To have a third element to ensure efficiency and stability.

The first element is an estimated value that gives the maximum step width when the estimated value of the error at the next step is within the allowable value of the error. The second element is a request of the principle of the present invention (h _j ≒
h _{j + 1} ) The estimated value of the error at the next step is within the allowable value of the error, and is the estimated value that gives the maximum step width. The third factor is a value that prevents the step width from becoming excessively large and ensures stability.

The lower limit of the step width is calculated based on Expression 18, and by setting it sufficiently larger than με and the allowable amount of error ε, the sphere used in the BDF method intersects with the solution curve. Guaranteed.

In step S26, once the process returns one step, the next step width is calculated based on expression 21 instead of expression 20. When the process returns one step at a time, the change of the solution curve is expected to be large, and a small step width is desired.

The predictor in the first step is calculated according to step S12. The predictor and the corrector are expected to be sufficiently close, since the solution curve is tracked in the reverse direction and the NR method is prevented from converging.

FIG. 7 shows that the simulation method of FIG.
4 shows how a predictor at the step is calculated. Intersection a is
This is the point where the straight line at x = 0 and the correct solution curve 2 intersect.
The sphere 1 having a step width is a sphere having the intersection a as a center and a radius having a step width according to the equation (10). The corrector c of the first step is obtained from the predictor by the NR method, and the correct solution curve 2
And a sphere 1 having a step width.

In the conventional simulation method, the predictor e is set to an appropriate value. The predictor e is, for example, x = x
0, t = step width, far away from the desired modifier c, the NR method using the predictor e as an initial value does not converge, and the desired modifier c may not be obtained. Also, even if the NR method converges, an undesirable modifier d with t <0 may be required.

In the simulation method shown in FIG.
The predictor b is obtained in accordance with Predictor b is expected to be sufficiently close to the desired modifier c. Predictor b is N
If used as the initial value of the R method, the NR method converges.
A desirable modifier c is sought.

In step S12, a maximum value is obtained from the coefficients of all the variables for the simultaneous equations of the circuit for each equation, the coefficients of all the variables are divided by the maximum value, and the homotopy equation is calculated by the NR method. solve. In this case, since the difference between the absolute values of the coefficients of the equations becomes smaller, the calculation error generated when solving the simultaneous equations by the NR method becomes smaller, so that the solution curve can be tracked stably.

FIG. 8 shows how the simulation method of FIG. 2 sets the optimum step width. Although the step width is reduced to some extent, the step width cannot be reduced without limit.

In the traced approximate analytical curve 6, when the step width is extremely reduced with respect to the allowable value E of the error, the intersection between the sphere 1 having the step width according to the equation (10) and the correct solution curve 2 disappears. So I can't ask for a modifier. In addition, the influence of the calculation error becomes large, and the correct corrector may not be obtained.

Since the lower limit h _min of the step width is calculated using Expression 18, the step width is suppressed so as not to be extremely small. The lower limit of the step width h _min = max in Expression 18
For {με, νh _-1 }, the following variables are used.

The constant μ is set to μ> 2. The first element [mu] [epsilon] prevents the step width from becoming extremely small and smoothes the change in the step width so as to find the corrector. The constant ν is set to 0 <ν <1, and the actual value is
0.1-0.5 is used. Lower limit of step width h _min
, The larger of the first element and the second element is selected.

FIGS. 9A and 9B show a verification circuit in which the simulation method of FIG. 2 is executed and a simulation result. As shown in FIG.
It is a closed circuit in which a voltage source, a resistor, and two tunnel diodes are connected in series. The circuit simulator obtains the potential at the node A where the two tunnel diodes are connected.

As shown in FIG. 9B, the simulation result is obtained by solving a circuit equation for DC operating point analysis,
This is a drawing of a solution curve for. The straight line of the homotopy variable t = 1 and the solution curve intersect at five points, indicating that the solution curve was correctly tracked.

FIG. 10 shows simulation contents calculated by the circuit simulator of FIG. The circuit simulator executes a simulation for four types of circuits A to D: an I / O buffer, a clock net, a PLL, and an A / D. In the circuits A to C, the solution of the circuit equation was obtained by the normal NR method, and in the circuit D, the solution of the circuit equation was not obtained by the normal NR method, but was obtained by the homotopy method.

FIG. 11 shows a simulation result of the circuit D of FIG. The simulation result is a part of the solution curve of the circuit D traced by the homotopy method. The circuit D is an 8-bit 8 MHz A / D converter including 3111 total nodes and 6839 MOSs. A commercially available circuit analysis simulator cannot analyze the circuit D because the NR method of DC operating point analysis for the circuit D does not converge.

As shown in the figure, the two solution curves respectively
It is a curve for calculating the node potentials of the node names RIN0 and TBO6, and the solution is traced from the homotopy variable t = 0 to t = 1. The node potential to be determined is the intersection of the solution curve and the straight line at t = 1. The same applies to all other node potentials not shown.

According to the above-described embodiment, a small step width is calculated based on the comparison result between the estimated value of the error and the allowable value of the error to control the optimum step width. Curve tracking becomes highly stable and accurate.

Further, since the step width in the vertical direction is made smaller than the traveling direction in which the solution curve is tracked, the number of processing steps is reduced, and the solution curve can be tracked efficiently.

Although the present invention has been described based on the preferred embodiment, the simulation method and the simulator of the present invention are not limited to the configuration of the above-described embodiment, and are not limited to the configuration of the above-described embodiment. Simulation methods and simulators with various modifications and changes from the configuration are also included in the scope of the present invention.

For example, the above simulation method is not limited to circuit design, but can be applied to all simultaneous equations for which a solution is obtained using a simulator.

[0098]

As described above, in the simulation method of the present invention, the tracking of the solution curve in the BDF method becomes highly stable and accurate, and the number of processing steps is reduced.
Since solution curve tracking can be performed efficiently, DC
Effective for analysis.

According to the present invention, in a simulation method for analyzing the operation of an electronic circuit, a solution can be obtained at higher speed and more reliably than a conventional simulation method using a homotopy method.

[Brief description of the drawings]

FIG. 1 shows an operation of a circuit simulator in which a simulation method according to an embodiment of the present invention is incorporated.

FIG. 2 is a flowchart showing details of step S5 in FIG.

FIG. 3 shows control for setting an optimum step width by the simulation method of FIG. 2;

FIG. 4 shows control for suppressing a corrector within an allowable range of an error by the simulation method of FIG. 2;

FIG. 5 shows how the simulation method of FIG. 2 effectively tracks a solution.

FIG. 6 illustrates how the simulation method of FIG. 2 tracks again from an erroneous solution curve to a correct solution curve.

FIG. 7 shows how the simulation method of FIG. 2 calculates a predictor in the first step.

FIG. 8 shows how the simulation method of FIG. 2 sets an optimum step width.

FIGS. 9A and 9B show a verification circuit in which the simulation method of FIG. 2 is executed and simulation results.

FIG. 10 shows simulation contents calculated by the circuit simulator of FIG. 1;

11 shows a simulation result of the circuit D of FIG.

FIG. 12 shows a simulation method according to the technique described in Japanese Patent Application Laid-Open No. 8-147267.

FIG. 13 shows how a solution curve is tracked by the simulation method of FIG.

[Explanation of symbols]

 1 Step width sphere 2 Correct solution curve 3 Approximate curve represented by modifier 4 Error tolerance 5 Error solution curve 6 Traced analysis curve

Claims (12)

[Claims] 1. A simulation method for solving simultaneous equations by tracking a solution curve of a homotopy equation using a variable step BDF method, wherein j on a solution curve obtained by tracking the solution curve in the BDF method is
A step of calculating a predictor from the homotopy equation based on the solutions up to the j-th solution, and calculating a corrector that forms an approximate solution based on the predictor; A second step of calculating an estimated value of an error with respect to a true solution of a corrector forming a solution; and a third step of determining whether the estimated value of the error is larger or smaller than an allowable value of the error. When it is determined that the predictor is large in the third step, a step width smaller than the step width when the predictor is obtained is calculated, and the first and second steps are performed based on the small step width. Is repeated to obtain a new corrector. If it is determined that the corrector is small in the third step, the corrector at that time is set as the (j + 1) th approximate solution, and then the process proceeds to the step of obtaining the (j + 2) th solution. Characteristic simulation method. 2. The error allowable value E is obtained by calculating a j-th solution and j +
The angle formed by the vector formed by the predictor of the first solution and the vector formed by the j-th solution and the modifier of the (j + 1) -th solution is θ, the constant ξ is 0 <ξ <1, and the constant ε is 0 <When ε, E = (1 + ξcos 2 θ) ε, E = (1 + ξ | cosθ |)
The simulation method according to claim 1, wherein the calculation is performed using ε or E = (1 + ξcos θ) ε. 3. The small step width h _small is a step width before being reduced, h is a constant α is 0 <α <1, η is a constant η <1, the estimated value of the error is T _{j + 1} , Assuming that the order of the BDF method is k and the function max ｛｝ is a function that selects the maximum value from all the elements in parentheses, h _small = max ｛αh, ηh (E / | T _{j + 1} |) ^{1 / (1 + k)}} in the calculation, simulation method according to claim 1 or 2. 4. The step width h _big adopted in the step of _obtaining the j + _2th solution after it is determined in the third step of obtaining the j + 1th solution that the size is small is h, where the j + 1th step width is h and the constant is β is 0 <β <1, constant γ is γ>
1, T _{j + 1} estimates of the error, the order of k of the BDF method, the step width of the maximum limit and h _max, the function min {} selects the minimum value from among all the elements in parentheses As a function, h _big = minＥβh (E / | T _{j + 1} |) ^{1 / (1 + k)} , γh, h
The simulation method according to claim 1, wherein the calculation is performed by _max ｝. 5. The lower limit h _min of the step width is as follows: h _min = max where h _{−1 is} the step width at the time of finding the immediately preceding solution, μ> 2 is the constant μ, and 1>ν> 0 is the constant ν. The simulation method according to claim 3, wherein the calculation is performed by {με, νh _-1 }. 6. When the small step width h _small calculated for the current step is smaller than the lower limit h _{min of the} step width, the process for finding the solution is returned one step before, and αh ₋₁
The simulation method according to claim 5, wherein a solution is obtained using the step width as a step width. 7. If the estimated value of the error is smaller than the allowable value of the error for the solution obtained one step before using the αh ₋₁ as the step width, the step width h of the current step, which is the next step, is used. _big is given by ξ <γ, h _big = min ｛βh (E / | T _{j + 1} |) ^{1 / (1 + k)} , ξh, h
The simulation method according to claim 6, wherein the calculation is performed by _max ｝. 8. In a first step, a tangent of a solution curve at a point of t = 0 and t = 0
The intersection of a sphere having a predetermined radius centered on the intersection of the plane and the solution curve is defined as a predictor where t is positive.
The simulation method according to any one of the above. 9. A method of solving the homotopy equation using the NR method, wherein for each of the simultaneous equations constituting the homotopy equation, the largest coefficient among the coefficients of the variables of each of the simultaneous equations is determined. The simulation method according to any one of claims 1 to 8, wherein coefficients of all variables in each simultaneous equation are divided by the following coefficient. 10. A simulator using the simulation method according to claim 1 to solve a simultaneous equation of an electric circuit. 11. A method for solving a homotopy equation using the NR method, wherein for each of the simultaneous equations constituting the homotopy equation, the largest coefficient among the variables of each simultaneous equation is determined. Dividing the coefficients of all the variables in each of the simultaneous equations by the following method. 12. The method according to claim 1, further comprising: a step of creating a system of simultaneous equations for the circuit based on netlist information representing the characteristics of each electronic component constituting the electronic circuit and the connection relationship thereof; and a solution step of solving the system of simultaneous equations. In a simulation method for analyzing the operation of an electronic circuit, the step of solving the simultaneous equations constructs a homotopy equation from the system of simultaneous equations using a homotopy method, and tracks a solution curve of the homotopy equation using a variable step BDF method. Calculating the solution of the system of simultaneous equations by performing j on the solution curve by tracking the solution curve in the BDF method.
A step of calculating a predictor from the homotopy equation based on the solutions up to the j-th solution, and calculating a corrector that forms an approximate solution based on the predictor; A second step of calculating an estimated value of an error with respect to a true solution of a corrector forming a solution; and a third step of determining whether the estimated value of the error is larger or smaller than an allowable value of the error. If it is determined that the predictor is large in the third step, a step width smaller than the step width when the predictor is obtained is calculated, and the first and second steps are performed based on the small step width. Is repeated to obtain a new corrector. If it is determined that the corrector is small in the third step, the corrector at that time is set as the (j + 1) th approximate solution, and then the process proceeds to the step of obtaining the (j + 2) th solution. Characteristic simulation method.