INTERNATIONAL JOURNAL OF CIRCUIT THEORY AND APPLICATIONS, VOL. 25, 95–105 (1997)

ANALYSIS OF NON-LINEAR TRANSMISSION LINES BYFREQUENCY DOMAIN PERTURBATION METHODAKIO USHIDA, YUICHI TANJI AND YOSHIFUMI NISHIO
Department of Electrical and Electronic Engineering, Faculty of Engineering, Tokushima University, Tokushima 770, JapanSUMMARY
We discuss a numerical method for solving non-linear transmission lines in the frequency domain. Such transmission lines
are widely used for communications such as in GaAs integrated circuits and varactor diode circuits. The circuit equations
are described by non-linear partial di erential equations, so their analysis is very complicated compared with that of linear
transmission lines.
In this paper we propose a frequency-domain perturbation method for weakly non-linear transmission lines where thewave-forms are approximated by Fourier expansions and each frequency component is calculated by a modi ed perturbation
method. To improve convergence, we introduce two new techniques, the compensation method and the homotopy method,which help to make the iteration stable and can be applied to a wide class of non-linear transmission lines. We haveanalysed shock wave phenomena in example. ? 1997 by John Wiley & Sons, Ltd. Int. J. Circ. Theor. Appl., vol. 25,
95–105 (1997)(No. of Figures: 6;

No. of Tables: 0;

No. of Refs: 20.)

1. INTRODUCTION
Recently the analysis and design of high-speed non-linear systems have become more important in the elds ofcommunications and LSI chips. On the other hand, fast-rise-time electrical devices are crucial for high-speed
digital and analogue applications and for wide-bandwidth electronic systems. In GaAs non-linear transmission
lines,1 varactor diode circuits and superconducting transmission lines2 the capacitors have non-linear voltage–
charge characteristics. The velocity of a travelling wave in these non-linear transmission lines is given by afunction of amplitudes as follows:
1dx=√
L(i)C(v)dt
If the capacitor and/or inductor are decreasing functions of voltages and currents, the velocity will increase
with their amplitudes. Therefore it has an interesting property that the higher part of the wave-form is fasterthan the lower part of it. This property is applied to many purposes such as picosecond pulse compression,broadband phase modulation,1 picosecond shock wave generation,3; 4 MESFET gate mixers5 and so on. Since
it is not very easy to analyse non-linear transmission lines, some papers discuss the phenomena from thephysical point of view.1 – 3; 6 – 8
On the other hand, the harmonic balance method is known as a powerful method for solving non-linear
circuits when the wave-form does not contain many higher harmonics. In Reference 9, shock waves in nonlinear transmission lines are analysed by the harmonic balance method, where the line is approximated bya nite number of discrete lumped models driven by a sinusoidal wave generator. Note that the harmonicbalance method can be eciently applied only when the number of discrete models is not too large. The
shock wave is one of the important phenomena in non-linear transmission lines. In Reference 10 it has been
analysed by applying the di erence approximation to the non-linear partial di erential equations. In Reference
11, non-linear wave propagation phenomena are analysed by two-dimensional FFT. On the other hand, device
models of IC chips are also described by non-linear partial di erential equations12; 13 and estimation of the
CCC 0098–9886/97/020095–11$17.50? 1997 by John Wiley & Sons, Ltd.

Received 7 November 1995
Revised 30 April 1996

96

A. USHIDA ET AL.

