2005 RISP International Workshop on Nonlinear
Circuit and Signal Processing (NCSP'05)
Hawaii, USA, Mar. 4-6, 2005.

Tracing Bifurcation Branches of High-Dimensional Nonlinear Algebraic Equations
Using SPICE
Seizo Hagino† , Yoshifumi Nishio† and Akio Ushida†††

Tokushima University
2-1 Minami-Josanjima, Tokushima, Japan
Phone:+81-88-656-7470, FAX:+81-88-656-7471
Email: hagino@ee.tokushima-u.ac.jp††
Tokushima Bunri University

abstract: It is very important to calculate the multiple
solutions of nonlinear equations, because there are manykind of practical problems in nonlinear electronic circuits
and the nonlinear equations in many scientiﬁc ﬁelds.
Therefore, there have been already proposed many algorithms for calculating multiple and/or all the solutions.
One of the powerful methods is a Newton homotopy
technique which can ﬁnd out the multiple and/or all
the solutions lying on the homotopy paths in the global
space. Unfortunately, it cannot trace all the solutions
when the homotopy path has bifurcation points on it,
because the path happens to have separated branches
at the point. In this paper, we propose an eﬃcient algorithm to ﬁnd the bifucation points at which a rank of theJacobian matrix is reduced more than one. After ﬁnding the points, we ﬁne the directions of branches, andtrace them with the arc-length method. Our algorithm
is based on the use of Spice and needs not to derive theJacobian matrix, so that it can be easily applied to manykinds of nonlinear equations.
1. Introduction
Numerical analysis of nonlinear algebraic equations is
one of the most fundamental and important problems inscientiﬁc and engineering ﬁelds. The problems in circuitdesigns are corresponding to ﬁnd the multiple DC operating points, the oscillator designs, and so on. Althoughin Spice simulator, Newton-Raphson method is mainly
used for this purpose, it can be only applied when theapproximate solutions are given. Thus, it cannot any
more use to ﬁnd the multiple and/or all the solutions.
Therefore, there have been published many papers [1]-[5]for ﬁnding the multiple and/or all the solutions of thefollowing nonlinear equation:
f1 (v1 , v2 , · · · , vn ) = 0 
f2 (v1 , v2 , · · · , vn ) = 0·········
fn (v1 , v2 , · · · , vn ) = 0

Some of them [2]-[3] can eﬃciently ﬁnd all the solutions of piecewise-linear resistive circuits, where the nonlinear elements are modeled by the piecewise-linear segments. Although the interval method [4] was proposed asa method for ﬁnding all the solutions, it is a rather timeconsuming for large scale systems. On the other hand,
the homotopy methods [5]-[10] have global convergenceproperty, so that they can be applied to wide classesof nonlinear equations. Especially, they are applied to
calculate the multiple solutions of nonlinear equations.
The Newton-homotopy method [7] coupled with the solution curve tracing algorithm [10] is possible to trace
all the solutions existing on the homotopy paths which
can start from any initial value v0 . However, it is notknown whether it has found all the solutions or not, because the homotopy paths may consist of the multiple
independent branches and the closed loops in the space.
For an example, many branches cross together at the
points called pitchfork bifurcation points. Therefore, it
is important to calculate the bifurcation points to ﬁnd
the multiple solutions, and to know the behaviors of thenonlinear eq.(1).
In this paper, we show the curve tracing algorithm
using SPICE in section 2. The Spice-oriented algorithms
to ﬁnd the pitchfork bifurcation points, and tracing thebranches starting at the points are shown in section 3
and 4.

2. Curve tracing algorithm
Nowadays, Spice is widely used for circuit simulationsof integrated circuits. Since nonlinear devices in ICs
are modeled by many kinds of the functions, Spice has
all kinds of controlled sources, and nonlinear functions
such as exponential functions, multiplications, additions
and subtractions and so on. Furthermore, we can easily
construct abitrary functions using ABMs(analog behavior models). An example of multiple-inputs and single(1)output function is shown in Fig.1, where the function
f (v1 , v2 , . . . , vn ) in the black box is realized by writing- 451 -

the function with the input variables {v1 , v2 , . . . , vn }. On
the other hand, a set of the nonlinear equations given byeq.(1) is realized by the equivalent circuit as shown in
Fig.2.

arc-length method eq.(4) is shown by Fig.3(b). That iscalled STC(solution curve tracing circuit)[8]. In many
practical problems [7], we can get all the DC solutions
from a suitable initial guess I0 in eq.(3).
ρoutput

input

v1
v2
...

Analog Behavior Model

vn

ρ-1
VCCS

-1
STC

+
f(v1,v2,...,vn)-.

I 0 i , i=1,2,...nnFigure 1: Analog Behavior Model.

Nonlinear(ρ−1)I0i

resistive circuit

i=1,2,...n.

ρ (=v n+1 )Initial guessvi i=1,2,...n.2 +v= Σ vi
i=1

.ρ2

v

-

ρ

R
1[A]
INTEG

vccs

(a)

vccs

(b)

Figure 3: (a)Equivalence circuit of Newton homotopy method.

+

+

+

v1 1[Ω]

v2 1[Ω]--(b)Solution curve tracing circuit.

... vn 1[Ω]-f1(v1,v2,...vn)+v1 f2(v1,v2,...vn)+v2 fn(v1,v2,...vn)+vn
Figure 2: Equivalence circuit of the nonlinear algebraic equation.

3. Calculations of bifurcation points
Now, consider an algorithm for calculating the bifurcation points of the nonlinear equations. In this section,
we show an eﬃcient Spice-oriented algorithm. For simplicity, we set vn+1 = ρ and rewrite the relation eq.(3)as follows:

Namely, if we assume the voltage at i-th resistor (1Ω)
by vi , we can realize the current source as follows:
Ii = fi (v1 , v2 , . . . , vn ) − vi(i = 1, 2, . . . , n)(2)Thus, we have a relation fi (v1 , v2 , . . . , vn ) = 0 with
the Kirchhoﬀ’s voltage law. Let us apply the Newton
homotopy method to solve the eq.(1). Then, we formulate
the following modiﬁed relation:F(v, ρ) = f (v) + (ρ − 1)I0 = 0(3)

where f (v) : Rn → Rn given by the relation eq.(1),and I0 = f (v0 ) is the initial guess. Thus, the relation
eq.(2) is consists “n” equations having “n + 1” variables,
so that the solution satisfying the relation eq.(3) is givenby the solution curves in “(n + 1)” dimensional space.
In order to trace the solution curve with the arc-length
method, we apply the following relation [10]:
2
K 

dvi
i=1

ds

=1

(4)

F(v) = 0

F : Rn+1 → Rn(5)

It is proved that our algorithm [10] can trace the solution curves only if the rank of the Jacobian matrix ofeq.(5) is n on the whole curve. Note that the rank at
the limit points (turning points) is equal to “n”, and therank is less than or equal to n − 1 at the pitchfork bifurcation points. Therefore, the limit points can be still
traced by our arc-length method. However, it is generally not so easy to calculate the Jacobian matrix and therank of eq.(5). Thus, let us to develop an algorithm tocalculate the bifurcation points without estimating theJacobian matrix.
At the pitchfork point of eq.(5), we consider a set ofthe following modiﬁed equations1 :
F1 (v) = F(v) = 0
F2 (v, ∆v) = F(v + ∆v) − F(v) = 0(6)which have 2n + 2 variables and 2n equations if we
assume ∆v as the additional variables. Since the rankof the Jacobian matrix of F(v) is equal or less than n − 1[11], if we set the variables ∆v1 and ∆v2 with small ﬁxed
constants in F2 (v, ∆v)), then, eq.(6) becomes the nonlinear equation with “2n” variables and “2n” equations.
Hence, we can also apply the Newton homotopy method
as follows:

