Harmonic Balance

Harmonic Balance Analysis¤

This notebook demonstrates the Harmonic Balance (HB) solver on two circuits:

Series LCR resonator — a linear benchmark. HB recovers the resonance peak and is validated against the analytical transfer function.
Diode half-wave clipper — a nonlinear circuit. The clipped waveform generates harmonics at $2f_0$, $3f_0$, … that are captured by HB.

Unlike transient simulation, HB finds the periodic steady state directly in ~10–20 Newton steps.

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np
import schemdraw
import schemdraw.elements as elm

from circulax import compile_circuit
from circulax.components.electronic import Capacitor, Diode, Inductor, Resistor, VoltageSourceAC

jax.config.update("jax_enable_x64", True)

Part 1: Series LCR Resonator¤

A sinusoidal voltage source drives a series $R$–$L$–$C$ circuit. We measure $V_\text{out}$ across the capacitor.

\[H(j\omega) = \frac{V_\text{out}}{V_\text{in}} = \frac{1/j\omega C}{R + j\omega L + 1/j\omega C} = \frac{1}{1 - \omega^2 LC + j\omega RC}\]

The resonant frequency is $f_0 = 1 / (2\pi\sqrt{LC})$ and the $Q$-factor is $Q = \sqrt{L/C}/R$. At resonance the capacitor voltage is amplified by $Q$ relative to the drive.

# Circuit parameters
R_val = 10.0  # Ω
L_val = 1e-6  # H  (1 µH)
C_val = 1e-9  # F  (1 nF)
V_amp = 1.0  # V (peak)

f_res = 1.0 / (2 * np.pi * np.sqrt(L_val * C_val))  # ≈ 5.033 MHz
Q = np.sqrt(L_val / C_val) / R_val  # ≈ 3.16
f_drive = f_res  # drive at resonance for maximum response

print(f"Resonant frequency : {f_res / 1e6:.3f} MHz")
print(f"Q-factor           : {Q:.3f}")
print(f"|H(jω₀)|           : {Q:.3f}  (capacitor voltage gain at resonance)")

Resonant frequency : 5.033 MHz
Q-factor           : 3.162
|H(jω₀)|           : 3.162  (capacitor voltage gain at resonance)

with plt.style.context(["default", {"axes.grid": True, "figure.facecolor": "white"}]), schemdraw.Drawing() as d:
    d.config(fontsize=13)
    Vs = d.add(elm.SourceSin().up().label("$V_s$\n1 V", loc="left"))
    d.add(elm.Line().right(d.unit * 0.5))
    d.add(elm.Resistor().right().label(f"$R$\n{R_val:.0f} Ω"))
    d.add(elm.Inductor2().right().label(f"$L$\n{L_val * 1e6:.0f} µH"))
    d.add(elm.Line().right(d.unit * 0.5))
    top_right = d.here
    d.add(elm.Capacitor().down().label(f"$C$\n{C_val * 1e9:.0f} nF", loc="right").label("$V_\\mathrm{out}$", loc="top"))
    d.add(elm.Line().left().tox(Vs.start))
    d.add(elm.Ground())
    d.add(elm.Line().right().tox(Vs.start))
    d.add(elm.Dot(open=True).at(top_right).label("$V_\\mathrm{out}$", loc="right"))

# Netlist: Vs -> R -> L -> C -> GND
lcr_net = {
    "instances": {
        "Vs": {"component": "vsrc", "settings": {"V": V_amp, "freq": f_drive}},
        "R1": {"component": "resistor", "settings": {"R": R_val}},
        "L1": {"component": "inductor", "settings": {"L": L_val}},
        "C1": {"component": "capacitor", "settings": {"C": C_val}},
        "GND": {"component": "ground"},
    },
    "connections": {
        "Vs,p1": "R1,p1",
        "R1,p2": "L1,p1",
        "L1,p2": "C1,p1",
        "C1,p2": "GND,p1",
        "Vs,p2": "GND,p1",
    },
    "ports": {"in": "Vs,p1", "out": "C1,p1"},
}

models = {
    "vsrc": VoltageSourceAC,
    "resistor": Resistor,
    "inductor": Inductor,
    "capacitor": Capacitor,
    "ground": lambda: 0,
}

circuit = compile_circuit(lcr_net, models, backend="dense")
num_vars = circuit.sys_size
print(f"System size: {num_vars} variables")
print(f"Named ports: {list(lcr_net['ports'])}")

