Skip to content

Chapter 2 — Companion models: discretizing R, L, and C

In Chapter 1 we saw the strategy: discretize time, then turn the network into a resistive circuit the computer can solve once per step. This chapter does the first half rigorously. We derive the companion model (also called the discrete or Norton-equivalent model) of the three lumped elements — resistor, inductor, capacitor — using the trapezoidal rule. By the end you will be able to compute, by hand, the very numbers the simulator produces.

Learning objectives

  • State why inductors and capacitors cannot be handled by ordinary nodal analysis directly.
  • Apply the trapezoidal rule to derive the companion conductance and history-current source of L and C.
  • Compute the companion parameters for given element and time-step values.
  • Hand-iterate the first steps of an R-L or R-C circuit and match the simulator's output.
  • Explain, qualitatively, why the trapezoidal rule is the EMT workhorse.

2.1 The obstacle: two of our elements are differential

Nodal analysis (which we build in Chapter 3) loves elements whose current is an instantaneous function of voltage — i.e. conductances. The resistor obliges:

iR(t)=1RvR(t)=GvR(t).

The inductor and capacitor do not. Their defining laws involve a derivative:

vL(t)=LdiLdt,iC(t)=CdvCdt.

The current now depends on the rate of change, not the present value alone. We cannot stamp a derivative into a linear algebraic system. The fix is to replace the derivative with a difference using a numerical integration rule, evaluated over one time step of length Δt.

2.2 The trapezoidal rule

Suppose a quantity x(t) obeys x˙=f(t), and we know x at the previous step tΔt. The exact update is

x(t)=x(tΔt)+tΔttf(τ)dτ.

The trapezoidal rule approximates the integral by the area of the trapezoid under f between the two endpoints:

tΔttf(τ)dτΔt2(f(t)+f(tΔt)).

It is implicit (the unknown endpoint f(t) appears on the right), second-order accurate (the error per step shrinks like (Δt)3, and globally like (Δt)2), and A-stable (it will not blow up for a stable system, no matter how large Δt is). Those three properties are exactly what an EMT integrator needs. We will apply it to each reactive element.

2.3 Companion model of the inductor

Start from vL=LdiL/dt, i.e. diL/dt=vL/L. Integrate over one step and apply the trapezoidal rule:

iL(t)=iL(tΔt)+Δt2L(vL(t)+vL(tΔt)).

Group the terms into "depends on the present unknown vL(t)" and "known from the previous step":

iL(t)=GLvL(t)+ILhist,GL=Δt2L,ILhist=iL(tΔt)+GLvL(tΔt).

Read this as a circuit: the inductor, over one step, behaves like a conductance GL=Δt/2L in parallel with a current source ILhist whose value is fixed by the previous step's current and voltage. That parallel combination is the inductor's companion model. Once ILhist is computed (it is just a number at solve time), the inductor is indistinguishable from a resistor-plus-source — and those, nodal analysis handles happily.

Companion model of the inductor: over one step it behaves as a conductance GL=Δt/2L in parallel with the history current source ILhist=iL(tΔt)+GLvL(tΔt).

Note the physics baked into GL=Δt/2L: a large inductance gives a tiny conductance (the inductor barely lets current change in one step), which is why an inductor looks almost like an open circuit at the first instant.

2.4 Companion model of the capacitor

Start from iC=CdvC/dt, i.e. dvC/dt=iC/C. Trapezoidal integration of vC gives

vC(t)=vC(tΔt)+Δt2C(iC(t)+iC(tΔt)).

Solve for the present current iC(t):

iC(t)=GCvC(t)+IChist,GC=2CΔt,IChist=(GCvC(tΔt)+iC(tΔt)).

Again a conductance GC=2C/Δt in parallel with a known history source. Here a large capacitance gives a large conductance (the capacitor easily passes current in one step) — the dual of the inductor, exactly as intuition demands: a capacitor looks almost like a short circuit to a sudden change.