Observe that the solutions satisfying eq.(1) are obtained at ρ = 1 on the solution curve of eq.(3). Now,
we show the equivalent circuit satisfying the Newton homotopy method eq.(3) and eq.(4). The block diagram isshown by Fig.3(a). Some of the “K” node voltages in
1 At the pitchfork points, the second order term of the Taylor
the eq.(1) should be chosen as the variables of the arc- expansions are assumed negligible values for the small variations
length method eq.(4) [7]. Thus, the circuit realizing the ∆v- 452 -


F1 (v) + (λ − 1)I01 = 0
F2 (v, ∆v) + (λ − 1)I02 = 0
K

dvi

K

d∆vi

2

2

2

method [12] using the Jacobian matrix of eq.(5), especially for the large scale systems.
(7)

4.
Tracing solution curves from bifurcation
pointsi=1
i=3Now, consider to calculate the branches starting fromTracing the homotopy path eq.(7), we can ﬁnd the bi- the pitchfork bifurcation points in eq.(5). We set a small
furcation points at λ = 1. Note that two of the variables sphere centered at the bifurcation points v̂, and get the∆v1 and ∆v2 should set to ﬁxed suﬃcient small values. intersections between the branches of eq.(5) and sphere.
Thus, we have
Example 3.1. Now, let us calculate bifurcationKpoints of the following equation:
F(v) = 0 ,(vi − v̂i )2 = r 2
(11)