System size: 6 variables
Named ports: ['in', 'out']

# DC operating point (zero for a purely AC circuit)
y_dc = circuit.dc()

# Harmonic Balance: 5 harmonics -> K = 11 time points per period
N_harmonics = 5
y_time, y_freq = circuit.hb(freq=f_drive, harmonics=N_harmonics, y0=y_dc)

print(f"y_time shape : {y_time.shape}  (K={2 * N_harmonics + 1} time points × {num_vars} variables)")
print(f"y_freq shape : {y_freq.shape}  ({N_harmonics + 1} harmonics × {num_vars} variables)")

y_time shape : (11, 6)  (K=11 time points × 6 variables)
y_freq shape : (6, 6)  (6 harmonics × 6 variables)

K = 2 * N_harmonics + 1
T = 1.0 / f_drive
t_ns = np.linspace(0, T * 1e9, K, endpoint=False)  # nanoseconds

vin_t = circuit.port(y_time, "in")
vout_t = circuit.port(y_time, "out")
vin_f = circuit.port(y_freq, "in")
vout_f = circuit.port(y_freq, "out")

fig, ax = plt.subplots(figsize=(8, 3.5))
ax.plot(t_ns, np.array(vin_t), "C0o-", ms=6, label=r"$V_\mathrm{in}$")
ax.plot(t_ns, np.array(vout_t), "C1s-", ms=6, label=r"$V_\mathrm{out}$")
ax.set_xlabel("Time (ns)")
ax.set_ylabel("Voltage (V)")
ax.set_title("LCR — time-domain waveform (one period)")
ax.legend()
ax.grid(True, alpha=0.4)
plt.tight_layout()
plt.show()

# Two-sided amplitude: multiply by 2 for k>=1 (rfft folds negative frequencies)
harmonics = np.arange(N_harmonics + 1)
scale = np.where(harmonics == 0, 1.0, 2.0)
vin_amp = scale * np.abs(np.array(vin_f))
vout_amp = scale * np.abs(np.array(vout_f))
print(f"|V_in  @ f_drive| = {vin_amp[1]:.4f} V  (expected {V_amp:.4f} V)")
print(f"|V_out @ f_drive| = {vout_amp[1]:.4f} V  (expected Q={Q:.4f} V at resonance)")

|V_in  @ f_drive| = 1.0000 V  (expected 1.0000 V)
|V_out @ f_drive| = 3.1623 V  (expected Q=3.1623 V at resonance)

Validation against the analytical transfer function¤

We sweep frequency and compare $|H(j\omega)|$ from the analytical formula with the HB result (a single point at $f_\text{drive}$).

freqs = np.logspace(5, 8, 500)  # 100 kHz → 100 MHz
w = 2 * np.pi * freqs
H = 1.0 / (1 - w**2 * L_val * C_val + 1j * w * R_val * C_val)

fig, ax = plt.subplots(figsize=(8, 3.5))
ax.semilogx(freqs / 1e6, 20 * np.log10(np.abs(H)), "C0", lw=2, label="Analytical |H(jω)|")
ax.scatter([f_drive / 1e6], [20 * np.log10(vout_amp[1] / vin_amp[1])], color="C1", zorder=5, s=80, label="HB result at $f_{drive}$")
ax.axvline(f_res / 1e6, color="gray", ls="--", lw=1, label=f"$f_0$ = {f_res / 1e6:.2f} MHz")
ax.set_xlabel("Frequency (MHz)")
ax.set_ylabel("Magnitude (dB)")
ax.set_title("Bode magnitude — HB vs analytical")
ax.legend()
ax.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.show()

Frequency Sweep with `jax.vmap`¤

circuit.hb(...) exposes the HB solve through the same high-level parameter-update API used by DC, transient, and AC. The solve is still vmappable over the drive frequency, and the source model's own freq setting is updated at the same time.

The pattern is:

def hb_solve_freq(freq):
    _, y_freq = circuit.hb(
        freq=freq,
        harmonics=N_harmonics,
        y0=y_dc,
        params={"Vs.freq": freq},
    )
    return circuit.port(y_freq, "out")

y_freq_sweep = jax.jit(jax.vmap(hb_solve_freq))(sweep_freqs)