delay time bound is also an important problem in the design of IC chips. In Reference 12 it is calculatedusing the RC discrete model.
In this paper we discuss an ecient numerical method for solving non-linear transmission lines. We approximate the wave-form by a Fourier expansion at each point on the line and the coecients are calculated
by a perturbation technique. Although the perturbation method is only applicable to weakly non-linear systems, we have improved the property by introducing two techniques. The rst one reduces the perturbed termby introducing a negative compensation component. The second one gradually increases the perturbed term
step-by-step from zero to a speci ed value, which is the same technique as the homotopy method.14 The ideas
will be applicable to the analysis of wide classes of non-linear circuits.
In Section 3 we apply our method to the analysis of shock waves in non-linear transmission lines. We ndfrom the example that the convergence ratio is greatly improved by the introduction of the compensation andhomotopy methods.
2. FREQUENCY DOMAIN PERTURBATION METHOD
2.1. Perturbation method applied to non-linear transmission lines
There are many kinds of non-linear elements such as Schottky barrier, metal–insulator–semiconductor (MIS)and GaAs integrated transmission lines1 – 4; 9 which are used for harmonic and shock wave generation at
microwave lengths. Their structures can be approximately described by a discrete model per small length(dx) in the spatial co-ordinate as shown in Figure 1. Thus the circuit is described by the non-linear partial
di erential equations
−

@L
@v
=
+ vR ;
@x
@t

−

@qC
@i
=
+ iG
@x
@t

(1)

where L and qC are the ux of the non-linear inductor and the charge of the non-linear capacitor respectively.
We assume the non-linear characteristics
iL = L + iˆL (L );vR = RiR + v̂R (iR )(2a)vC = SqC + v̂C (qC );iG = GvG + iˆG (vG )(2b)where  is a small constant and iˆL (L ); v̂R (iR ); v̂C (qC ) and iˆG (vG ) are non-linear terms. From the circuit modelin Figure 1, we have the relations i ≡ iL = iR and v ≡ vC = vG .

Fig. 1. Discrete model of non-linear transmission line
? 1997 by John Wiley & Sons, Ltd.

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

ANALYSIS OF NONLINEAR TRANSMISSION LINES BY FREQUENCY DOMAIN PERTURBATION METHOD

To analyse shock waves, consider the input impulse wave-formEm sin[(2k=T )t]; 06t6T=2ke(t) =
0;
T=2k6t6T

97(3)

where T is the period of the impulse and k is an integer. An example of this wave-form is shown in Figure 5
where T = 1 ns; Em = 0·25 V and k = 16 (see Section 3).
In order to apply the frequency domain approach, we describe it by a Fourier expansion in the complex
domain ase(t) =MP

Ek ejk!t ;! = 2=T

(4)

k=−M

for suciently large M . Let us describe the steady-state responses in Fourier expansions asv(x; t) =MP

Vk (x)e jk!t ;k=−ML (x; t) =

M
P

M
P

i(x; t) =Ik (x)e jk!t

(5a)

k=−M

k (x)e jk!t ;MP

qC (x; t) =k=−M

Qk (x)e jk!t

(5b)

k=−M

Substituting (5) into (1) and (2), we get the relations−

dVk= jk!k + RIk + V̂ R; k (i)
dx

(6a)

dIk
= jk!Qk + GVk + IˆG; k (v)
dx

(6b)

Ik = k + IˆL; k (L )(6c)

Vk = SQk + V̂ C; k (qC )

(6d)

−

where V̂ R; k (i); IˆG; k (v); IˆL; k (L ) and V̂ C; k (qC ) are the kth frequency components obtained from the nonlinear terms. Observe that since the non-linear terms are functions of 2M + 1 Fourier coecients, the solutionsof (6) satisfying the terminal relations are calculated by solving the non-linear two-point boundary value problem in 4(2M + 1) dimensions.15 This is really time-consuming when the frequency components are increased.
Fortunately, if the non-linear terms are small enough compared with the linear terms, the perturbation method16
can be eciently applied to solve (6).
Now let us assume the solutions at the mth iteration asvm (x; t) =

M
P
k=−M

m
L (x; t) =

M
P
k=−M

Vkm (x)e jk!t ;im (x; t) =

jk!t
m
;
k (x)eqCm (x; t) =

M
P
k=−M

Ikm (x)e jk!t

(7a)

Qkm (x)e jk!t

(7b)

M
P
k=−M