i=1
F1 (v) = (v12 + v22 + v32 − 3)(v1 − v2 ) = 0(8)
for a suﬃcient small positive value r. The relations
F2 (v) = v2 − v3 = 0
eq.(11) are also solved by the Newton homotopy method.
We have the following modiﬁed relations:
We assume that the values (p1 , p2 , p3 , p4 ) obtained fromthe eq.(11) are the initial values of the branches startingfrom the bifurcation point as shown in Fig.5.

f1 (v) = (v12 + v22 + v32 − 3)(v1 − v2 ) = 0f2 (v) = v2 − v3 = 0
f3 (v, ∆v) = F1 (v1 + ∆v1 , v2 + ∆v2 , v3 + ∆v3 )(9)− F1 (v1 , v2 , v3 ) = 0p1
f4 (v, ∆v) = F2 (v1 + ∆v1 , v2 + ∆v2 , v3 + ∆v3 ) 
p

+

ds

dλ
ds

+

ds

=1 



4

− F2 (v1 , v2 , v3 ) = 0

Note that eq.(8) consists of three variables and two∂F1 (v)=0∂v1

equation, and that it satisﬁes both

r

and

p

∂F1 (v)= 0 at the pitchfork points. Hence, the variables∂v2∆v1 , ∆v2 and ∆v3 can be set as shown in eq.(9). Applying the Newton homotopy method, we have








f1 (v) + (λ − 1)I01 = 0
f2 (v) + (λ − 1)I02 = 0
f3 (v, ∆v) + (λ − 1)I03 = 0
f4 (v, ∆v) + (λ − 1)I04 = 0
3

dvi

2

+

ds

i=1

d∆v3
ds

2

+

dλ
ds

2





=1 


(10)

3