jax.vmap batches freq across the call. jax.jit then compiles the vectorised computation as a single XLA program.

# A thin wrapper that exposes freq as a JAX argument, making the function vmappable.
def hb_solve_freq(freq):
    _, y_freq = circuit.hb(
        freq=freq,
        harmonics=N_harmonics,
        y0=y_dc,
        params={"Vs.freq": freq},
    )
    return y_freq  # shape: (N_harmonics+1, num_vars)


# Five log-spaced frequencies spanning the LCR passband
sweep_freqs = jnp.geomspace(f_res * 0.35, f_res * 5.0, 10)
print("Sweep frequencies:", [f"{float(f) / 1e6:.3f} MHz" for f in sweep_freqs])

# jax.vmap maps the HB solve over all 10 frequencies in a single XLA compilation.
y_freq_sweep = jax.jit(jax.vmap(hb_solve_freq))(sweep_freqs)
# y_freq_sweep: shape (10, N_harmonics+1, num_vars)

# |H(jw)| = |V_out| / |V_in| at the fundamental (k=1)
vin_sweep = circuit.port(y_freq_sweep, "in")
vout_sweep = circuit.port(y_freq_sweep, "out")
H_hb_sweep = jnp.abs(vout_sweep[:, 1]) / jnp.abs(vin_sweep[:, 1])

# --- Plot against the analytical Bode curve ---
fig, ax = plt.subplots(figsize=(8, 3.5))
ax.semilogx(freqs / 1e6, 20 * np.log10(np.abs(H)), "C0", lw=2, label="Analytical |H(jω)|")
ax.scatter(
    np.array(sweep_freqs) / 1e6,
    20 * np.log10(np.array(H_hb_sweep)),
    color="C2",
    zorder=5,
    s=100,
    marker="D",
    label=f"HB vmap ({len(sweep_freqs)} freqs)",
)
ax.axvline(f_res / 1e6, color="gray", ls="--", lw=1, label=f"$f_0$ = {f_res / 1e6:.2f} MHz")
ax.set_xlabel("Frequency (MHz)")
ax.set_ylabel("Magnitude (dB)")
ax.set_title("Bode magnitude — vmapped HB sweep vs analytical")
ax.legend()
ax.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.show()

print("\n|H(jω)| at sweep points:")
for f, Hv in zip(sweep_freqs, H_hb_sweep):
    w_rad = 2 * np.pi * float(f)
    H_exact = abs(1.0 / (1 - w_rad**2 * L_val * C_val + 1j * w_rad * R_val * C_val))
    print(f"  {float(f) / 1e6:.3f} MHz:  HB = {float(Hv):.4f},  analytical = {H_exact:.4f}")

Sweep frequencies: ['1.762 MHz', '2.367 MHz', '3.181 MHz', '4.274 MHz', '5.744 MHz', '7.718 MHz', '10.371 MHz', '13.936 MHz', '18.727 MHz', '25.165 MHz']

|H(jω)| at sweep points:
  1.762 MHz:  HB = 1.1306,  analytical = 1.1306
  2.367 MHz:  HB = 1.2612,  analytical = 1.2612
  3.181 MHz:  HB = 1.5799,  analytical = 1.5799
  4.274 MHz:  HB = 2.5834,  analytical = 2.5834
  5.744 MHz:  HB = 2.1242,  analytical = 2.1242
  7.718 MHz:  HB = 0.6964,  analytical = 0.6964
  10.371 MHz:  HB = 0.3020,  analytical = 0.3020
  13.936 MHz:  HB = 0.1487,  analytical = 0.1487
  18.727 MHz:  HB = 0.0775,  analytical = 0.0775
  25.165 MHz:  HB = 0.0416,  analytical = 0.0416

Part 2: F-domain equivalents of Capacitor and Inductor¤

Time-domain components encode reactive behaviour through the $Q$-term in the DAE:

Component	Time-domain	Frequency-domain admittance
Capacitor	$q = C \cdot V$, so $i = C\,dV/dt$	$Y_C(f) = j2\pi f C$
Inductor	state $i_L$, flux $\phi = -L\,i_L$, $v = L\,di/dt$	$Y_L(f) = \dfrac{1}{j2\pi f L}$

Both descriptions are mathematically equivalent in the frequency domain. To show this we repeat the Part 1 LCR analysis after replacing Capacitor and Inductor with @fdomain_component versions and verify the $V_\text{out}$ waveform and harmonic spectrum are numerically identical.