Then from (6) and (7) the perturbed equations are written in the formsdVkm= −(R + jk!L)Ikm −  V̂ R; k (im−1 ) − jk!LIˆL; k (m−1
)
L
dx


dIkm
= −(G + jk!C)Vkm −  IˆG; k (vm−1 ) − jk!C V̂ C; k (qCm−1 )dxwhere L = 1= ; C = 1=S and
 m

m−1
ˆ
m
)
k = L Ik − I L; k (L
Qkm = C Vkm − V̂ C; k (qCm−1 )? 1997 by John Wiley & Sons, Ltd.

(8a)
(8b)

(8c)
(8d)

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

98

A. USHIDA ET AL.

Observe that relations (8a) and (8b) are ordinary di erential equations w.r.t. the distance x whose secondmterms on the right-hand side correspond to the forced terms obtained at the (m − 1)th iteration. mk and Qk
can be decided from (8c) and (8d) after solving (8a) and (8b).
The solutions of (8a) and (8b) are written asmk x+ Bkm e−k x + V̂ k (x)Vkm (x) = Am
ke

1
m
k x
+ Bkm e−k x + Iˆk (x)−Am
Ikm (x) =
ke
Zk;0

(9a)
(9b)

where the propagation constant k and characteristic impedance Zk;0 for the kth frequency component are
given byr
√
jk!L + R
Zk;0 =
k = [(jk!L + R)(jk!C + G)];
jk!C + Gmm
The rst two terms of (9) correspond to the zero-input responses; V̂ k (x) and Iˆk (x) are the zero-state17
responses, which can be calculated by a numerical intergration technique such as backward di erence.18
Here the constant parameters Ak and Bk in (9) are evaluated by the boundary conditions of the transmission
line. When it is terminated by the impedance ZL at the far end, they are given as
h m
i


m
 Iˆk (l)ZL (jk!) − V̂ k (l) Zk;0 − Zk;0 − ZL (jk!) Ek e−k l




(10a)
Am
k =
Zk;0 + ZL (jk!) ek l − Zk;0 − ZL (jk!) e−k l
h m
i


m
 −Iˆk (l)ZL (jk!) + V̂ k (l) Zk;0 + Zk;0 + ZL (jk!) Ek ek l




(10b)
Bkm =
Zk;0 + ZL (jk!) ek l − Zk;0 − ZL (jk!) e−k lThus the steady state wave-forms at the mth iteration are estimated by (5a) and (5b).
The iteration will be continued until the variation
!
s
M 
M 
2
2
P
P
m−1
m−1
m
m
Vk (l) − Vk (l) +Ik (l) − Ik (l)
m ≡
k=−M

(11)

k=−M

becomes suciently small for a given constant . Note that we cannot say whether the solution wave-form
is suciently accurate or not. Thus let us de ne the following residual error. Assume that the solutions for
M frequency components are vM (l; t) and iM (l; t) and that the exact solutions are v̂(l; t) and iˆ(l; t). Then theresidual error is de ned asr Z T
1{[vM (l; t) − v̂(l; t)]2 + [iM (l; t) − iˆ(l; t)]2 }dt
(12)
”M ≡
T 0
Since it is impossible to get the exact solutions, we assume vM 0 (x; t) and iM 0 (x; t) as the exact solutions which
are obtained with a larger number of frequency components (M 0  M ). Then we have
!
s
M0
M0
P
P
M
M0
2
M
M0
2
[Vk (l) − Vk (l)] +[Ik (l) − Ik (l)]
”M ≈
|k|=M

|k|=M

If the residual error M is not small enough, we need to choose a much larger M for the approximation andagain try our perturbation method.
2.2. Improvements in perturbation method
It is known that the perturbation method16 can be applied only to weakly non-linear circuits. In thissubsection we improve the method such that it can be applied to a wider class of non-linear circuits.
Assume that the non-linear characteristics given in (2) are monotone increasing functions. For simplicity
we describe the situation as
u1 = Hu2 + ĥ(u2 )? 1997 by John Wiley & Sons, Ltd.

(13)
Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

ANALYSIS OF NONLINEAR TRANSMISSION LINES BY FREQUENCY DOMAIN PERTURBATION METHOD

where u1 = (iL ; vR ; vC ; iG )T ;u2 = (L ; iR ; qC ; vG )T and

H = diag( ; R; S; G);99

Tĥ(u2 ) =  iˆL (L ); v̂R (iR ); v̂C (qC ); iˆG (vG )To improve the convergence ratio of our perturbation method, we try to reduce the perturbed term ĥ(u2 ).
Let us introduce a small constantH ≡ diag( ; R; S; G)
as the compensation parameter in the form
u1 = (H + H )u2 + [ĥ(u2 ) − Hu2 ](14)The second term corresponds to the perturbed term, which can be greatly reduced by choosing a suitableH . We found from our numerical example that the convergence ratio is greatly improved even for a small
compensation parameter H .
The homotopy method14 is sometimes used when the iteration method does not guarantee convergence.
It has the property of global convergence. We will apply it to our perturbation method. Let us introduce aparameter  for the non-linear elements as follows:
u1 = (H + H )u2 + [ĥ(u2 ) − Hu2 ];: 0 −→ 1(15)

It is clear that relation (15) for  = 0 corresponds to the linear case and that the non-linearity is gradually
increased by choosing {: 0 → 1}. At  = 1 it is reduced to the original equation. If we choose a small
variation  per iteration such asm = m−1 + 
our perturbation method can get the solution stably.
2.3. Convergence conditions
Now let us consider the convergence conditions of our perturbation method. For simplicity let us introduce
the following notation for the variables in (8) and (9).
(a) For the solutions we set


 m
 m
Vk (x)
k (x)mm
X1; k (x) =;X2; k (x) =(16a)
Ikm (x)Qkm (x)(b) For the perturbed terms we setV̂ R; k (im−1 )m−1X̂ 1; k (w1 ; x) = ˆ;
I G; k (vm−1 )X̂ 2; k (w2m−1 ; x) =where w1 = (v; i)T and w2 = (L ; qC )T .
We describe the coecient matrices of (8) as0
R + jk!LAk =;G + jk!C
0


0 L
C0 =
C 0


Bk =

IˆL; k (m−1)
L
m−1
V̂ C; k (qC )jk!L 00 jk!C


(16b)


(17a)
(17b)

Now we describe the zero-state response of (9) in the integrations. Then we have from (9) and (10)Z xmme−Ak (x−s) X̂ 1; k (w1m−1 ; s)ds
X1; k (x) = C1; k (x)Ek + C2; k (x)X1; k (l) − 0Z xe−Ak (x−s) X̂ 2; k (w2m−1 ; s)ds
(18a)
−Bk
0

X2;mk (x) = C0 X1;mk (x) − C0 X̂ 2; k (w2m−1 ; x)(18b)where C1; k (x) and C2; k (x) are obtained from (9) and (10).
? 1997 by John Wiley & Sons, Ltd.

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

100

A. USHIDA ET AL.

Let us adopt the norm19rkw(t)k ≡

Z T

1
T


w2 (t)dt(19)

0Assume that w(t) is a periodic function of the period T as follows:MP

w(t) =

Wk e jk!tk=−M

Then we have the Euclidean norm
kW k ≡√|W−M |2 + |W−M +1 |2 + · · · + |WM |2



(20)

Thus the two norms coincide with each other, i.e.
k w(t) k = k W k
Now we have the following convergence conditions.
Theorem 1.
Assume that for an approximate solution (w10 ; w20 )T there are constants D1 ; D2 ; K0 ; L1 and L2 satisfying thefollowing:(i) De ne the solution domain by
(21a)≡ w1 ; w2 | k w1 − w10 k6D1 ; kw2 − w20 k 6D2(ii) De ne the norms of the coecient matrices by
K0 ≡ k C0 k;

K1 ≡

k A−1
k k L1 l;−M 6k6M
1 − C2kmax

K2 ≡k Bk A−1
k k L2 l−M 6k6M
1 − C2kmax(21b)

k X̂ 1 (w1 0 ; x) − X̂ 1 (w1 00 ; x) k 6L1 k w1 0 − w1 00 k

(21c)

k X̂ 2 (w2 0 ; x) − X̂ 2 (w2 00 ; x) k 6L2 k w2 0 − w2 00 k

(21d)

where
C2k ≡ max k C2; k (x) k06x6l(iii) Assume that the perturbed terms satisfy the Lipschitz conditions:infor all 06x6l, where L1 and L2 are Lipschitz constants.
(iv) Set the maximum values of the variables in x = [0; l] as
X1;mk = max k X1;mk (x) k;

X2;mk = max k X2;mk (x) k

06x6l

Furthermore, if it satis es the condition ≡  k P k¡ 1;

06x6l


for P =K2
K1
K0 K1 K0 (K2 + L2 )

(21e)


(22)

our perturbation method (18) will converge to the unique solution (w1 ; w2 )T .
Proof. For non-linear transmission lines we can assume k e−Ak x k 61 for 06x6l. Hence we have from(18a) and (21b)17m−1m
k X1;mk (x) − X1;m−1
k (x) k 6 C2k k X1; k (l) − X1; k (l) k? 1997 by John Wiley & Sons, Ltd.

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

101

ANALYSIS OF NONLINEAR TRANSMISSION LINES BY FREQUENCY DOMAIN PERTURBATION METHOD

Z x+ k A−1
k k

0+ k Bk A−1
k k

k X̂ 1; k (w1m−1 ; s) − X̂ 1; k (w1m−2 ; s) k ds

Z x0

k X̂ 2; k (w2m−1 ; s) − X̂ 2; k (w2m−2 ; s) k ds

Applying the Lipschitz conditions (21c) and (21d), we have
M
P
k=−M

k X1;mk − X1;m−1
k k 6K1

M
P
k=−M

m−2
k X1;m−1
k − X1; k k +K2

M
P
k=−M

m−2
k X2;m−1
k − X2; k kThus we have from (20)k X1m − X1m−1 k 6K1 k X1m−1 − X1m−2 k +K2 k X2m−1 − X2m−2 k(23a)

In the same manner we have from (18b)k X2m − X2m−1 k 6K0 k X1m − X1m−1 k +K0 L2 k X2m−1 − X2m−2 k(23b)

Set X = (X1 ; X2 )T . Then we have from (22) and (23)k X m − X m−1 k 6 k X m−1 − X m−2 k

(24)

Hence
k X m − X m−1 k 6m−1 k X 1 − X 0 k
Thus our perturbation method will converge to the solution X and the error bound of the approximate solution
X 0 is estimated ask X m − X 0 k 6 k X m − X m−1 k + k X m−1 − X m−2 k + · · · + k X 1 − X 0 k
6 (m−1 + m−2 + · · · +  + 1) k X 1 − X 0 k
1 − m=k X1 − X0 k
1−Thus we have proved the convergence condition.

Q.E.D

We found from the Theorem 1 that our perturbation method will converge to an approximate solution X if
the non-linear term  is suciently small.

3. ILLUSTRATIVE EXAMPLEAs an example of non-linear transmission lines, let us calculate the shock wave due to an impulse. Assume
that the line is terminated by a resistor RL = 10 . For simplicity, we assume that the inductor is linear andthe other characteristics in (2) are given as
iL = 350 × 1011 L ;? 1997 by John Wiley & Sons, Ltd.

vR = 0·015iR + 0·1iR3 ;vC = 1011 qC + 0·5 × 1034 qC3 ;3iG = 0·015vG + 0·1vG

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

102

A. USHIDA ET AL.

Fig. 2. Characteristic curve of non-linear capacitor: vC = 1011 qC + 0·5 × 1034 qC3

Fig. 3. Convergence ratios for various compensation parameters, where variations are de ned by m =k X m (l) − X m−1 (l) k

The q – v characteristic curve is shown in Figure 2. We assume Em = 0·25 V; T = 1 ns and k = 16 for the input
impulse in (3). For the calculation of zero-state responses in (8) we applied the rst- and/or second-orderbackward di erence formulae18 , where the step size is chosen as x = l=60.
? 1997 by John Wiley & Sons, Ltd.

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

ANALYSIS OF NONLINEAR TRANSMISSION LINES BY FREQUENCY DOMAIN PERTURBATION METHOD

103

Fig. 4. Perturbed term by compensation method for  = 0·1. v̂C = 0·5 × 1034 qC3 − 0·1 × 1011 qC

The convergence ratios are shown in Figure 3 for various compensation parameters H . In this case thevariations are estimated by
128 
P

m =

k=−1282
Vkm (l) − Vkm−1 (l) +

128 
2
P
Ikm (l) − Ikm−1 (l)!1=2

k=−128

We found from the results that the convergence will be greatly improved even for a small  in (14). Figure 4
shows the perturbed term for  = 0·1:v̂C = 0·5 × 1034 qC3 − 0·1 × 1011 qC
When we applied both the compensation method with  = 0·03 and the homotopy method withm = m−1 + 0·1;

for m610;

=1

for m ¿ 10

our method could get the solution stably. The wave-forms are shown in Figures 5 and 6, where we assumed
128 frequency components. Wave-form (a) in Figure 5 is the response of the non-linear transmission line,
while wave-form (b) shows the linear response neglecting non-linear terms. The features of the shock wave
are apparent in the non-linear response. Wave-form (c) is the input impulse wave-form. The residual error
de ned by (12) for M 0 = 256 is given by128 = 0·97 × 10−2
The error in the amplitude is about 2·4%, which seems to be accurate enough.
We remark that although the perturbation method never converges for the stronger non-linear term
vC = 1011 qC + 0·8 × 1034 qC3(25)

with  = 0, our perturbation method with  = 0·1 can get the solution stably.
Thus we found from the example that our perturbation method is greatly improved by incorporating thetwo techniques of compensation and homotopy.
? 1997 by John Wiley & Sons, Ltd.

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

104

A. USHIDA ET AL.

Fig. 5. Impulse responses: (a) response of non-linear transmission line; (b) response of linear transmission line neglecting non-linear
terms; (c) input impulse, (Em = 0·25 V; T = 1 ns, and k = 16 in (3))

Fig. 6. Response of non-linear transmission line in (x; t) domain (l = 10 cm; T = 1 ns, Em = 0·25 V, RL = 10 )4. CONCLUSIONS AND REMARKS
We have presented a frequency domain perturbation algorithm for calculating transient responses of non-linear
transmission lines. We can greatly improve the convergence by introducing two techniques, the compensation method and the homotopy method. They have the properties of improving the convergence ratios andstabilizing the iteration by reducing the perturbed terms. Thus our perturbation method becomes much more
powerful for the analysis of non-linear transmission lines.
Note that since the algorithm in this paper is based on the perturbation method, convergence may not
be guaranteed even if we choose any  and  for strong non-linearity. The non-linearity of the element depends on the amplitude of the input variable and can be estimated by the ratio (non-linear term)/? 1997 by John Wiley & Sons, Ltd.

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)