Companion model of the capacitor: a conductance GC=2C/Δt in parallel with the history source IChist=(GCvC(tΔt)+iC(tΔt)). It is the exact dual of the inductor.

2.5 The resistor (for completeness)

No derivative, no history term:

iR(t)=GvR(t),G=1R.

2.6 The three companion models at a glance

ElementConductance GeqHistory source Ihist (known from tΔt)
Resistor R1/R— (none)
Inductor LΔt/2Li(tΔt)+GLv(tΔt)
Capacitor C2C/Δt(GCv(tΔt)+i(tΔt))

The history source is the memory of the element: it carries stored energy forward from one step to the next. Because it depends only on already-known quantities, it is a constant by the time the present step's network is solved.

2.7 Worked example: the first steps of the R-L step response

Let us reproduce, by hand, what the simulator does on the R-L step response sample (V=10 V DC, R=10 Ω, L=0.1 H, Δt=50 μs). The inductor's lower terminal is grounded, so the inductor voltage equals the mid-node voltage, vL=Vmid.

The R-L step-response circuit: an ideal 10 V source, series R=10 Ω, and L=0.1 H whose lower terminal is grounded, so vL=Vmid.

Companion conductance.

GL=Δt2L=50×1062(0.1)=2.5×104 S.

Step 1 (t=Δt). The simulator starts every branch current and the history term at zero, so ILhist=0. Kirchhoff's current law at the mid-node says the current through R equals the inductor current:

VVmidR=GLVmid+ILhist.10Vmid10=2.5×104Vmid+0Vmid=10.10025=9.975 V,iL(Δt)=GLVmid=2.5×104(9.975)=2.494×103 A.

Now open the sample, run it, and read the very first sample on the scope (or in the exported CSV): Vmid=9.975 V, IL=0.002494 A. Our hand calculation matches the simulator to four significant figures — because we performed exactly the arithmetic the engine performs.

Two physical sanity checks fall out immediately:

  • The inductor current is almost zero at the first instant — the inductor resists the sudden change, behaving like a near-open circuit (its companion conductance GL is tiny).
  • Almost the entire source voltage appears across the inductor (Vmid10 V), again because it is momentarily open.

Onward in time. Each step updates ILhist from the previous step and re-solves. Marched forward, the current traces the exponential rise iL(t)=(V/R)(1et/τ) with time constant

τ=LR=0.110=0.01 s=10 ms.

At t=τ the analytic value is iL=(1)(1e1)=0.632 A; the simulator reports 0.631 A — agreement to about a tenth of a percent. The two values are not identical because the engine seeds its history from zero current at t=0 (a one-step initialization choice) and because trapezoidal integration carries its own small, second-order error. Both effects are tiny here and shrink as you reduce Δt.

2.8 Lab: confirm the companion model

Run both first-order samples and verify the theory you just derived.

Open the R-L step response in app →

Open the R-C charging circuit in app →

For the R-L circuit:

  1. Confirm the final current settles near V/R=1 A.
  2. Find the time at which iL reaches 0.632 A — it should be one time constant, ~10 ms.
  3. Read the first scope sample and check it against the hand calc above.
The R-C charging circuit: an ideal 10 V source, series R=1 kΩ, and C=10 μF to ground; the capacitor-node voltage is Vcap.

For the R-C circuit (R=1 kΩ, C=10 μF), repeat the derivation yourself (Problem 2.3) and confirm:

  1. GC=2C/Δt=0.4 S.
  2. The first sample is Vcap=0.0249 V, iC=0.009975 A.
  3. Vcap reaches 63% of 10 V (i.e. 6.32 V) after τ=RC=10 ms, and the charging current starts at V/R=10 mA and decays to zero.

2.9 Why the trapezoidal rule?

EMT engines overwhelmingly use the trapezoidal rule because it is:

  • A-stable — large steps cannot make a stable circuit's simulation diverge (an explicit rule like forward Euler can).
  • Second-order accurate — halving Δt quarters the error, cheap accuracy.
  • Symmetric, so for a fixed Δt the companion conductances GL, GC are constant and the system matrix can be factorized once and reused (Chapter 3).