The key point for HB: the time-domain capacitor contributes $$j k\omega_0 \cdot C(V_{1,k} - V_{2,k})$$ to $R_k$ via the $jk\omega Q_k$ path, while the f-domain capacitor contributes $$Y_C(k f_0)(V_{1,k}-V_{2,k}) = jk\omega_0 C\,(V_{1,k}-V_{2,k})$$ directly — these are identical. The same identity holds for the inductor.

from circulax import fdomain_component


@fdomain_component(ports=("p1", "p2"))
def CapacitorFD(f: float, C: float = 1e-12):
    """Capacitor via frequency-domain admittance: Y_C(f) = j·2πf·C.

    - At DC (f=0): Y=0, so the capacitor is an open circuit.
    - At frequency f: Y = jωC, matching the time-domain reactive term C·dV/dt.
    """
    omega = 2.0 * jnp.pi * f
    Y = 1j * omega * C
    return jnp.array([[Y, -Y], [-Y, Y]], dtype=jnp.complex128)


@fdomain_component(ports=("p1", "p2"))
def InductorFD(f: float, L: float = 1e-9):
    """Inductor via frequency-domain admittance: Y_L(f) = 1/(j·2πf·L).

    A 1 nΩ series resistance regularises the DC singularity — the inductor
    remains a near-short at DC while being exact at all nonzero harmonics.
    """
    omega = 2.0 * jnp.pi * f
    Z = 1j * omega * L + 1e-9  # 1 nΩ prevents 1/0 at f=0
    Y = 1.0 / Z
    return jnp.array([[Y, -Y], [-Y, Y]], dtype=jnp.complex128)


print(f"CapacitorFD._is_fdomain = {CapacitorFD._is_fdomain}")
print(f"InductorFD._is_fdomain  = {InductorFD._is_fdomain}")

# Spot-check admittances at the resonant frequency
f_check = f_res
Y_C = 1j * 2 * np.pi * f_check * C_val
Y_L = 1.0 / (1j * 2 * np.pi * f_check * L_val)
print(f"\nAt f_res = {f_res / 1e6:.3f} MHz:")
print(f"  Y_C = {Y_C:.4e} S  (Im > 0 → capacitive)")
print(f"  Y_L = {Y_L:.4e} S  (Im < 0 → inductive)")

CapacitorFD._is_fdomain = True
InductorFD._is_fdomain  = True

At f_res = 5.033 MHz:
  Y_C = 0.0000e+00+3.1623e-02j S  (Im > 0 → capacitive)
  Y_L = 0.0000e+00-3.1623e-02j S  (Im < 0 → inductive)

# Identical topology to the time-domain LCR — only the C and L models change.
lcr_fd_net = {
    "instances": {
        "Vs": {"component": "vsrc", "settings": {"V": V_amp, "freq": f_drive}},
        "R1": {"component": "resistor", "settings": {"R": R_val}},
        "L1": {"component": "inductor_fd", "settings": {"L": L_val}},
        "C1": {"component": "capacitor_fd", "settings": {"C": C_val}},
        "GND": {"component": "ground"},
    },
    "connections": {
        "Vs,p1": "R1,p1",
        "R1,p2": "L1,p1",
        "L1,p2": "C1,p1",
        "C1,p2": "GND,p1",
        "Vs,p2": "GND,p1",
    },
    "ports": {"in": "Vs,p1", "out": "C1,p1"},
}

models_fd = {
    "vsrc": VoltageSourceAC,
    "resistor": Resistor,
    "inductor_fd": InductorFD,
    "capacitor_fd": CapacitorFD,
    "ground": lambda: 0,
}

circuit_fd = compile_circuit(lcr_fd_net, models_fd, backend="dense")
num_vars_fd = circuit_fd.sys_size

print(f"Time-domain LCR  →  system size = {num_vars} variables  (nodes + i_L + i_src)")
print(f"F-domain LCR     →  system size = {num_vars_fd} variables  (nodes + i_src, no i_L state)")
print()

# DC operating point — trivially zero for a purely AC source
y_dc_fd = circuit_fd.dc()
print(f"DC operating point (f-domain): max|y_dc| = {float(jnp.max(jnp.abs(y_dc_fd))):.2e} V")

