2.4 Summary of This Chapter
3.1.2 Time Domain Solution Algorithm
In order to simulate the circuit shown in Fig. 2.2, we formulate a system of ODEs and algebraic equations by using the KCL for this circuit (see Chapter 2, equations (2.28)-(2.30) and (2.31)-(2.34)). For a specified time period T , these nonlinear ODEs needs to be solved in the time domain. The following steps are taken to accomplish this [35, 82]:
• Step 1: The initial time step t, total time period T , and time step 4t. are determined first.
• Step 2: The decoupling method is used to decouple all equations (2.28)-(2.30) and (2.31)-(2.34).
• Step 3: Each decoupled ODE is solved sequentially with the MI and Runge-Kutta methods.
• Step 4: Each MI loop is convergence tested.
• Step 5: Overall outer loop is convergence tested.
• Step 6: If the specified stopping criterion is not met for the outer loop, all data is updated with newer results and recalculate fromStep 3.
• Step 7: Step 3 - 6 are repeated until the time step meets the specified time period T .
First, all coupled ODEs are decoupled by the waveform relaxation (WR) method [30]-[32]. The WR method starts with a general nonlinear system of
D ODEs with associated initial conditions, dV
dt = f (V, t), V (0) = V0 , t ∈ [0, T ] , (3.14) where T > 0, f : RD× [0, T ] −→ RD, V0 = [V1,0, V2,0, · · · , VD,0] ∈ RDis the initial vector of V , and V (0) = [V1(t), V2(t), · · · , VD(t)] ∈ RDis the solution vector. The system can be written as follows:
The WR method for solving Eq. (3.14) is a continuous-time iterative method.
Therefore, given a function which approximates the solution, it calculates a new approximation along the whole time-interval of interest. Clearly, it differs from most standard iterative techniques in that its iterations are functions in time instead of a scalar value. The iteration formula is chosen in such a way that one avoids having to solve a large system ODEs. A particularly simple, but often very effective iteration scheme which maps the old iterate V(n−1) is written here:
This is similar to the Gauss-Seidel method for iteratively solving linear and nonlinear systems of algebraic equations; the so-called Gauss-Seidel wave-form relaxation scheme. This scheme converts the task of solving a differen-tial equation in D variables into a task of solving a sequence of differendifferen-tial equations in a single variable. A closely related iteration is the Jacobi wave-form relation scheme, the iteration wave-formula is given by:
Note that the Jacobi algorithm is fully parallel. The equations can be solved simultaneously. The Gauss-Seidel waveform relaxation algorithm is formu-lated in the following.
Algorithms:
n = 0; choose Vx0(t)for t ∈ [0, T ], x = 1, 2, · · · , D do
n = n + 1
for x = 1, 2, · · · , D
solve dtdVx(n)= fx(V1(n), V2(n), · · · , Vx−1(n), Vx(n), Vx+1(n−1), · · · , VD(n−1), t) with Vxn(0) = Vx,0
end for
while |Vxn− Vxn−1| < TOL, for x = 1, 2, · · · , D
Each decoupled ODE is solved with the MI algorithm. To clarify the MI algorithm for the numerical solution of the decoupled nonlinear ODEs, we write the ODEs (2.28)-(2.30) and (2.31) in the following form
dVx
dt = f (Vx, t), (3.18)
where Vx is the unknown to be solved and f is the collection of nonlinear functions from each node. We define the MI parameter λ = ∂V∂fx and insert λ into Eq. (3.18), arriving at the MI equation
dVx
dt = f (η, t) − λ(Vx− η) , (3.19) where v0 ≤ η ≤ ω0 is a value in [0, T ], and v0 and ω0 are the lower and upper solutions of Eq. (3.18), respectively. Based on nonlinear behaviors of each decoupled circuit ODEs,we will show mathematically that the solution algorithm exhibits monotone convergence.
To clarify the MI algorithm [35, 36, 81] for the numerical solution of the decoupled nonlinear ODEs, we rewrite the decoupled ODEs as the following form:
dVX(g)
dt = f (VX(g), t) VX(g)(0) = VX(g)0
, (3.20)
where VX(g) is the unknown to be solved, g is the decoupling index (g = 0, 1, 2, · · · ). We note that the f is a collection of nonlinear functions with f ∈ C[I × R, R] and I = [0, T ]. For a fixed index g and X, since the upper and lower solutions, V(g)X and V(g)X ,exist in the circuit and V(g)X ≥ V(g)X , we can prove the solution existence in the set Ω = {(t, VX(g)) | V(g)X ≥ VX(g) ≥ V(g)X , ∀t ∈ I} for each decoupled circuit ODE.
Theorem 3.1.1 Let V(g)X and V(g)X be the upper and lower solutions of Eq. (3.20) in C1[I × R, R] such that V(g)X ≥ V(g)X in the time interval I and f ∈ C[I × R, R]. Then there exists a solution VX(g) of Eq. (3.20) such that V(g)X ≥ V(g)X in the time interval I.
Proof 3.1.1 It is a direct result with the continuous property of f, here the comparison theorem is applied [35].
Remark 3.1.1 We note that for each decoupled ODE, the nonlinear function f is nonincreasing function of the unknown VX(g) and the upper and lower solutions V(g)X (0)and V(g)X (0)of Eq. (3.20) in I can be found. We can further prove there exists a unique solution VX(g) of Eq. (3.20) in I and V(g)X (0) ≥ VX(g) ≥ V(g)X (0).
We see that theTheorem 3.1.1 provides an existence condition of the problem, and we can now describe a monotone constructive method for the computer simulation of the circuit ODEs. The constructed sequences will converge to the solution of Eq. (3.20) for all decoupled ODEs in the circuit simulation. In this condition, instead of original nonlinear ODE to be solved, a transformed ODE is solved with such as the RK method. Now we state the main result for the solution of each decoupled circuit ODEs.
Theorem 3.1.2 Let the f ∈ C[I × R, R], V(g)X (0) and V(g)X (0) be the upper and lower solutions of Eq. (3.20) in I. Because f(t, VX(g)) − f(t, eVX(g)) ≥
−λ(VX(g) − eVX(g)), V(g)X (0) ≥ VX(g) ≥ eVX(g) ≥ V(g)X (0) and λ ≥ 0, sequences V(g)Xn unif.−→ V(g)X and V(g)Xn unif.−→ V(g)X as n −→ ∞ monotonically in I [35].
Proof 3.1.2 For V ∈ C[I × R, R] such that V(g)X ≥ V ≥ V(g)X , we consider the following transformed ODE equation for the fixed g and X
dVX
dt = f (V, t) − λ(VX(g)− V) VX(g)(0) = VX(g)0
, (3.21)
then ∀V, ∃!VX(g) of Eq. (3.20) in I. Define ΘV = VX(g) and we can prove: (i) ΘV(g)X (0) ≥ V(g)X (0) and ΘV(g)X (0) ≤ V(g)X (0); (ii) Θ is a monotone operator in [V(g)X (0), V(g)X (0)] ≡ [VX(g) ∈ C[I, R] | V(g)X (0) ≥ VX(g) ≥ V(g)X (0)].
Now we construct two sequences by using the mapping Θ : ΘV(g)Xn = V(g)Xn+1 and ΘV(g)Xn = V(g)Xn+1 and by above observations, the following rela-tion holds
V(g)X0 ≥ · · · ≥ V(g)Xn ≥ V(g)Xn ≥ · · · ≥ V(g)X0
in I. Hence V(g)Xn unif.−→ V(g)X and V(g)Xn unif.−→ V(g)X as n −→ ∞ monotonically in
respectively. Thus V(g)X and V(g)X are the solution of Eq. (3.20).
Theorem 3.1.3 For decoupled ODEs. the nonlinear function f is nonincreas-ing in VX(g) and f(t, VX(g)1) − f(t, VX(g)2) ≥ −λ(VX(g)1 − VX(g)2), ∀VX(g)1 ≥ VX(g)2. Thus n
V(g)Xno∞
n=1 and n
V(g)Xno∞
n=1 converge uniformly and monotonically to the unique solution VX(g) of Eq. (3.20) [35].
Proof 3.1.3 By using Theorem 2 and note the nonincreasing property of f, the result follows directly.
Fig. 3.3 shows the flowchart of the proposed time domain simulation method-ology. In order to solve the equations (2.28)-(2.30) and (2.31)-(2.34) in the time domain, we solve VC, VE, VB, VBX, VCX, and VEX sequentially by the decoupled scheme of the system ODEs. Here we state the computational steps of the methodology:
• Step 1: For a given time interval T , we obtain the DC solution as an initial guess and start the time evolution processes.
• Step 2: Giving a time step (t), we can apply the MI and Runge-Kutta (RK) method [20, 21] to solve the decoupled ODE in Eq. (2.28) with the solved solution (V(t−1)C , V(t−1)E , V(t−1)B , V(t−1)BX , V(t−1)CX ). The notation V(t−1)X represents the solution of voltage VX at time step (t − 1). In the monotone loop, V(t−1)C has been solved with the RK formula (fourth-order RK method is used in practical case). The estimation time step of RK method is adaptive to fit the error tolerance cri-terion in this work [20].
• Step 3: If the norm error |V(t)C − V(t)CO| is less than the given error tolerance value (TOL), the solution algorithm will exit the monotone loop and enter the next decoupled ODE Eq. (2.29).
Otherwise we update V(t)COwith V(t)C and continue the mono-tone loop described inStep 2 until the tolerance condition has been achieved.
• Step 4: With the newly solved V(t)C , we can solve the ODE Eq. (2.29) with the same method as described inStep 2.
• Step 5: When the norm error |V(t)E −V(t)EO| is less than TOL, we solve the next ODE Eq. (2.30) with the new results (V(t)C ,V(t)E ).
Otherwise we return to theStep 4.
• Step 6: Similar iterative loops, such as that formed by Step 2 and Step 3, are executed for all unknown variables. As (V(t)C , V(t)E , V(t)B , V(t)BX, V(t)CX) converge in their respective inner loops, we apply the convergence test for the overall outer loop. If the outer loop error is less than the TOL, this calculated set of (V(t)C , V(t)E , V(t)B , V(t)BX, V(t)CX) is the solution for this given time step (t). The process is then repeated for the next time step (t + 1). If the TOL is not satisfied, we update all (V(t)CO, V(t)EO, V(t)BO, V(t)BXO, V(t)CXO) and continue the inner loop for each ODE until the outer loop tolerance condition is met.
• Step 7: When t exceeds the given time interval T , the simulation will terminate and the results be ready for further processes, such as fast Fourier transform (FFT).
The convergence properties of this time-domain solution algorithm are presented in Appendix C.
For a time step t < T, let VX(t-1)and VXO(t), be given where X = C, E, B, BX, CX, EX , etc.
Solve the corresponding equation with the RK formula VC(t) = VC(t-1)+ hF1(h, VCO(t), VEO(t), VBO(t), ..., VB(t-1), ....)
Solve the corresponding equation with the RK formula VE(t)= VE(t-1) + hF2(h, VC(t), VEO(t), VBO(t), ..., VB(t-1), ....)
Collect the computed solutions at the present time step t, and set t t + t Solve the Eq. (2.28) Eq.(2.34) in the DC condition
Let a finite time interval be given and start the time evolution processes with the DC solution
Apply the MI method to Eq. (2.28) for unknown
Apply the MI method to Eq. (2.29) for unknown
Obtain all computed data for post processes T
VC
VE
Figure 3.2: A flowchart of the proposed time-domain simulation methodology.