ANALYSIS OF NONLINEAR TRANSMISSION LINES BY FREQUENCY DOMAIN PERTURBATION METHOD

105(linear term). In the example of (2) we have about 0·5 at vC = 0·4. We found from our numerical experiencesthat we cannot hope for convergence when the ratio is greater than unity.
Most practical circuits such as GaAs transmission lines1 belong to the class of non-linear RC transmission
lines. They are special cases of our example where our algorithm can be eciently applied to much stronger
non-linearity.
REFERENCES
1. M. J. W. Rodwell, ‘Nonlinear transmission line for picosecond pulse compression and broadband phase modulation’, Electron. Lett.,
23, 109–110 (1987).
2. R. K. Arora and J. L. Thaker, ‘Electromagnetic pulse transmission on superconducting interconnects’, Int. J. Electron., 71, 117–125(1991).
3. C. J. Madden, R. A. Marsland, M. J. W. Rodwell, D. M. Bloom and Y. C. Pao, ‘Hyperabrupt-doped GaAs nonlinear transmission
line for picosecond shock wave generation’, Appl. Phys. Lett., 54, 1019–1021 (1989).
4. K. S. Champlin and D. R. Singh, ‘Small-signal second-harmonic generation by a non-linear transmission line’, IEEE Trans.
Microwave Theory Tech., MTT-34, 351–353 (1986).
5. C. Camacho-Penalosa and C. S. Aitchison, ‘Analysis and design of MESFET gate mixers’, IEEE Trans. Microwave Theory Tech.,
MTT-35, 643–652 (1987).
6. F. A. Benson and J. D. Last, ‘Non-linear-transmission-line harmonic generator’, Proc. IEE, 112, 635–643 (1965).
7. R. Landauer, ‘Shock waves in nonlinear transmission lines and their e ect on parametric ampli cation’, IBM J., October, pp.
391– 401 (1960).
8. D. Jager and F. -J. Tegude, ‘Non-linear wave propagation along periodic-loaded transmission line’, Appl. Phys., 15, 393–397 (1978).
9. C. Camacho-Penalosa and I. Molina-Fernandez, ‘Harmonic balance analysis of non-linear transmission lines’, Electron. Lett., 24,1235–1236 (1988).
10. R. H. Freeman and A. E. Karbowiak, ‘An investigation of nonlinear transmission lines and shock waves’, J. Phys. D: Appl. Phys.,
10, 633–643 (1977).
11. I. Molina-Fernandez, C. Camacho-Penalosa and J. Ramos, ‘Application of the two-dimensional Fourier transform to non linear wave
propagation phenomena’, IEEE Trans. Microwave Theory Tech., 42, 1079–1085 (1994).
12. C. A. Marinov and P. Peittaanmaki, ‘A delay time bound for distributed parameter circuits with bipolar transistors’, Int. j. cir. theor.
appl., 18, 99–106 (1990).
13. J. L. Wyatt Jr., ‘Monoton sensitivity of nonlinear nonuniform RC transmission lines with application to timing analysis of digital
MOS integrated circuits’, IEEE Trans. Circuits and Systems, CAS-32, 28–83 (1985).
14. C. B. Garcia and W. I. Zangwill, Pathways to Solutions, Fixed Points, and Equilibria, Prentice-Hall, Englewood Cli s, NJ, 1981.
15. S. M. Roberts and J. S. Shipman, Two-Point Boundary Value Problems: Shooting Methods, American Elsevier, New York, 1972.
16. J. Kevorkian and J. D. Cole, Perturbation Methods in Applied Mathematics, Springer, New York, 1981.
17. L. A. Zadeh and C. A. Desoer, Linear System Theory: the State Space Approach, McGraw-Hill, New York, 1963.
18. L. O. Chua and P.-M. Lin, Computer-Aided Analysis of Electrical Circuits: Algorithms and Computational Techniques, PrenticeHall, Englewood Cli s, NJ, 1975.
19. M. Urabe,‘Galerkin’s procedure for nonlinear periodic systems’, Arch. Rat. Mech. Anal., 20, 120–152 (1965).
20. A. Ushida, T. Adachi and L. O. Chua, ‘Steady-state analysis of non-linear circuits based on hybrid method’, IEEE Trans. Circuits
and Systems—I, CAS 39, 649–661 (1992).

? 1997 by John Wiley & Sons, Ltd.

Int. J. Circ. Theor. Appl., Vol. 25, 95–105 (1997)