# Harmonic Balance with same settings as Part 1
y_time_fd, y_freq_fd = circuit_fd.hb(freq=f_drive, harmonics=N_harmonics, y0=y_dc_fd)
print(f"y_time_fd shape: {y_time_fd.shape}")

Time-domain LCR  →  system size = 6 variables  (nodes + i_L + i_src)
F-domain LCR     →  system size = 5 variables  (nodes + i_src, no i_L state)



DC operating point (f-domain): max|y_dc| = 0.00e+00 V


y_time_fd shape: (11, 5)

vin_td = np.array(circuit.port(y_time, "in"))
vout_td = np.array(circuit.port(y_time, "out"))
vin_fd = np.array(circuit_fd.port(y_time_fd, "in"))
vout_fd = np.array(circuit_fd.port(y_time_fd, "out"))

max_err_in = float(jnp.max(jnp.abs(jnp.array(vin_td) - jnp.array(vin_fd))))
max_err_out = float(jnp.max(jnp.abs(jnp.array(vout_td) - jnp.array(vout_fd))))
print(f"Max |ΔV_in|  (time-domain vs f-domain) = {max_err_in:.2e} V")
print(f"Max |ΔV_out| (time-domain vs f-domain) = {max_err_out:.2e} V")

fig, axes = plt.subplots(1, 2, figsize=(12, 3.8))

# ── Left: Vin comparison ──────────────────────────────────────────────────────
ax = axes[0]
ax.plot(t_ns, vin_td, "C0-", lw=2.5, label="Time-domain  (L, C)")
ax.plot(t_ns, vin_fd, "C1--", lw=2, label="F-domain (InductorFD, CapacitorFD)")
ax.set_xlabel("Time (ns)")
ax.set_ylabel("Voltage (V)")
ax.set_title(r"$V_\mathrm{in}$ — time-domain vs f-domain")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.4)

# ── Right: Vout comparison ────────────────────────────────────────────────────
ax = axes[1]
ax.plot(t_ns, vout_td, "C0-", lw=2.5, label="Time-domain  (L, C)")
ax.plot(t_ns, vout_fd, "C1--", lw=2, label="F-domain (InductorFD, CapacitorFD)")
ax.set_xlabel("Time (ns)")
ax.set_ylabel("Voltage (V)")
ax.set_title(r"$V_\mathrm{out}$ — time-domain vs f-domain")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.4)

plt.suptitle(
    f"LCR at resonance — max waveform error: ΔVin={max_err_in:.1e} V,  ΔVout={max_err_out:.1e} V",
    fontsize=10,
    y=1.01,
)
plt.tight_layout()
plt.show()

Max |ΔV_in|  (time-domain vs f-domain) = 6.42e-18 V
Max |ΔV_out| (time-domain vs f-domain) = 3.16e-10 V

Part 3: Diode Half-Wave Clipper¤

A series resistor and diode form a half-wave clipper: the diode blocks the negative half-cycle, so $V_\text{out}$ is a rectified sinusoid. This strongly nonlinear waveform contains harmonics at $2f_0$, $3f_0$, $4f_0$, … — harmonic balance captures all of them simultaneously.

A pure transient simulation would need to run for many cycles before the diode's junction capacitance settles; HB finds the periodic state directly.

with plt.style.context(["default", {"axes.grid": True, "figure.facecolor": "white"}]), schemdraw.Drawing() as d:
    d.config(fontsize=13)
    Vs2 = d.add(elm.SourceSin().up().label("$V_s$\n2 V", loc="left"))
    d.add(elm.Resistor().right().label("$R$\n1 kΩ"))
    top_mid = d.here
    d.add(elm.Diode().right().label("$D_1$"))
    top_right = d.here
    d.add(elm.Resistor().down().label("$R_L$\n10 kΩ", loc="right").label("$V_\\mathrm{out}$", loc="top"))
    d.add(elm.Line().left().tox(Vs2.start))
    d.add(elm.Ground())
    d.add(elm.Dot(open=True).at(top_right).label("$V_\\mathrm{out}$", loc="right"))

f_clip = 1e3  # 1 kHz — low enough that diode junction sees quasi-static operation
V_clip = 2.0  # V (peak) — enough to forward-bias the diode
R_clip = 1e3  # Ω series resistor
RL_clip = 10e3  # Ω load resistor