Its one notorious weakness is a tendency to produce small sustained numerical oscillations right after an abrupt event (for example, opening a switch that carries inductor current). Production engines damp these with techniques such as critical-damping adjustment or a brief switch to backward Euler around discontinuities; we return to this in Module 2 and the appendices.

2.10 Summary

  • Resistors are conductances; inductors and capacitors are differential and need discretizing.
  • The trapezoidal rule turns each reactive element into a conductance in parallel with a known history current source — its companion model: GL=Δt/2L and GC=2C/Δt.
  • The history source encodes the element's stored energy and is a known constant when the present step is solved.
  • Hand-iterating the companion recursion reproduces the simulator's output exactly; over a time constant it tracks the analytic exponential to a fraction of a percent.

2.11 Problems

Problem 2.1. A 5 mH inductor is simulated with Δt=20 μs. Compute its companion conductance GL. Repeat for Δt=50 μs. What happens to GL as you refine the step?

Solution 2.1

GL=Δt/2L. For Δt=20 μs: GL=20×106/(25×103)=2×103 S. For Δt=50 μs: GL=50×106/0.01=5×103 S. GL is proportional to Δt, so a smaller step gives a smaller companion conductance — the inductor looks "more open" over a shorter interval, consistent with v=Ldi/dt.

Problem 2.2. A 100 μF capacitor is simulated with Δt=50 μs. Compute GC. If at step n1 the capacitor had v=4 V and i=0.2 A, compute the history source IChist to be used at step n.

Solution 2.2

GC=2C/Δt=2(100×106)/(50×106)=4 S. IChist=(GCv(tΔt)+i(tΔt))=(44+0.2)=16.2 A.

Problem 2.3. Derive the single-node update equation for the R-C charging circuit (ideal source V, series R, capacitor C to ground) using the capacitor companion model, then compute Vcap and iC at the first step (V=10, R=1000, C=10 μF, Δt=50 μs, zero initial history).

Solution 2.3

KCL at the capacitor node: (VVcap)/R=iC=GCVcap+IChist. With GC=2C/Δt=0.4 S and IChist=0 at step 1:

10Vcap1000=0.4Vcap0.01=(0.001+0.4)VcapVcap=0.010.401=0.02494 V.iC=GCVcap=0.4(0.02494)=0.009975 A.

These match the simulator's first sample exactly. Note iC starts near V/R=10 mA — the capacitor looks momentarily like a short.

Problem 2.4. Trapezoidal integration is second-order accurate. If a study has unacceptable integration error at Δt=40 μs, roughly what step would cut the error by a factor of 16, and what is the cost?

Solution 2.4

Global error scales like (Δt)2. To cut error by 16 you need (Δt)2(Δt)2/16, i.e. ΔtΔt/4=10 μs. The cost is 4× as many time steps (4× the run time), since the number of steps scales like 1/Δt.

Problem 2.5. Explain physically why GC=2C/Δt grows as the step shrinks, while GL=Δt/2L shrinks. Relate each to how the element responds to a fast change.

Solution 2.5

Over a very short interval a capacitor can absorb a large current for a small voltage change (i=Cdv/dt is large when dt is small), so it looks like a large conductance — hence GC1/Δt. An inductor resists current change over a short interval (v=Ldi/dt is large for small dt), so very little current flows for a given voltage: it looks like a small conductance — hence GLΔt. The two are duals.

2.12 References

  • H. W. Dommel, Electromagnetic Transients Program (EMTP) Theory Book, Bonneville Power Administration — the companion-model derivations for lumped elements.
  • J. Arrillaga and N. R. Watson, Power Systems Electromagnetic Transients Simulation, IET Power and Energy Series 39 — numerical integrator substitution and the discretization of R, L, C.

Previous: Chapter 1 — What EMT simulation is · Next: Chapter 3 — Building and solving the network each step.