^ )
p.p ( v^ 1 ,v2

p2

Figure 5: Tracing branches from a bifurcation point.

Example 4.1.
We trace the solution curve of the following equation.
f (v1 , v2 ) = (v12 − v2 )(v13 + v12 − v22 ) = 0(12)

We set the ﬁxed small variables ∆v1 and ∆v2 . Thus,
where we chose ∆v1 =0, ∆v2 =0.01.
Then, we we have the following modiﬁed relations:have 2 pitchfork bifurcation points (1.00, 1.00, 1.00) and(−1.00, −1.00, −1.00) at λ = 1 as shown in Fig.4.

f1 (v) = (v12 − v2 )(v13 + v12 − v22 ) = 0
f2 (v, ∆v) = f (v1 + ∆v1 , v2 + ∆v2 ) − f (v1 , v2 ) = 0
v
v
v
(13)
Applying the Newton homotopy method, we have
4

1

3

2

2

ρ=1

v

P1

P2

ρ

0


f1 (v) + (ρ − 1)I01 = 0
f2 (v, ∆v) + (ρ − 1)I02 = 0
2 
2  2

dv2
dρ
dv1


+
+
=1 
ds
ds
ds

-2

-4

-6
0

5

10

15

20

s

(14)

where we chose ∆v1 =0, ∆v2 =0.01. Then, we have 3
pitchfork bifurcation points (0.42 × 10−9, −0.32 × 10−6),(1.62, 2.62) and (−0.62, 0.38) at λ = 1 as shown in
The algorithm is much more simple compared to the Fig.6(a).
Figure 4: Solution of the pitch fork point.

- 453 -

Applying the eq.(11) and Newton homotopy method, Furthermore, we have used Newton homotopy method
to obtain the initial value of the branchies. Thus, we for ﬁnding the multiple solutions of nonlinear equations.
obtain the following equation:
Therefore, we can trace the solution curves with thetransient algorithm of SPICE. After ﬁnding the singular
points, we restart to trace the solutions curves in the dif(v12 − v2 )(v13 + v12 − v22 ) + (ρ − 1)I03 = 0ferent directions. We could trace all the solution curves(v − v̂ )2 + (v − v̂ )2 − r2 + (ρ − 1)I04 = 0 
 1 12  2 2 2  2(15) and ﬁnd all the solutions. But, although it is impossible
dv1
dv2
dρ


to prove whether all the solutions of the nonlinear alge
+
+
=1
ds
ds
ds
braic equations are found with our method or not, all
The results for the above example are shown by the solutions can be found in many practical problems.
In the future problems, we need to extend the algorithm
Fig.6(b),(c) and (d), respectively.
to large dimensional cases.
References
ρ=1
1

v2
v1

2

1

Bifurcation-1

1.05

3

0.95

ρ=1

0.9

ρ

0.85

[1] L.Kolev, “An interval method for global nonlinear analysis,”IEEE Trans. on Circuits and Systems-I, vol.47, pp.675-683,
2000.

ρ

0.8

0

0.75

[2] K.Yamamura and T.Ohshima, “Finding all solutions ofpiecewise-linear resistive circuits using linear programming,”
IEEE Trans. on Circuits and Systems-I, vol.45, pp.434-445,
1998.

-1

0.7
0.65
0.8

-2
6

8

10

12

14

16

18

20

1.2

1.4

1.6

1.8

2

(b) (0.42 × 10−9 ,-0.32 × 10−6 )Bifurcation-2

1.2

Bifurcation-3

1.1

ρ=1

1

2.2

s

s

(a) Bifurcation points.

[3] K.Yamamura, M.Sato, O.Nakamura and T.Kumakura, “Aneﬃcient algorithm for ﬁnding all DC solutions of piecewiselinear circuits,” IEICE Trans. Fundamentals, vol.E85-A,
PP.2459-2468, 2002.

ρ=1

1

0.8

0.9

ρ

0.6

0.8

0.4

ρ

0.7

0.2

0

1

[4] L.Kolev and V.Mladenov, “An interval method for ﬁnding
all operating points of nonlinear resistive circuits,” Int. Jour.
Circuit Theory and Applications, vol.18, pp.257-267, 1990.

0.6

0

1

2

3

4

5

6

7

8

9

0.5

s

0

0.5

(c) (1.62,2.62)

1

s

1.5

2

(d) (-0.62,0.38)

Figure 6: Solving of bifurcation points and initial value.
Starting from these intersections of bifurcation points,
we can trace the solution curves of eq.(12) with arclength method [9].(Fig.7)[5] W.I.Zangwill and C.B.Garcia, Pathways to Solutions, Fixed
Points and Equilibria, Prentice-Hall, Inc. Englewood Cliﬀs,N.J., 1981.
[6] L.T. Watson, “Globally convergent homotopy algorithms fornonlinear systems of equations,”Nonlinear Dynamics, vol.1,
pp.143-191, 1990.
[7] A.Ushida, Y.Yamagami,Y.Nishio, I.Kinouchi and Y.Inoue,“An eﬃcient algorithm for ﬁnding multiple DC solutions
based on SPICE-oriented Newton homotopy method,”IEEETrans. on CAD, vol.21, pp.337-348, 2002.

3

2[8] Y.Inoue, “A practical algorithm for dc operating-point analysis of large scale circuits,”IEICE Trans. Fundamentals,
vol.J77-A, pp.388-398, 1994.

Bifurcation Point-1
1

Bifurcation Point-2

v2

Bifurcation Point-3

0[9] R.C.Melville, Lj.Trajkovic, S.C.Fang and L.T.Watson, “Artiﬁcial parameter homotopy methods for the dc operating
point problem,” IEEE Trans. on CAD, vol.12, pp.861-877,
1993.

-1

-2
-2

-1

0

1

2

3

[10] A.Ushida and L.O.Chua, “Tracing solution curves of nonlinear equations with sharp turning points,” Int. Jour. Circuit
Theory and Applications, vol.12, pp.1-21, 1984.

v1

Figure 7: Solution curves of eq.(12).
5. Conclusions and remarks
It is very important to ﬁnd all the solutions of nonlinear equations. In this paper, we proposed an eﬃcient
Spice-oriented algorithm for ﬁnding bifurcation points,
where we never use the Jacobian matrices and only use
the set of algebraic equations having auxiliary variables.

[11] M.Kubicek and M.Marek, Computational Methods in Bifurcation Theory and Dissipative Structure, Springer-Vrlag,
New-York, 1983.
[12] N.Yamamoto, “Bifurcations of solutions of nonlinear equations involving parameters” Theoretical and Applied Mechanics, University of Tokyo press, vol.33, pp.435-444, 1985.

- 454 -