clipper_net = {
    "instances": {
        "Vs": {"component": "vsrc", "settings": {"V": V_clip, "freq": f_clip}},
        "Rs": {"component": "resistor", "settings": {"R": R_clip}},
        "D1": {"component": "diode", "settings": {}},
        "RL": {"component": "resistor", "settings": {"R": RL_clip}},
        "GND": {"component": "ground"},
    },
    "connections": {
        "Vs,p1": "Rs,p1",
        "Rs,p2": "D1,p1",
        "D1,p2": "RL,p1",
        "RL,p2": "GND,p1",
        "Vs,p2": "GND,p1",
    },
    "ports": {"in": "Vs,p1", "out": "RL,p1"},
}

clip_models = {"vsrc": VoltageSourceAC, "resistor": Resistor, "diode": Diode, "ground": lambda: 0}

circuit_clip = compile_circuit(clipper_net, clip_models, backend="dense")
clip_dc = circuit_clip.dc()

N_clip = 10  # 10 harmonics to capture the rectified waveform
yt_clip, yf_clip = circuit_clip.hb(freq=f_clip, harmonics=N_clip, y0=clip_dc)

print(f"Converged. y_time shape: {yt_clip.shape}")

Converged. y_time shape: (21, 5)

K_clip = 2 * N_clip + 1
T_clip = 1.0 / f_clip
t_ms = np.linspace(0, T_clip * 1e3, K_clip, endpoint=False)

vin_clip_t = circuit_clip.port(yt_clip, "in")
vout_clip_t = circuit_clip.port(yt_clip, "out")
vin_clip_f = circuit_clip.port(yf_clip, "in")
vout_clip_f = circuit_clip.port(yf_clip, "out")

fig, axes = plt.subplots(1, 2, figsize=(12, 3.8))

# Left: time-domain
ax = axes[0]
ax.plot(t_ms, np.array(vin_clip_t), "C0o-", ms=5, label=r"$V_\mathrm{in}$")
ax.plot(t_ms, np.array(vout_clip_t), "C1s-", ms=5, label=r"$V_\mathrm{out}$")
ax.axhline(0, color="gray", lw=0.8, ls="--")
ax.set_xlabel("Time (ms)")
ax.set_ylabel("Voltage (V)")
ax.set_title("Clipper — time-domain waveform")
ax.legend()
ax.grid(True, alpha=0.4)

# Right: frequency spectrum
ax = axes[1]
harms = np.arange(N_clip + 1)
sc = np.where(harms == 0, 1.0, 2.0)
vin_spec = sc * np.abs(np.array(vin_clip_f))
vout_spec = sc * np.abs(np.array(vout_clip_f))

w = 0.35
ax.bar(harms - w / 2, vin_spec, width=w, color="C0", alpha=0.8, label=r"$V_\mathrm{in}$")
ax.bar(harms + w / 2, vout_spec, width=w, color="C1", alpha=0.8, label=r"$V_\mathrm{out}$")
ax.set_xticks(harms)
ax.set_xticklabels([f"{int(k)}$f_0$" if k > 0 else "DC" for k in harms])
ax.set_xlabel("Harmonic")
ax.set_ylabel("Amplitude (V)")
ax.set_title("Clipper — harmonic spectrum")
ax.legend()
ax.grid(True, axis="y", alpha=0.4)

plt.tight_layout()
plt.show()

print("Output harmonic amplitudes:")
for k in range(N_clip + 1):
    label = "DC" if k == 0 else f"{k}f0"
    print(f"  {label:4s}: {vout_spec[k]:.4f} V")

Output harmonic amplitudes:
  DC  : 0.3823 V
  1f0 : 0.6375 V
  2f0 : 0.3454 V
  3f0 : 0.0752 V
  4f0 : 0.0471 V
  5f0 : 0.0368 V
  6f0 : 0.0089 V
  7f0 : 0.0203 V
  8f0 : 0.0029 V
  9f0 : 0.0118 V
  10f0: 0.0070 V

Component	Time-domain	Frequency-domain admittance
Capacitor	\(q = C \cdot V\), so \(i = C\,dV/dt\)	\(Y_C(f) = j2\pi f C\)
Inductor	state \(i_L\), flux \(\phi = -L\,i_L\), \(v = L\,di/dt\)	\(Y_L(f) = \dfrac{1}{j2\pi f L}\)

Harmonic Balance