July 4, 2003 8:40 WSPC/103-M3AS
00290
Mathematical Models and Methods in Applied Sciences Vol. 13, No. 9 (2003) 1–18 c World Scientific Publishing Company
ELLIPTIC PARTIAL DIFFERENTIAL-ALGEBRAIC MULTIPHYSICS MODELS IN ELECTRICAL NETWORK DESIGN
G. AL`I Istituto per le Applicazioni del Calcolo “M. Picone”, Via Pietro Castellino 111, 80131 Naples, Italy
[email protected] ¨ A. BARTEL and M. GUNTHER Institut f¨ ur Wissenschaftliches Rechnen und Mathematische Modellbildung, Universit¨ at Karlsruhe (TH), Germany C. TISCHENDORF Institut f¨ ur Mathematik, Humboldt-Universit¨ at zu Berlin, Germany Received 29 April 2002 Revised 21 November 2002 Communicated by P. A. Markowich In refined network analysis, a compact network model is combined with distributed models for semiconductor devices in a multiphysics approach. For linear RLC networks containing diodes as distributed devices, we construct a mathematical model that combines the differential-algebraic network equations of the circuit with elliptic boundary value problems modeling the diodes. For this mixed initial-boundary value problem of partial differential-algebraic equations a first existence result is given. Keywords:
1. Introduction Spatial dimensions of semiconductor technology shrink steadily, so that former secondary effects tend to influence or even conduct the general behavior of electric circuits. Hence full models for semiconductor devices coupled to electric circuit equations have to be investigated. In this paper, we concentrate on RLC networks with bipolar semiconductor devices, such as diodes. Due to the very different time scales related to the relaxation of diodes to equilibrium and to the electric current in the network, it is appropriate to model the devices by stationary drift-diffusion equations. For the sake of simplicity, we derive appropriate coupling conditions for one-dimensional diodes and linear RLC networks, set up by Modified Nodal Analysis (MNA). This multiphysics approach yields a coupled system of (elliptic) partial differential 1
July 4, 2003 8:40 WSPC/103-M3AS
2
00290
G. Al`ı et al.
equations (PDEs) and differential-algebraic equations (DAEs), for short, a system of partial differential-algebraic equations (PDAEs): the node potentials of the network define boundary conditions for the diode model and this causes each diode to produce a current flow, such that each diode acts as a voltage-defined current source for the electric network. In this paper we aim at proving the existence of solutions to this PDAE. As a by-product of the existence proof, we show that the perturbation index of the coupled system is one. Although, from a physical point of view, we expect uniqueness of solutions, we are not able to prove such a result at the moment. Nevertheless, in the framework of the present paper, we believe that it is possible to prove a uniqueness result following the analysis of the stationary drift-diffusion model7 at least when the boundary conditions for the diodes remain close to equilibrium. The drift-diffusion model and the RLC network equations, considered as non-coupled systems, have been object of a conspicuous number of papers and books. From a mathematical point of view, the drift-diffusion equations have been extensively studied with reference to the existence and uniqueness of solutions, and to discretization procedures. Existence results have been established for several physical models of mobility, diffusivity and generation-recombination rate. In general, uniqueness can be guaranteed only around equilibrium. Non-uniqueness of solutions is sometimes physically required by specific devices. For these mathematical oriented topics, and for their relevance to numerical procedures, we refer to the books5,7,8,10 of the bibliography and to the references therein. The network approach, which aims at an automatic assembly of a mathematical model, has been established in computer-aided analysis of complex technical systems of various kinds. In that way, a set of redundant variables is inevitable and, indeed, for circuit simulation the first DAE solver has been developed by Gear;1 many theoretical investigations followed. For the theory of DAEs, we refer to the book of Hairer and Wanner6 (and references therein) and for the set up of network equations to the MNA.4,9 The concept of perturbation index is crucial for the numerical analysis.2,6,12 It describes the sensitivity of the solution to data. Therefore it is a measurement of the numerical difficulties which have to be expected and resolved. The main numerical aspects and problems of the coupled system will be addressed in a forthcoming paper. To the authors’ knowledge, there is only one appearance in literature of an existence analysis of a coupled system of DAEs and PDEs, i.e. RLC networks containing uniform lossy transmission lines.2,3 The latter are described by the telegrapher’s equations, which form a system of linear hyperbolic equations. The model presented in this paper exhibits the same coupling structure as the network with transmission lines: currents generated by the diodes act as source to the network equation, and boundary conditions for the elliptic system are posed by the node potential of the electric network. But, it differs from the previous PDAEs in the respect of type (elliptic instead of hyperbolic) and nonlinearity. In particular,
July 4, 2003 8:40 WSPC/103-M3AS
00290
Elliptic Partial Differential-Algebraic Multiphysics Models
3
the existence proof requires a different approach, namely fixed point arguments instead of Galerkin-type methods. We stress that this is the first mathematical result concerning coupled DAEs and nonlinear PDEs. This paper is organized as follows: Sec. 2 covers the modeling of the coupled system; both subsystems are described in details and, in particular, the coupling terms are defined. The following two sections (Secs. 3 and 4) are devoted to the analysis. Specifically, we state and prove the existence result by applying Schauder’s fixed point theorem in a non-standard fashion. Finally, in the last section we draw some conclusions and give an outlook on a few open problems. 2. Modeling of Diodes and Networks First, we summarize the set of equations for the diode (drift-diffusion)7,8 and the electric circuit (network equations)4,9 and we specify the mathematical coupling terms. Then we present the full coupled system. 2.1. The one-dimensional diode model We model the diode as a line segment of length l, characterized by a doping profile N (x), x ∈ (0, l). We neglect all thermal effects, and assume that two carriers are responsible for the diode’s output current, i.e. electrons with negative charge −q, and holes with positive charge q. The behavior of the diode is described in terms of number densities of electrons and holes, denoted by n = n(x, t), p = p(x, t), current densities for electrons and holes, denoted by jn = jn (x, t), jp = jp (x, t), and electrostatic potential, denoted by V = V (x, t). These variables satisfy the following drift-diffusion system,7 ∂x jn = qR ,
(2.1a)
∂x jp = −qR ,
(2.1b)
jn = −q{−Dn∂x n + µn n∂x V } ,
(2.1c)
jp = q{−Dp ∂x p − µp p∂x V } ,
(2.1d)
−ε∂x2 V = q(N + p − n) ,
(2.1e)
where (x, t) ∈ (0, l) × (0, ∞). The time dependency accounts for the dynamics of the electric network containing the diode and affecting its boundary condition, as we will see. In (2.1), Dn , Dp are the diffusivities and µn , µp are the mobilities for electrons and holes, respectively. They are bounded, strictly positive functions depending on x, n, p and ∂x V . The dependency on the particle densities is usually neglected in most physical models.7,10 Around equilibrium, diffusivities and mobilities are linked by the Einstein relations, Dn = U T µ n ,
Dp = U T µ p ,
(2.2)
July 4, 2003 8:40 WSPC/103-M3AS
4
00290
G. Al`ı et al.
where UT is the thermal potential. Since we neglect thermal effects, we assume at once the validity of the relations (2.2). The generation-recombination term R = ˜ p, jn , jp ) can be modeled as R(n, ˜ p, jn , jp ) = RSRH (n, p) + RAU (n, p) + RII (jn , jp ) , R(n, using Shockley–Read–Hall RSRH and Auger RAU generation-recombination, and the impact ionization rate RII , which is the emission rate of electrons and holes. These terms are given by RSRH (n, p) =
np − n2i , τp (n + ni ) + τn (p + ni )
RAU (n, p) = (Ccn n + Ccp p)(np − n2i ) , RII (jn , jp ) = −αn
|jn | |jp | − αp , q q
where ni is the intrinsic carrier concentration, τp and τn are the lifetimes of electrons and holes, respectively, Ccn and Ccp are the Auger constants, and αn and αp are temperature-related quantities.10 In the following we neglect the impact ionization rate, and assume αn = αp = 0. More generally, we assume R as a function of the type ˜ p, V ) = F˜ (n, p, V ) · np − 1 , (2.3) R(n, n2i with F˜ ≥ 0. It is convenient to introduce the quasi-Fermi potentials n p φn = V − UT ln , φp = V + UT ln . ni ni Then Eqs. (2.1c) and (2.1d) can be written as jn = −qµn n∂x φn ,
(2.4a)
jp = −qµp p∂x φp ,
(2.4b)
where n = ni exp
V − φn UT
,
V − φp p = ni exp − UT
.
The generation-recombination term becomes R = R(φn , φp , V ), with φp − φ n −1 R(φn , φp , V ) = F (φn , φp , V ) · exp UT
(2.5)
(2.6)
(by a slight abuse of notation). The system (2.1) is supplemented with the boundary conditions n(0, t) = n1 ,
p(0, t) = p1 ,
(2.7a)
n(l, t) = n2 ,
p(l, t) = p2 ,
(2.7b)
July 4, 2003 8:40 WSPC/103-M3AS
00290
Elliptic Partial Differential-Algebraic Multiphysics Models
5
φn (0, t) = φp (0, t) = u1 (t) ,
(2.7c)
φn (l, t) = φp (l, t) = u2 (t) ,
(2.7d)
with q q 1 1 N (0) + N (0)2 + 4n2i , p1 = − N (0) + N (0)2 + 4n2i , 2 2 q q 1 1 N (l) + N (l)2 + 4n2i , p2 = − N (l) + N (l)2 + 4n2i . n2 = 2 2 The conditions (2.7a) and (2.7b) come from the requirement that the semiconductor is in equilibrium at the boundary, i.e. n1 =
n(x)p(x) = n2i ,
N (x) + p(x) − n(x) = 0 ,
x = 0, l .
The potentials u1 , u2 in (2.7c), (2.7d) represent the external electric potentials applied to the devices. They are not independent functions to be assigned, but they are to be determined by the equations for the electrical network. The conditions (2.7c) and (2.7d) imply the following conditions for the electric potential, V (0, t) = Vbi (0) + u1 (t) , V (l, t) = Vbi (l) + u2 (t) , where the built-in potential is defined by s 2 N (x) N (x) + + 1 . Vbi (x) = UT ln 2ni 2ni For later use, we introduce the electric currents ! ! jn (0, t) + jp (0, t) λ1 (t) , = −jn (l, t) − jp (l, t) λ2 (t)
(2.8)
which represent the diode’s output to the applied potentials u1 and u2 . In the case of a one-dimensional diode model, we have def
λ1 = −λ2 = I .
(2.9)
In general, if the circuit comprises d diodes, the above introduced scalars have to be replaced by the corresponding vectors: u1 (t), u2 (t) ∈ Rd , and λ(t) = (λ1 (t), λ2 (t))T = (I(t), −I(t))T ∈ R2d . The drift-diffusion equations (2.1) are coupled to the electric network by ! u1 (t) = AT (2.10) λ u(t) , u2 (t) where u(t) ∈ Rn are the node potentials of the DAE model for the electric network, and Aλ is an appropriate incidence matrix, which describes the topology of the coupling (see below). In the next section, we discuss the coupling of the input potential (u1 , u2 )T and the output current λ = (I, −I)T , due to the electric network containing the diodes.
July 4, 2003 8:40 WSPC/103-M3AS
6
00290
G. Al`ı et al.
2.2. Network models for electric circuits An electric circuit is modeled as a network of basic elements: by combining Kirchhoff’s laws and characteristic equations for the circuit elements, the network approach (which aims at preserving the topological structure of the network, and does not aim at models with a minimal set of unknowns) yields a system of differential-algebraic equations (DAEs) in a redundant set of network variables. In the following we consider a linear RLC network with diodes, i.e. our circuit is composed of diodes, linear capacitors, inductors and resistors, and independent voltage and current sources, v(t) ∈ RnV and ı(t) ∈ RnI , only. In the classical Modified Nodal Analysis4,9 (MNA) the vector of unknowns x comprises all node potentials u(t) ∈ Rn , and the currents L (t) ∈ RnL and V (t) ∈ RnV through inductors and voltage sources, respectively: thus x = (u, L , V )T , for t ∈ [t0 , t1 ], say. Rendering the coupling idea for RLC networks and interconnects,3 this set of unknowns is supplemented by currents λ ∈ R2d at the boundaries of all d diodes. It is precisely through λ that the diodes are coupled as non-basic elements to the RLC network. Then the DAE network equation for the RLC part is given by AR GAT AL AV AC CAT 0 0 R C dx + −AT 0 0 x 0 L 0 L dt −AT 0 0 0 0 0 V AI ı Aλ λ (2.11a) + 0 + 0 = 0, 0
v
with consistent initial data
x(t0 ) = x0 .
(2.11b)
The consistency of the initial data (2.11b) with Eq. (2.11a) will be discussed later. The capacitance, inductance and conductance matrices C ∈ RnC ×nC , L ∈ RnL ×nL and G ∈ RnG ×nG are assumed to be positive-definite and symmetric. The incidence matrices AC ∈ Rn×nC , AL ∈ Rn×nL and AR ∈ Rn×nG describe the branchnode relationships for capacitors, inductors and resistors; the incidence matrices AV ∈ Rn×nV and AI ∈ Rn×nI describe this relationship for voltage and current sources, respectively. Finally, Aλ ∈ Rn×2d matches the diodes’ boundaries with the corresponding network nodes. Using a well-established procedure,12 we denote by QC a projector onto the kernel of AT C , and set PC = Id − QC , such that PC QC = QC PC = 0. Then def
the network variables can be split into a differential component, y = (y1 , y2 )T = def
(PC u, L )T , and an algebraic component, z = (z1 , z2 )T = (QC u, V )T . In terms of these components, the network equations read
July 4, 2003 8:40 WSPC/103-M3AS
00290
Elliptic Partial Differential-Algebraic Multiphysics Models
H
0
0
L +
!
y1
d dt
y2
!
+
AR GAT R
AV
−AT L
0
AR GAT R
AL
−AT L
0
!
T QT C AR GAR QC
QT C AV
AT V QC
0
+
QT C (AI ı + Aλ λ) −v
z1 z2 !
!
!
+
z1 z2
!
!
y1 y2
!
(AI ı + Aλ λ) 0 +
7
!
= 0,
(2.12a)
T QT C AR GAR PC
QT C AL
AT V PC
0
= 0.
!
y1 y2
!
(2.12b)
T Here H = AC CAT C + QC QC is a positive-definite, symmetric matrix. This system can be written in the compact form
AP
dy + BP y + CP z + FP (ı, λ) = 0 , dt
(2.13a)
BQ z + CQ y + FQ (ı, v, λ) = 0 ,
(2.13b)
with obvious notation. The functions FP and FQ depend nonlinearly on u = y1 +z1 through the terms Aλ λ and QT C Aλ λ. In fact, the current vector λ is a function of u, through the diode equations. AT λ Next, we discuss the consistency of the initial data (2.11b). We write x0 = (u0 , L0 , V 0 )T and define T y0 = (PT C u0 , L0 ) ,
T z0 = (QT C u0 , V 0 ) .
We supplement Eq. (2.13a) with the initial data y(t0 ) = y0 .
(2.14)
Then the algebraic part z0 of the initial data (2.11b) must satisfy the consistency condition BQ z0 + CQ y0 + FQ (ı(t0 ), v(t0 ), λ0 ) = 0 .
(2.15)
We assume the following topological index-1 conditions to hold,3,12 ker(AC , AR , AV )T = {0} ,
(2.16)
ker QT C AV = {0} .
(2.17)
Then, for the matrix BQ =
T QT C AR GAR QC
QT C AV
AT V QC
0
!
,
July 4, 2003 8:40 WSPC/103-M3AS
8
00290
G. Al`ı et al.
due to the positive definiteness of G, we have ker BQ = ker((AC , AR , AV )T QC ) × ker QT C AV = ker QC × {0} , and z is given as a linear function of y, ı, v and λ, since, by definition, QC z1 = z1 . Moreover, λ does not depend on z provided that the following topological condition holds: AT λ QC = 0 ,
(2.18)
i.e. any diode terminal is connected to the ground by a path of capacitors. In turn, condition (2.18) implies that z is given as a linear function of y, ı and v. Remark 1. Recalling (2.9) and its extension to d diodes, we have λ = (Id, −Id)T I. Then, the topological condition (2.18) can be relaxed by requiring that FQ does not depend on λ, i.e. T QT C Aλ (Id, −Id) = 0 .
(2.19)
Electrically, this means that we require the terminals of each diode to be connected by a path of capacitors. All the results in this paper can be easily generalized after replacing (2.18) with (2.19). For simplicity, we will only refer to (2.18). 2.3. The coupled system Without loss of generality, we consider a network which contains exactly one diode. All the arguments below generalize in a straightforward way to networks containing d diodes. Combining the stationary drift-diffusion model and the RLC-network equations discussed in the previous subsections, we can formulate the full coupled problem, which is summarized in Box 1. Explicitly, Eqs. (2.20) are derived from (2.1), (2.7), (2.8), (2.10), (2.11) in the following way: using dn (x, n, ∂x V ) = µn (x, n, ∂x V ) · n and dp (x, p, ∂x V ) = µp (x, p, ∂x V ) · p, Eqs. (2.1) together with (2.4a), (2.4b) become system (2.20a) with x varying in the interval (0, l); recalling (2.5) and (2.6), n, p and R are regarded as a functions of V , φn and φp ; the elliptic BVP is supplemented with boundary conditions (2.20b); at last, the DAE-IVP (2.20d), (2.20e) is characterized by a symmetric, positive-definite matrix AP , where ı, v are given functions. The coupling structure is as follows. On the one hand, the boundary data u1 , u2 in (2.20b) is identified via y1 = PC u with the corresponding node potentials of the electric network equations (2.20c). On the other hand, the branch current λ entering the network equations (2.20d) is defined by the elliptic BVP via (2.20c). In conclusion, for given initial data y0 and independent current and voltage sources ı(t) and v(t), we wish to find functions φn (x, t), φp (x, t), V (x, t) and y(t) which solve the elliptic partial differential-algebraic system (2.20). A positive answer to this problem is given in the following two sections.
July 4, 2003 8:40 WSPC/103-M3AS
00290
Elliptic Partial Differential-Algebraic Multiphysics Models
9
Mixed system of elliptic BVP and DAE-IVP.
Box 1 Elliptic BVP for the diode:
∂x (dn ∂x φn ) = −R , ∂x (dp ∂x φp ) = R ,
(2.20a)
ε 2 ∂x V = n − p − N , with q φp − φ n −1 , R = F (φn , φp , V ) · exp UT V − φp V − φn n = ni · exp , p = ni · exp − ; UT UT φn (0, t) = φp (0, t) = u1 (t) , φn (l, t) = φp (l, t) = u2 (t) , Vbi = UT ln Coupling interface: ! u1 = AT λ y1 , u2
nbi , ni
λ=
V (0, t) = Vbi (0) + u1 (t), V (l, t) = Vbi (l) + u2 (t) , with q nbi = N/2 + (N/2)2 + n2i ;
−[qdn ∂x φn + qdp ∂x φp ]|x=0 [qdn ∂x φn + qdp ∂x φp ]|x=l
!
;
(2.20b)
(2.20c)
DAE-IVP problem for the network: AP
dy + BP y + CP z + FP (ı, λ) = 0 , dt BQ z + CQ y + FQ (ı, v) = 0 , y(t0 ) = y0 .
(2.20d) (2.20e)
3. Existence of Solutions In this section we prove the existence of a solution to the coupled problem (2.20), by using a fixed point argument. Before entering into the details of the proof, in the first subsection we prove a priori estimates for the network equations, which gives us some necessary restrictions for the network variables. This estimate permits to define a subset of an appropriate function space, where the sought solution may be found. The existence theorem is stated and proved in the subsequent subsection. Finally, we discuss the perturbation index of the system. 3.1. A priori estimate for network variables For any finite-dimensional vector space Rm , let | · | denote the Euclidean vector norm, i.e. for w ∈ Rm holds |w|2 = wT w ≡ w · w.
July 4, 2003 8:40 WSPC/103-M3AS
10
00290
G. Al`ı et al.
Lemma 1. Let C, L, G be symmetric, positive-definite matrices, and let the topological conditions (2.16)–(2.18) be satisfied. Furthermore, let x = (u, L , V )T ∈ (C([t0 , t1 ]))n+nL +nV be a solution to the network equation (2.11a), with consistent initial value x0 = (u0 , L0 , V 0 )T . We assume that ı ∈ (L2 ([t0 , t1 ]))nI , v ∈ (L2 ([t0 , t1 ]))nV , and that λ satisfies the condition T (AT λ u) λ ≥ 0 .
(3.1)
Then, for all t ∈ [t0 , t1 ], the differential part y = (PC u, L )T and the algebraic part z = (QC u, V )T of the solution satisfy the estimates |y|2 (t) ≤ Cy ec1 (t−t0 ) (|y0 |2 + kık2(L2 ([t0 ,t1 ]))nI + kvk2(L2 ([t0 ,t1 ]))nV ) ,
(3.2)
|z|2 (t) ≤ Cz (|y|2 (t) + |ı|2 (t) + |v|2 (t)) ,
(3.3)
for some positive constants Cy , Cz and c1 . Proof. We multiply (2.11a) from the left by xT . After some algebra, we obtain 1 d T T T T {u AC CAT C u + L LL } + u AR GAR u 2 dt + u T Aλ λ + u T AI ı + v T V = 0 .
(3.4)
Recalling the definition of y and z, we have T uT AC CAT C u = y1 Hy 1 .
(3.5)
Since the matrices H and L are symmetric and positive definite, there exist two positive constants cy ≤ Cy such that T 2 cy |y|2 ≤ uT AC CAT C u + L LL ≤ Cy |y| .
(3.6)
Using the positive definiteness of G, and the Schwarz inequality, (3.4) leads to the estimate Z t uT (s)Aλ λ(s) ds cy |y|2 (t) + 2 t0
≤ Cy |y|2 (t0 ) + +
Z
t t0
Z
t t0
2 2 {kAT I k∗ |ı| + |v| }(s) ds
2 T 2 2 {kAT I k∗ |y1 | + kAI k∗ |z1 | + |z2 | }(s) ds ,
(3.7)
where k · k∗ denotes the operator norm. Next, using the topological conditions (2.16)–(2.18), the algebraic system (2.12b) can be solved for z in terms of y, ı and v. Therefore, there exists a constant Cz such that |z|2 ≤ Cz (|y|2 + |ı|2 + |v|2 ) .
(3.8)
July 4, 2003 8:40 WSPC/103-M3AS
00290
Elliptic Partial Differential-Algebraic Multiphysics Models
We fix a time t, with t0 < t < t1 . Using (3.8) into (3.7), we find Z t Z t 2 T cy |y| (t) + u (s)Aλ λ(s) ds ≤ c1 |y|2 (s) ds + d1 , t0
11
(3.9)
t0
where c1 and d1 are constants, with d1 = d1 (y0 , ı, v, t0 , t1 ) given by Z t1 2 d1 = Cy |y0 | + c2 {|ı|2 + |v|2 }(s) ds . t0
Using (3.1) in (3.9), we find cy |y|2 (t) ≤ c1
Z
t
|y|2 (s) ds + d1 .
(3.10)
t0
We can now apply Gronwall’s lemma, and obtain the a priori estimate |y|2 (t) ≤ d1 /cy · e(t−t0 )c1 /cy
for all t ∈ (t0 , t1 ) .
(3.11)
It follows the boundedness of the differential part y of a solution of (2.11a) and, using (3.8), of its algebraic part z. Next, we prove that hypothesis (3.1) is fulfilled by the current λ defined in (2.20c). Lemma 2. For a fixed t > 0, let (φn , φp , V )(x, t) be a solution of the drift-diffusion equation (2.20a) with boundary condition (2.20b). Then, for the current λT (t) = (−[qdn ∂x φn + qdp ∂x φp ]|(x=0,t) , [qdn ∂x φn + qdp ∂x φp ]|(x=l,t) ) , we indeed have T (AT λ u) λ ≥ 0 .
Proof. We multiply the first two equations in (2.20a) by φn and φp , respectively. Summing up and integrating by parts, we obtain Z l l [φn dn ∂x φn + φp dp ∂x φp ]0 − (dn (∂x φn )2 + dp (∂x φp )2 ) dx 0
=
Z
l
F (φn , φp , V ) · exp 0
φp − φ n UT
− 1 · (φp − φn ) dx ≥ 0 ,
since F ≥ 0 and (exp(x)−1)·x ≥ 0. Moreover, by definition, the functions dn = µn n and dp = µp p are strictly positive for any electric potential V and quasi-Fermi potentials φn , φp . Thus, the integral term on the left-hand side is negative or zero. The claim holds immediately by investing the boundary conditions (2.20b) and the coupling conditions (2.20c).
July 4, 2003 8:40 WSPC/103-M3AS
12
00290
G. Al`ı et al.
3.2. Main theorem We consider the Banach space X = C([t0 , t1 ], L2 ([0, l])) × C([t0 , t1 ], L2 ([0, l])) × C([t0 , t1 ], Rny ) . By taking into account the lemmas of the previous section, we define M ⊂ X as the set of all functions (φn , φp , y) ∈ X which satisfy for all x ∈ [0, l], t ∈ [t0 , t1 ] the inequalities |y|2 (t) ≤ Cy ec1 (t−t0 ) (|y0 |2 + kık2(L2 ([t0 ,t1 ]))nI + kvk2(L2 ([t0 ,t1 ]))nV ) ,
(3.12a)
u1 (t) ∧ u2 (t) ≤ φn (x, t) , φp (x, t) ≤ u1 (t) ∨ u2 (t) ,
almost everywhere ,
and fulfill the following initial value and coupling conditions: ! u1 y(t0 ) = y0 and = AT λ y1 , u2
(3.12b)
(3.12c)
with y = (y1 , y2 )T , y2 (t) ∈ RnL . In (3.12b), the symbols ∧ and ∨ connecting two functions denote the minimum and the maximum, respectively, of the two functions. It is simple to verify that M is a nonempty, bounded, closed, convex subset of X. Theorem 1. (Existence) Let the source functions ı and v be continuous, the network matrices be symmetric, positive definite and the topological conditions (2.16)–(2.18) be fulfilled. Then, problem (2.20) admits a solution contained in M. Proof. The proof is based on a non-standard application of Schauder’s theorem. To define a fixed point map, we fix (φ0n , φ0p , y0 ) ∈ M , with y0 (t) = (y10 (t), y20 (t))T ,
y20 (t) ∈ RnL ,
and proceed in two steps. Firstly, we solve the semilinear, elliptic equation stated as Problem A (Box 2). In (3.13a), the functions n(V, φn ) and p(V, φp ) are defined by (2.5). We denote by V 0 the solution to Problem A. Problem A
Box 2 ε 2 ∂ V = n(V, φ0n ) − p(V, φ0p ) − N , q x BC
V (0, ·) = Vbi (0) + u01 (·), (u01 , u02 ) = y10T Aλ .
V (l, ·) = Vbi (l) + u02 (·) ,
(3.13a) (3.13b) (3.13c)
Secondly, we solve the coupled system of two semilinear elliptic equations and the electric network, which is referred to as Problem B (Box 3). In Problem B,
July 4, 2003 8:40 WSPC/103-M3AS
00290
Elliptic Partial Differential-Algebraic Multiphysics Models
13
Problem B
Box 3
φp − φ n −1 , UT φp − φ n ∂x (d0p ∂x φp ) = F 0 · exp −1 , UT
∂x (d0n ∂x φn ) = −F 0 · exp
BC
φn (0, ·) = φp (0, ·) = u1 (·) ,
AP
(3.14b)
φn (l, ·) = φp (l, ·) = u2 (·) , (u1 , u2 ) = y1T Aλ ,
DAE
(3.14a)
(3.14c)
dy ˜ = 0, + BP y + CP z + FP (ı(t), λ) dt
(3.14d)
BQ z + CQ y + FQ (ı(t), v(t)) = 0 ,
(3.14e)
y(t0 ) = y0 (t0 ) ! −qd0n ∂x φn − qd0p ∂x φp . qd0n ∂x φn + qd0p ∂x φp IV
coupling current
˜= λ
(3.14f) (3.14g)
˜ (3.14g) is primes for d and f denote evaluation at (V 0 , φ0n , φ0p ). The current λ well-defined, since (3.14a) and (3.14b) imply d0n ∂x φn + d0p ∂x φp = constant in [0, l] .
(3.15)
We denote by (φ00n , φ00p , y00 ) the solution to Problem B, with the initial data satisfying y00 (t0 ) = y0 (t0 ) = y0 . Next, we consider the map T :M ⊂X →X,
(φ0n , φ0p , y0 ) 7→ (φ00n , φ00p , y00 ) .
(3.16)
In the next section, we show that both Problems A and B are uniquely solvable (Lemmas 2 and 4). Then, T is well-defined and T (M ) ⊂ M (Lemma 4). Moreover, T is continuous due to the well-posedness of uniformly elliptic problems. To apply Schauder’s fixed point theorem, we only need to show that T (M ) is precompact in M ⊂ X.13 The set M1 = {y00 ∈ C([t0 , t1 ], Rny ) | (φ00n , φ00p , y00 ) ∈ T (M )} a– is equicontinuous, and therefore it is precompact in C([t0 , t1 ], Rny ) by the Arzel` Ascoli theorem. In a similar way, we can show that the set M2 = {(φ00n , φ00p ) ∈ C([t0 , t1 ], (L2 ([0, l]))2 ) | (φ00n , φ00p , y00 ) ∈ T (M )} is precompact, since Lipschitz continuity of (φ00n , φ00p ) with respect to the boundary (Lemma 4) implies k(φ00n , φ00p )(·, t) − (φ00n , φ00p )(·, s)k(L2 )2 ≤ c3 ky00 (t) − y00 (s)kRny ,
July 4, 2003 8:40 WSPC/103-M3AS
14
00290
G. Al`ı et al.
for any t, s ∈ [t0 , t1 ], with c3 independent of the choice of (φ00n , φ00p ). Finally, Tichonov’s theorem implies that T (M ) = M2 × M1 is precompact, and Schauder’s fixed point theorem applies. Next, we discuss a consequence of Theorem 1 concerning the perturbation index 2 of the PDAE system (2.20). Corollary 1. The mixed PDAE system (2.20) equipped with the hypothesis of Theorem 1 has perturbation index 1 with respect to perturbations of the DAE system. Proof. The mixed problem (2.20) can be rewritten as AP
dy + BP y + CP z + FP (ı, λ(AT λ y1 )) = 0 , dt BQ z + CQ y + FQ (ı, v) = 0
(3.17a) (3.17b)
with consistent initial values y(t0 ) = y0 . In (3.17a), the current λ at the interconnects of the diode can be regarded as an operator which maps the boundary voltages AT λ y1 onto the coupling branch current (2.20c) by solving the elliptic BVP (2.20a), (2.20b). Thus we can regard λ as a voltage-controlled current source, and, together with the assumed network topology, this implies12 that the perturbation index of (3.17) is 1. 4. On the Definition of the Fixed-Point Map To show that the fixed-point map T (3.16) is well-defined, we consider separately Problems A and B. 4.1. Existence and uniqueness for Problem A Problem A has the following structure: We wish to determine a function w(x, t) such that ∂x (a(x, t)∂x w) = f (x, t, w) , w(0, t) = w1 (t) ,
for (x, t) ∈ (0, l) × [t0 , t1 ] ,
w(l, t) = w2 (t) ,
for t ∈ [t0 , t1 ] ,
(4.1a) (4.1b)
where a(x, t) and f (x, t, w) satisfy, for all (x, t) ∈ (0, l) × [t0 , t1 ], the conditions: a(x, t) ≥ a > 0 , ∂f (x, t, w) ≥ 0 , for all w ∈ R , ∂w ¯ w) , for all w ∈ R , f(t, w) ≤ f (x, t, w) ≤ f(t,
(4.2) (4.3) (4.4)
¯ there exist w(t), w(t) for some constant a and some functions f, f; ¯ ∈ R, such that f (t, w(t)) ¯ = 0,
¯ w(t)) = 0 . f(t,
(4.5)
July 4, 2003 8:40 WSPC/103-M3AS
00290
Elliptic Partial Differential-Algebraic Multiphysics Models
15
As a consequence of (4.3), we note that the function f (x, t, ·) : R → R is continuous. The assumptions (4.2)–(4.5) are sufficient to prove the following lemma. Lemma 1. Under the hypothesis (4.2)–(4.5), for all t ∈ [t0 , t1 ] there exists a unique solution w(·, t) of problem (4.1) in H 1 ([0, l]). Moreover, this solution satisfies ¯ . min{w1 (t), w2 (t), w(t)} ≤ w(x, t) ≤ max{w1 (t), w2 (t), w(t)}
(4.6)
This result (without the occurrence of the parameter t) has been proved by Markowich7 (see also Taylor’s book11 ). Using standard arguments,7,8 Lemma 1 implies immediately the existence and uniqueness of a solution to Problem A. Lemma 2. Under the assumption φ0p (x, t) ≤ u1 (t) ∨ u2 (t) .
u1 (t) ∧ u2 (t) ≤ φ0n (x, t) ,
(4.7) 0
1
Problem A has for any time t ∈ [t0 , t1 ] a unique solution V = V (·, t) in H ([0, l]). Moreover, the solution satisfies inf Vbi + u1 (t) ∧ u2 (t) ≤ V 0 (x, t) ≤ sup Vbi + u1 (t) ∨ u2 (t) .
[0,l]
(4.8)
[0,l]
4.2. Existence and uniqueness of Problem B First, we investigate the elliptic part of Problem B. This problem has the following structure: We wish to find a pair of functions (wn , wp )(x, t) satisfying, for all (x, t) ∈ (0, l) × [t0 , t1 ], the elliptic system ∂x (an (x, t)∂x wn ) = −f (x, t, wp − wn ) ,
(4.9a)
∂x (ap (x, t)∂x wp ) = f (x, t, wp − wn ) ,
(4.9b)
and the boundary conditions wn (0, t) = wp (0, t) = w1 (t) ,
wn (l, t) = wp (l, t) = w2 (t) .
(4.9c)
Lemma 3. Under the assumptions (4.3)–(4.4), and f (t, 0) = f¯(t, 0) = 0 , an (x, t) ≥ an > 0 ,
ap (x, t) ≥ ap > 0 ,
(4.10) (4.11)
for all (x, t) ∈ (0, l) × [t0 , t1 ], the problem (4.9) has a unique solution (wn , wp ) in C([t0 , t1 ], (H 1 ([0, l]))2 ) for time continuous w1 and w2 with kw1 k∞ , kw2 k∞ ≤ K, K a given constant. This solution satisfies w1 (t) ∧ w2 (t) ≤ wn (x, t) ,
wp (x, t) ≤ w1 (t) ∨ w2 (t) . def
(4.12)
Moreover, the solution (wn , wp ) and the flux I = −an ∂x wn − ap ∂x wp are Lipschitz continuous functions of the data (w1 , w2 ).
July 4, 2003 8:40 WSPC/103-M3AS
16
00290
G. Al`ı et al.
Proof. For each t in the interval [t0 , t1 ], the existence of a solution in (H 1 ([0, l]))2 , which satisfies the estimate (4.12), can be proved by using Schauder’s fixed point theorem in a standard way, with the help of Lemma 1. Next, we prove that a solution (wn , wp ) to (4.9) and the corresponding flux I are Lipschitz continuous with respect to the data (w1 , w2 ). Let (wn0 , wp0 ), (wn00 , wp00 ) be two solutions to the problem (4.9a), (4.9b), corresponding to the data (4.9c), with (w1 , w2 ) equal to (w10 , w20 ), (w100 , w200 ), respectively. We introduce the notation δg(w) = g(w0 ) − g(w00 ) for any suitable functional g, e.g. δwn = wn0 − wn00 , δI = −an ∂x δwn − ap ∂x δwp . Then we have ∂x (an (x, t)∂x δwn ) = −δf (wp − wn ) , ∂x (ap (x, t)∂x δwp ) = δf (wp − wn ) . Multiplying the first and second equation by δwn and δwp , respectively, summing the result and integrating over [0, l] with respect to x, we obtain Z l (an |∂x δwn |2 + ap |∂x δwp |2 ) dx [δwn an ∂x δwn + δwp ap ∂x δwp ]l0 − 0
=
Z
l 0
δf (wp − wn ) (δwp − δwn )2 dx ≥ 0 , δ(wp − wn )
because f is monotone increasing in w. Since δwn = δwp at the boundary and δI is constant with respect to x, the above inequality yields Z l δI(δw2 − δw1 ) + (an |∂x δwn |2 + ap |∂x δwp |2 ) dx ≤ 0 , (4.13) 0
where δw1 = w10 − w100 , δw2 = w20 − w200 . Applying (4.11), there exist two positive constants c1 , c2 , such that δI(δw2 − δw1 ) + c1 |δI|2 ≤ 0 ,
(4.14)
δI(δw2 − δw1 ) + c2 k∂x (δwn , δwp )k2 ≤ 0 .
(4.15)
It follows that |δI| ≤
1 |δw2 − δw1 | c1
(4.16)
and k∂x (δwn , δwp )k2 ≤
1 |δw2 − δw1 |2 , c1 c2
(4.17)
which imply Lipschitz continuity of I and (wn , wp ) with respect to (w1 , w2 ). This result implies that problem (4.9) is uniquely solvable and well-posed. Finally, for time continuous w1 and w2 , the solution (wn , wp ) depends continuously on time. We can immediately apply Lemma 3 to Problem B.
July 4, 2003 8:40 WSPC/103-M3AS
00290
Elliptic Partial Differential-Algebraic Multiphysics Models
17
Lemma 4. Problem B has a unique solution (φ00n , φ00p , y00 ), satisfying u1 (t) ∧ u2 (t) ≤ φ00n (x, t) ,
φ00n (x, t) ≤ u1 (t) ∨ u2 (t) .
(4.18)
˜ are Lipschitz continuous functions of the boundary data Moreover, φ00n , φ00p and λ u1 , u2 , and we have T˜ (AT λ u) λ ≥ 0 .
(4.19)
Proof. Using Lemma 3, the functional I˜ : R2 → R ,
˜ 1 , u2 ) := −qd0 ∂x φu − qd0 ∂x φu , (u1 , u2 ) 7→ I(u n n p p
is well-defined. Here, (φun , φup ) is the solution of (3.14a) and (3.14b) corresponding ˜ = (I, ˜ −I) ˜ T as a Lipschitz continuous to the data (u1 , u2 ). Then we can regard λ function of y. Therefore the differentialnetwork equation (3.14d) is Lipschitz in y and the Picard–Lindel¨ of theorem yields the existence of a unique solution y 00 to the problem (3.14d), (3.14e) corresponding to the initial data (3.14f). Moreover, the triple (φ00n , φ00p , y00 ) solves Problem B — using the notation, for fixed t, φ00 (·, t) = 00 φu (·, t), with (u001 , u002 )T = AT λ y (t). The condition (4.19) follows as in Lemma 2. 5. Conclusions and Outlook In this paper, we have investigated one example of a PDAE, which arises in network analysis: the coupled system of elliptic boundary value problems and differential algebraic network equations, which describe one-dimensional semiconductor devices (diodes) plugged into a discrete modeled electric circuit. In the general case of index-1 networks, we have shown the existence of solutions to the resulting PDAE model (2.20), which can be generalized to nonlinear RLC networks. This analysis presents the first result of this kind for a nonlinear PDE and DAE. We note that the numerical simulation of (2.20)can be realized somewhat straightforwardly: numerical time integration is used for the DAE initial value problem, and an elliptic boundary value problem solver provides the evaluation ˜ at the node voltages required by the time integration scheme. of the operator λ The result on the perturbation index (Corollary 1) enables this straight arithmetic, while the behavior to perturbations of the PDE is still open. The numerics will be addressed in an upcoming paper. And there are a few other open problems: first of all, we can only expect uniqueness of a solution close to equilibrium. Generally, we are faced with bifurcations. Further questions are related to higher-index network equations where some or even all of the topological conditions (2.16)–(2.18) are violated. Here we can expect that the PDAE problem inherits the weakly ill-posedness of ill-conditioned higher-index DAE systems. Acknowledgment The first author wishes to thank Prof. Peter Rentrop for the opportunity to spend two months at IWRMM, where this paper has been conceived. Also, the first author
July 4, 2003 8:40 WSPC/103-M3AS
18
00290
G. Al`ı et al.
wishes to thank Prof. Angelo Marcello Anile for many fruitful discussions and for his encouragement. The subsequent authors express their gratitude towards the BMBF (German Ministry for Education and Research) for supporting this work within the project consortium “Numerische Simulation elektrischer Schaltungen im Zeitbereich” (Karlsruhe/Berlin). References 1. C. W. Gear, Simultaneous numerical solution of differential-algebraic equations, IEEE Trans. Circuit Theory 18 (1971) 89–95. 2. M. G¨ unther, Partielle Differential-Algebraische Systeme in der Numerischen Zeitbereichsanalyse Elektrischer Schaltungen (VDI, 2001). 3. M. G¨ unther, A PDAE model for interconnected linear RLC networks, Math. Comput. Model. Dynam. Systems 7 (2001) 189–203. 4. C. W. Ho, A. E. Ruehli and P. A.Brennan, The modified nodal approach to network analysis, IEEE Trans. Circuits Systems 22 (1975) 505–509. 5. J. W. Jerome, Analysis of Charge Transport. A Mathematical Study of Semiconductor Devices (Springer, 1995). 6. E. Hairer and G. Wanner, Solving Ordinary Differential Equation: II, Stiff and Differential-Algebraic Problems (Springer, 2002), 2nd revised. 7. P. A. Markowich, The Stationary Semiconductor Device Equations (Springer, 1986). 8. P. A. Markowich, C. A. Ringhofer and C. Schmeiser, Semiconductor Equations (Springer, 1990). 9. W. J. McCalla, Fundamentals of Computer Aided Circuit Simulation (Kluwer, 1988). 10. S. Selberherr, Analysis and Simulation of Semiconductor Devices (Springer, 1984). 11. M. E. Taylor, Partial Differential Equations III — Nonlinear Equations (Springer, 1996). 12. C. Tischendorf, Topological index calculation of differential-algebraic equations in circuit simulation, Surv. Math. Ind. 8 (1999) 187–199. 13. E. Zeidler, Nonlinear Functional Analysis and its Applications I : Fixed-Point Theorems (Springer, 1986).
