From 2827c13e974f14805eee59e780bd33c9356f1dde Mon Sep 17 00:00:00 2001
From: William PIAT <william.piat3@gmail.com>
Date: Fri, 16 Oct 2026 18:13:03 +0200
Subject: [PATCH 1/6] add robust init system

---
 src/ElecSolver/TemporalSystemBuilder.py | 141 +++++++++++++++++++++++-
 src/ElecSolver/utils.py                 |  28 ++++-
 zensical.toml                           |   1 +
 3 files changed, 162 insertions(+), 8 deletions(-)

diff --git a/src/ElecSolver/TemporalSystemBuilder.py b/src/ElecSolver/TemporalSystemBuilder.py
index 0353517..f00457a 100644
--- a/src/ElecSolver/TemporalSystemBuilder.py
+++ b/src/ElecSolver/TemporalSystemBuilder.py
@@ -1,10 +1,10 @@
 import numpy as np
 from scipy.sparse import coo_matrix, block_diag, coo_array
-from .utils import SolutionTemporal,GradientsParametersTemporal,compute_graph_components
+from .utils import SolutionTemporal,GradientsParametersTemporal,compute_graph_components,get_roots
 import warnings
 
 class TemporalSystemBuilder():
-    def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,inductive_mutual_coords,inductive_mutual_data,res_mutual_coords,res_mutual_data):
+    def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,inductive_mutual_coords,inductive_mutual_data,res_mutual_coords,res_mutual_data,size=None):
         """Initialize a temporal electrical system builder.
 
         Parameters
@@ -35,6 +35,8 @@ def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_dat
         res_mutual_data : numpy.ndarray, shape (P,)
             Mutual resistance values. Each sign follows the node order in
             ``coil_coords`` or ``res_coords``.
+        size : int, optional
+            The size of the system. If not provided, it will be determined during graph analysis.
         """
 
         self.coil_coords=coil_coords
@@ -58,6 +60,7 @@ def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_dat
         self.rhs = (np.array([]),(np.array([],dtype=int),))
 
         self.analysed=False
+        self.size = size
 
     def graph_analysis(self):
         """Analyze connected components and initialize system dimensions.
@@ -84,7 +87,10 @@ def graph_analysis(self):
         self.all_impedences = np.concatenate([self.coil_data,self.res_data,self.capa_data,self.voltage_source_data],axis=0)
 
         # Actual number of nodes in the system.
-        self.size = np.max(self.all_coords)+1
+        if self.size is None:
+            self.size = np.max(self.all_coords)+1
+        else:
+            self.size = max(np.max(self.all_coords)+1, self.size)
         # Number of currents.
         self.number_intensities = self.all_impedences.shape[0]
         # Keep the connected subgraphs.
@@ -455,6 +461,135 @@ def get_init_system(self,sparse_rhs=True):
             np.add.at(rhs, nodes, data_rhs)
         return sys,rhs
 
+    def get_robust_init_system(self,sparse_rhs=True):
+        """Fonction to return a more robust initial system and the right hand side
+        in case the initial system is singular, this comes at the cost of additionnal
+        computational cost.
+        It merges the capacitors nodes
+        It removes the coils wires, the mutuals for the coils
+        Returns
+        -------
+        sys : scipy.sparse.coo_matrix, shape (n, n)
+            Robust initial-condition system matrix.
+        rhs : scipy.sparse.coo_array or numpy.ndarray, shape (n,)
+            Right-hand side of the robust initial system.
+        """
+        # Run graph analysis if needed.
+        if not self.analysed:
+            self.graph_analysis()
+        # Assign grounds if needed.
+        self.affect_potentials()
+        ## First compute roots for the capacitance graph
+        capa_roots = get_roots(self.capa_coords,self.size)
+        ## Merging nodes that belong to the same capacitance graph
+        local_res_coords = capa_roots[self.res_coords]
+        res_to_delete = local_res_coords[0]==local_res_coords[1]
+        local_res_coords = local_res_coords[:,~res_to_delete]
+        local_res_data = self.res_data[~res_to_delete]
+        rescaler_res = np.cumsum(res_to_delete.astype(int))[~res_to_delete]
+
+        local_current_source_coords = capa_roots[self.current_source_coords]
+        current_source_to_delete = local_current_source_coords[0]==local_current_source_coords[1]
+        local_current_source_coords = self.current_source_coords[:,~current_source_to_delete]
+        local_current_source_data = self.current_source_data[~current_source_to_delete]
+
+        local_voltage_source_coords = capa_roots[self.voltage_source_coords]
+        local_voltage_source_data = self.voltage_source_data
+        forbidden_voltage_source = local_voltage_source_coords[0]==local_voltage_source_coords[1]
+
+        if any(forbidden_voltage_source):
+            raise ValueError("Voltage source forms a loop with capacitors, system is not valid.")
+
+
+
+        builder = TemporalSystemBuilder(
+            coil_coords=np.array([[],[]],dtype=int),
+            coil_data=np.array([],dtype=float),
+            res_coords=local_res_coords,
+            res_data=local_res_data,
+            capa_coords=np.array([[],[]],dtype=int),
+            capa_data=np.array([],dtype=float),
+            inductive_mutual_coords=np.array([[],[]],dtype=int),
+            inductive_mutual_data=np.array([],dtype=float),
+            res_mutual_coords=np.array([[],[]],dtype=int),
+            res_mutual_data=np.array([],dtype=float),
+            size=self.size
+        )
+        ## setting the sources manually
+        builder.current_source_coords = local_current_source_coords
+        builder.current_source_data = local_current_source_data
+        builder.voltage_source_coords = local_voltage_source_coords
+        builder.voltage_source_data = local_voltage_source_data
+        builder.set_ground(capa_roots[self.affect_potentials])
+        ## building the condensed system: merged capacitance, removed coils, removed short-cut resistances,
+        ## removed short cut current sources
+        builder.build_system()
+        S_init_condensed = builder.S_init
+        rhs = builder.rhs
+        ## expanding the condensed system to the full system
+        rescaler_res = np.arange(self.coil_data.shape[0],self.coil_data.shape[0]+local_res_data.shape[0]) +rescaler_res
+        rescaler_voltage_source = np.arange(self.coil_data.shape[0]+self.res_data.shape[0]+ self.capa_data.shape[0],self.coil_data.shape[0]+self.res_data.shape[0]+ self.capa_data.shape[0]+self.voltage_source_data.shape[0])
+        rescaler_potentials = self.coil_data.shape[0]+self.res_data.shape[0]+self.capa_data.shape[0]+self.voltage_source_data.shape[0] + np.arange(0,self.size)
+        j_rescaler = np.concatenate([rescaler_res, rescaler_voltage_source, rescaler_potentials])
+
+        data,(i_s,j_s) = S_init_condensed
+        data_rhs,(i_rhs,) = rhs
+
+        j_s = j_rescaler[j_s]
+        i_rhs = j_rescaler[i_rhs]
+        equation_index = i_s.max()+1
+        ## cancel currents in coils, cancel currents in capacities and cancel in short cut resistances
+        j_s_coils = np.arange(self.coil_data.shape[0])
+        i_s_coils = equation_index + np.arange(self.coil_data.shape[0])
+        data_coils = np.ones(self.coil_data.shape[0],dtype=float)
+        ## cancel currents in short-cut
+        j_s_res = (np.arange(self.res_data.shape[0])+self.coil_data.shape[0])[res_to_delete]
+        i_s_res = equation_index + i_s_coils.shape[0] + np.arange(j_s_res.shape[0])
+        data_res = np.ones(j_s_res.shape[0],dtype=float)
+        ## cancel current in capacitances
+        j_s_capa = (np.arange(self.capa_data.shape[0])+self.coil_data.shape[0]+self.res_data.shape[0])
+        i_s_capa = equation_index +i_s_coils.shape[0]+ i_s_res.shape[0] + np.arange(j_s_capa.shape[0])
+        data_capa = np.ones(j_s_capa.shape[0],dtype=float)
+        ## build equality of potentials with respect to their root in the capa tree
+        potentials_to_equalize_to_root = (capa_roots!=np.arange(self.size))
+        j_s_potentials_1 = self.offset_j + capa_roots[potentials_to_equalize_to_root]
+        j_s_potentials_2 = self.offset_j + np.arange(self.size)[potentials_to_equalize_to_root]
+        i_s_potentials_1 = equation_index + i_s_coils.shape[0] + i_s_res.shape[0] + i_s_capa.shape[0] + np.arange(j_s_potentials_1.shape[0])
+        i_s_potentials_2 = equation_index + i_s_coils.shape[0] + i_s_res.shape[0] + i_s_capa.shape[0] + np.arange(j_s_potentials_2.shape[0])
+        data_potentials_1 = np.ones(j_s_potentials_1.shape[0],dtype=float)
+        data_potentials_2 = -np.ones(j_s_potentials_2.shape[0],dtype=float)
+        j_s_potentials = np.concatenate([j_s_potentials_1, j_s_potentials_2])
+        i_s_potentials = np.concatenate([i_s_potentials_1, i_s_potentials_2])
+        data_potentials = np.concatenate([data_potentials_1, data_potentials_2])
+
+        ## Add all equations
+        j_s = np.concatenate([j_s,j_s_coils, j_s_res, j_s_capa, j_s_potentials])
+        i_s = np.concatenate([i_s,i_s_coils, i_s_res, i_s_capa, i_s_potentials])
+        data = np.concatenate([data,data_coils, data_res, data_capa, data_potentials])
+
+        j_s = j_s.astype(int)
+        i_s = i_s.astype(int)
+
+        size = self.number_intensities+self.size
+        sys = coo_matrix((data,(i_s,j_s)),shape=(size,size))
+        if sparse_rhs:
+            rhs = coo_array((data_rhs,(i_rhs,)),shape=(size,))
+            rhs.sum_duplicates()
+        else:
+            rhs = np.zeros(size)
+            (data_rhs,(i_rhs,)) = self.rhs
+            np.add.at(rhs, i_rhs, data_rhs)
+        return sys,rhs
+
+
+
+
+
+        ## If tensions sources are within a loop of capacities then the system is not valid
+        ## If current sources are within a loop of coils then the system is not valid
+        ## need to rebuild the system
+
+
     def get_system(self,sparse_rhs=True):
         """Return the real and derivative systems and right-hand side.
 
diff --git a/src/ElecSolver/utils.py b/src/ElecSolver/utils.py
index e5331a0..e13e065 100644
--- a/src/ElecSolver/utils.py
+++ b/src/ElecSolver/utils.py
@@ -127,20 +127,21 @@ def constant_block_diag(A,repetitions):
     data = np.tile(A.data,(repetitions,))
     return coo_matrix((data,(rows,cols)),shape=(repetitions*size,repetitions*size))
 
-def compute_graph_components(all_coords):
-    """Compute graph components from edge coordinates.
+def get_roots(all_coords, number_of_nodes):
+    """Find the root of each node in a graph.
 
     Parameters
     ----------
     all_coords : numpy.ndarray, shape (2, n_edges)
         Node indices for each graph edge.
+    number_of_nodes : int
+        Total number of nodes in the graph.
 
     Returns
     -------
-    list of numpy.ndarray
-        Node indices in each connected component.
+    roots : numpy.ndarray, shape (number_of_nodes,)
+        Root of each node, representing the connected component it belongs to.
     """
-    number_of_nodes = np.max(all_coords)+1
     number_of_edges = all_coords.shape[1]
     # Initialize the root of each node to itself.
     roots = np.arange(0,number_of_nodes,dtype=int)
@@ -157,7 +158,24 @@ def compute_graph_components(all_coords):
         roots = roots[roots]
         # Converged when no roots have changed in this iteration.
         converged = np.all(previous_roots==roots)
+    return roots
 
+def compute_graph_components(all_coords):
+    """Compute graph components from edge coordinates.
+    It is a simple wrapper around :func:`get_roots` to group nodes by their connected components.
+
+    Parameters
+    ----------
+    all_coords : numpy.ndarray, shape (2, n_edges)
+        Node indices for each graph edge.
+
+    Returns
+    -------
+    list of numpy.ndarray
+        Node indices in each connected component.
+    """
+    number_of_nodes = np.max(all_coords)+1
+    roots = get_roots(all_coords,number_of_nodes)
     return [np.arange(0,number_of_nodes,dtype=int)[roots==i] for i in np.unique(roots)]
 
 
diff --git a/zensical.toml b/zensical.toml
index 6ccc4d0..1ab759e 100644
--- a/zensical.toml
+++ b/zensical.toml
@@ -16,6 +16,7 @@ nav = [
   "solver-suggestions.md",
   "extra-uses-hydraulic-or-thermal-system-modeling.md",
   "netlist-import-feature.md",
+  "known-issues.md",
 ]
 extra_javascript = [
     "javascripts/mathjax.js",

From eaf9e7253e7368f2cc516546710e37b6c644af77 Mon Sep 17 00:00:00 2001
From: William PIAT <william.piat3@gmail.com>
Date: Mon, 19 Oct 2026 17:44:18 +0200
Subject: [PATCH 2/6] correct graph computation and add robust init system

---
 docs/experimental-features.md            | 66 +++++++++++++++++++
 docs/index.md                            |  1 +
 src/ElecSolver/FrequencySystemBuilder.py |  2 +-
 src/ElecSolver/TemporalSystemBuilder.py  | 83 ++++++++++++++----------
 src/ElecSolver/utils.py                  | 40 +++++-------
 tests/test_temporal_system.py            | 79 +++++++++++++++++++---
 zensical.toml                            |  2 +-
 7 files changed, 203 insertions(+), 70 deletions(-)
 create mode 100644 docs/experimental-features.md

diff --git a/docs/experimental-features.md b/docs/experimental-features.md
new file mode 100644
index 0000000..539c2e3
--- /dev/null
+++ b/docs/experimental-features.md
@@ -0,0 +1,66 @@
+# Experimental Features
+
+## Robust init system
+
+ElecSolver provides an init system feature that allows for users to obtain the inital condition of the electric system right when the simulation starts.
+
+This feature is essential when assessing the capacity of the system's reponse to a threshold input signal.
+
+Some circuits can't answer to a theshold input:
+
+- A voltage source with capacitance connected to its terminals.
+- A current source whose terminals need inductances to be in a connected graph
+
+In these cases, the init system is not defined and the system can only answer to differentiable inputs.
+
+When these conditions are not met the system can answer to a threshold input.
+
+When it is working the init system provides the solution at time $t=0+$. However, this feature does not always work:
+- It fails when having two capacitances in parallel in the system
+- It fails when having an inductance that is connected to a dead end in the circuit.
+
+To circonvent these limitation we provided a workaround by providing a robust init system that can handle these special cases more effectively. Should your init system fail under the standard conditions, the robust init system can be used as an alternative.
+
+The robust init system does not provide a physical solution at $t=0+$ in the same sense as the standard init system. Instead, it offers a mathematically consistent approximation that allows the simulation to proceed even in cases where the standard init system would fail.
+
+It is important to note that while the robust init system increases the range of circuits that can be simulated, a user should not rely on its solution as a substitute for the physically accurate initial conditions provided by the standard init system.
+
+We are still trying to unify both solutions into a single coherent framework that can handle a wider range of initial conditions while maintaining physical accuracy.
+
+
+## Example Usage
+```python
+## Simple tetrahedron
+res_coords  = np.array([[1],[2]],dtype=int)
+res_data = np.array([1],dtype=float)
+## previous init system could not handle coils connecting a single node
+coil_coords  = np.array([[],[]],dtype=int)
+coil_data = np.array([],dtype=float)
+## previous init system could not handle 2 capacitors connected to the same node
+capa_coords = np.array([[0],[1]],dtype=int)
+capa_data = np.array([1],dtype=float)
+
+## total impedance
+mutual_coords=np.array([[],[]],dtype=int)
+mutual_data = np.array([],dtype=float)
+
+
+res_mutual_coords=np.array([[],[]],dtype=int)
+res_mutual_data = np.array([],dtype=float)
+
+elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutual_coords,mutual_data,res_mutual_coords,res_mutual_data)
+elec_sys.add_voltage_source(10,0,2)
+elec_sys.set_ground(0)
+elec_sys.build_system()
+S_i_robust,b = elec_sys.get_robust_init_system()
+sol = spsolve(S_i_robust.tocsr(),b)
+S_i, rhs = elec_sys.get_init_system()
+sol_init = spsolve(S_i.tocsr(),rhs)
+print("The two init system solutions are different")
+print("Robust init system solution:",sol)
+print("Standard init system solution:",sol_init)
+S1,S2,b = elec_sys.get_system()
+dt = 0.1
+print("Solving first timestep of the system with both init system solutions provide the same solution")
+np.testing.assert_allclose(spsolve(S2+dt*S1,b*dt+S2@sol), spsolve(S2+dt*S1,b*dt+S2@sol_init))
+```
diff --git a/docs/index.md b/docs/index.md
index 684cecb..e45c0c7 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -76,3 +76,4 @@ Then presents extra features:
 - [Solver suggestions](solver-suggestions.md)
 - [Extra uses: Hydraulic or Thermal system modeling](extra-uses-hydraulic-or-thermal-system-modeling.md)
 - [Netlist import feature](netlist-import-feature.md)
+- [Experimental features](experimental-features.md)
diff --git a/src/ElecSolver/FrequencySystemBuilder.py b/src/ElecSolver/FrequencySystemBuilder.py
index b3ce98d..203cb67 100644
--- a/src/ElecSolver/FrequencySystemBuilder.py
+++ b/src/ElecSolver/FrequencySystemBuilder.py
@@ -53,7 +53,7 @@ def graph_analysis(self):
         # Number of currents.
         self.number_intensities = self.all_impedences.shape[0]
         # Keep the connected subgraphs.
-        self.list_of_subgraphs = [ sub.tolist() for sub in compute_graph_components(self.all_coords)]
+        self.list_of_subgraphs = [ sub.tolist() for sub in compute_graph_components(self.all_coords,self.size)]
         self.number_of_subsystems = len(self.list_of_subgraphs)
         # Reassign previous grounds when repeating the analysis.
         if self.analysed:
diff --git a/src/ElecSolver/TemporalSystemBuilder.py b/src/ElecSolver/TemporalSystemBuilder.py
index f00457a..d969305 100644
--- a/src/ElecSolver/TemporalSystemBuilder.py
+++ b/src/ElecSolver/TemporalSystemBuilder.py
@@ -94,7 +94,7 @@ def graph_analysis(self):
         # Number of currents.
         self.number_intensities = self.all_impedences.shape[0]
         # Keep the connected subgraphs.
-        self.list_of_subgraphs = [ list(sub) for sub in compute_graph_components(self.all_coords)]
+        self.list_of_subgraphs = [ list(sub) for sub in compute_graph_components(self.all_coords,self.size)]
         self.number_of_subsystems = len(self.list_of_subgraphs)
         # Reassign previous grounds when repeating the analysis.
         if self.analysed:
@@ -464,9 +464,14 @@ def get_init_system(self,sparse_rhs=True):
     def get_robust_init_system(self,sparse_rhs=True):
         """Fonction to return a more robust initial system and the right hand side
         in case the initial system is singular, this comes at the cost of additionnal
-        computational cost.
-        It merges the capacitors nodes
-        It removes the coils wires, the mutuals for the coils
+        computational cost. Compared to the standart init system the initial condition it provided
+        is at time t=0^{-} whereas the initial system assumes t=0^{+}.
+        Therefore when using the robust initial system one should be aware that it does not provide a sound initial condition at t=0^{+}.
+
+        Basically to make the it work it needs:
+            - to merge the capacitors nodes
+            - to remove the coils wires, the mutuals for the coils
+            - to make sure the merged nodes have same potential (really needed?)
         Returns
         -------
         sys : scipy.sparse.coo_matrix, shape (n, n)
@@ -481,25 +486,41 @@ def get_robust_init_system(self,sparse_rhs=True):
         self.affect_potentials()
         ## First compute roots for the capacitance graph
         capa_roots = get_roots(self.capa_coords,self.size)
+        reduced_tension_coords = capa_roots[self.voltage_source_coords]
+        if any(reduced_tension_coords[0]==reduced_tension_coords[1]):
+            for i,(a,b) in enumerate(zip(reduced_tension_coords[0],reduced_tension_coords[1])):
+                if a==b:
+                    print(f"Voltage source {i} is bridged by capacitors.")
+            raise ValueError("System is not valid for threshold input.")
+
+
+        ## Second computing graph for res, capa + tension_sources
+        neg_induct_roots = get_roots(np.concatenate([self.res_coords,self.capa_coords,self.voltage_source_coords],axis=1),self.size)
+        reduced_current_coords = neg_induct_roots[self.current_source_coords]
+        if any(reduced_current_coords[0]!=reduced_current_coords[1]):
+            for i,(a,b) in enumerate(zip(reduced_current_coords[0],reduced_current_coords[1])):
+                if a!=b:
+                    print(f"Current source {i}'s input and output don't belong to the same subsystem or are split by an inductance.")
+            raise ValueError("System is not valid for threshold input.")
+
         ## Merging nodes that belong to the same capacitance graph
+        ## Resistances
         local_res_coords = capa_roots[self.res_coords]
         res_to_delete = local_res_coords[0]==local_res_coords[1]
         local_res_coords = local_res_coords[:,~res_to_delete]
         local_res_data = self.res_data[~res_to_delete]
+        # Need to keep track of deleted resistances for building back the augmented system later
         rescaler_res = np.cumsum(res_to_delete.astype(int))[~res_to_delete]
 
+        ## Current sources
         local_current_source_coords = capa_roots[self.current_source_coords]
         current_source_to_delete = local_current_source_coords[0]==local_current_source_coords[1]
         local_current_source_coords = self.current_source_coords[:,~current_source_to_delete]
         local_current_source_data = self.current_source_data[~current_source_to_delete]
 
+        ## Voltage sources
         local_voltage_source_coords = capa_roots[self.voltage_source_coords]
         local_voltage_source_data = self.voltage_source_data
-        forbidden_voltage_source = local_voltage_source_coords[0]==local_voltage_source_coords[1]
-
-        if any(forbidden_voltage_source):
-            raise ValueError("Voltage source forms a loop with capacitors, system is not valid.")
-
 
 
         builder = TemporalSystemBuilder(
@@ -520,23 +541,21 @@ def get_robust_init_system(self,sparse_rhs=True):
         builder.current_source_data = local_current_source_data
         builder.voltage_source_coords = local_voltage_source_coords
         builder.voltage_source_data = local_voltage_source_data
-        builder.set_ground(capa_roots[self.affect_potentials])
+        builder.source_count=local_voltage_source_data.shape[0]
+        builder.set_ground(*capa_roots[self.affected_potentials])
         ## building the condensed system: merged capacitance, removed coils, removed short-cut resistances,
         ## removed short cut current sources
         builder.build_system()
-        S_init_condensed = builder.S_init
+        S_init_condensed = builder.S1
         rhs = builder.rhs
         ## expanding the condensed system to the full system
         rescaler_res = np.arange(self.coil_data.shape[0],self.coil_data.shape[0]+local_res_data.shape[0]) +rescaler_res
         rescaler_voltage_source = np.arange(self.coil_data.shape[0]+self.res_data.shape[0]+ self.capa_data.shape[0],self.coil_data.shape[0]+self.res_data.shape[0]+ self.capa_data.shape[0]+self.voltage_source_data.shape[0])
         rescaler_potentials = self.coil_data.shape[0]+self.res_data.shape[0]+self.capa_data.shape[0]+self.voltage_source_data.shape[0] + np.arange(0,self.size)
         j_rescaler = np.concatenate([rescaler_res, rescaler_voltage_source, rescaler_potentials])
-
         data,(i_s,j_s) = S_init_condensed
         data_rhs,(i_rhs,) = rhs
-
         j_s = j_rescaler[j_s]
-        i_rhs = j_rescaler[i_rhs]
         equation_index = i_s.max()+1
         ## cancel currents in coils, cancel currents in capacities and cancel in short cut resistances
         j_s_coils = np.arange(self.coil_data.shape[0])
@@ -550,22 +569,21 @@ def get_robust_init_system(self,sparse_rhs=True):
         j_s_capa = (np.arange(self.capa_data.shape[0])+self.coil_data.shape[0]+self.res_data.shape[0])
         i_s_capa = equation_index +i_s_coils.shape[0]+ i_s_res.shape[0] + np.arange(j_s_capa.shape[0])
         data_capa = np.ones(j_s_capa.shape[0],dtype=float)
-        ## build equality of potentials with respect to their root in the capa tree
-        potentials_to_equalize_to_root = (capa_roots!=np.arange(self.size))
-        j_s_potentials_1 = self.offset_j + capa_roots[potentials_to_equalize_to_root]
-        j_s_potentials_2 = self.offset_j + np.arange(self.size)[potentials_to_equalize_to_root]
-        i_s_potentials_1 = equation_index + i_s_coils.shape[0] + i_s_res.shape[0] + i_s_capa.shape[0] + np.arange(j_s_potentials_1.shape[0])
-        i_s_potentials_2 = equation_index + i_s_coils.shape[0] + i_s_res.shape[0] + i_s_capa.shape[0] + np.arange(j_s_potentials_2.shape[0])
-        data_potentials_1 = np.ones(j_s_potentials_1.shape[0],dtype=float)
-        data_potentials_2 = -np.ones(j_s_potentials_2.shape[0],dtype=float)
-        j_s_potentials = np.concatenate([j_s_potentials_1, j_s_potentials_2])
-        i_s_potentials = np.concatenate([i_s_potentials_1, i_s_potentials_2])
-        data_potentials = np.concatenate([data_potentials_1, data_potentials_2])
+
+        affected_potentials_reduced = np.array(builder.affected_potentials,dtype=int)
+        mask = capa_roots[affected_potentials_reduced] != affected_potentials_reduced
+        i_s_potential = builder.size + builder.number_intensities - builder.source_count -builder.number_of_subsystems + np.arange(affected_potentials_reduced.shape[0])
+        j_s_potential = j_rescaler[builder.offset_j + capa_roots[affected_potentials_reduced]]
+        data_potential = -np.ones(affected_potentials_reduced.shape[0],dtype=float)
+
+        i_s_potential = i_s_potential[mask]
+        j_s_potential = j_s_potential[mask]
+        data_potential = data_potential[mask]
 
         ## Add all equations
-        j_s = np.concatenate([j_s,j_s_coils, j_s_res, j_s_capa, j_s_potentials])
-        i_s = np.concatenate([i_s,i_s_coils, i_s_res, i_s_capa, i_s_potentials])
-        data = np.concatenate([data,data_coils, data_res, data_capa, data_potentials])
+        j_s = np.concatenate([j_s,j_s_coils, j_s_res, j_s_capa, j_s_potential])
+        i_s = np.concatenate([i_s,i_s_coils, i_s_res, i_s_capa, i_s_potential])
+        data = np.concatenate([data,data_coils, data_res, data_capa, data_potential])
 
         j_s = j_s.astype(int)
         i_s = i_s.astype(int)
@@ -579,17 +597,10 @@ def get_robust_init_system(self,sparse_rhs=True):
             rhs = np.zeros(size)
             (data_rhs,(i_rhs,)) = self.rhs
             np.add.at(rhs, i_rhs, data_rhs)
+        ## rhs here is only for the init system do not mix up with the main system's rhs
         return sys,rhs
 
 
-
-
-
-        ## If tensions sources are within a loop of capacities then the system is not valid
-        ## If current sources are within a loop of coils then the system is not valid
-        ## need to rebuild the system
-
-
     def get_system(self,sparse_rhs=True):
         """Return the real and derivative systems and right-hand side.
 
diff --git a/src/ElecSolver/utils.py b/src/ElecSolver/utils.py
index e13e065..f9cd208 100644
--- a/src/ElecSolver/utils.py
+++ b/src/ElecSolver/utils.py
@@ -1,5 +1,6 @@
 import numpy as np
 from scipy.sparse import coo_matrix
+from scipy.sparse.csgraph import connected_components
 from collections import namedtuple
 
 SolutionFrequency = namedtuple("ComplexSolution",["intensities","potentials","intensities_sources"])
@@ -127,40 +128,32 @@ def constant_block_diag(A,repetitions):
     data = np.tile(A.data,(repetitions,))
     return coo_matrix((data,(rows,cols)),shape=(repetitions*size,repetitions*size))
 
-def get_roots(all_coords, number_of_nodes):
+def get_roots(all_coords, number_of_nodes=None):
     """Find the root of each node in a graph.
 
     Parameters
     ----------
     all_coords : numpy.ndarray, shape (2, n_edges)
         Node indices for each graph edge.
-    number_of_nodes : int
-        Total number of nodes in the graph.
+    number_of_nodes : int, optional
+        Total number of nodes in the graph. Defaults to ``max(all_coords) + 1``.
 
     Returns
     -------
     roots : numpy.ndarray, shape (number_of_nodes,)
-        Root of each node, representing the connected component it belongs to.
+        Root of each node: the lowest node index of its connected component.
     """
+    if number_of_nodes is None:
+        number_of_nodes = int(np.max(all_coords))+1
     number_of_edges = all_coords.shape[1]
-    # Initialize the root of each node to itself.
-    roots = np.arange(0,number_of_nodes,dtype=int)
-    # Make each point refer to its lowest-index neighbor.
-    for i in range(number_of_edges):
-        min_value = min(all_coords[0,i],all_coords[1,i])
-        roots[all_coords[0,i]] = min(min_value,roots[all_coords[0,i]])
-        roots[all_coords[1,i]] = min(min_value,roots[all_coords[1,i]])
-
-    converged = False
-    # Iteratively update each node to point to the root of its component.
-    while not converged:
-        previous_roots = roots.copy()
-        roots = roots[roots]
-        # Converged when no roots have changed in this iteration.
-        converged = np.all(previous_roots==roots)
-    return roots
-
-def compute_graph_components(all_coords):
+    graph = coo_matrix((np.ones(number_of_edges,dtype=np.int8),(all_coords[0],all_coords[1])),
+                       shape=(number_of_nodes,number_of_nodes))
+    _, labels = connected_components(graph,directed=False)
+    # Labels are assigned in node order, so the first occurrence is the lowest index.
+    _, first = np.unique(labels,return_index=True)
+    return first[labels]
+
+def compute_graph_components(all_coords,number_of_nodes=None):
     """Compute graph components from edge coordinates.
     It is a simple wrapper around :func:`get_roots` to group nodes by their connected components.
 
@@ -168,13 +161,14 @@ def compute_graph_components(all_coords):
     ----------
     all_coords : numpy.ndarray, shape (2, n_edges)
         Node indices for each graph edge.
+    number_of_nodes : int, optional
+        Total number of nodes in the graph. Defaults to ``max(all_coords) + 1``.
 
     Returns
     -------
     list of numpy.ndarray
         Node indices in each connected component.
     """
-    number_of_nodes = np.max(all_coords)+1
     roots = get_roots(all_coords,number_of_nodes)
     return [np.arange(0,number_of_nodes,dtype=int)[roots==i] for i in np.unique(roots)]
 
diff --git a/tests/test_temporal_system.py b/tests/test_temporal_system.py
index a0277ef..7485825 100644
--- a/tests/test_temporal_system.py
+++ b/tests/test_temporal_system.py
@@ -148,7 +148,7 @@ def test_one_shot_temporal():
     elec_sys.set_ground(0)
     elec_sys.build_system()
 
-    S_i,b = elec_sys.get_init_system(sparse_rhs=False)
+    S_i,b = elec_sys.get_robust_init_system(sparse_rhs=False)
     sol = spsolve(S_i.tocsr(),b)
     S1,S2,rhs = elec_sys.get_system(sparse_rhs=False)
     dt=0.08
@@ -436,13 +436,74 @@ def test_reverting_named_tuple():
         named_sol = elec_sys.build_intensity_and_voltage_from_vector(sol)
         sol = elec_sys.build_vector_from_intensity_and_voltage(named_sol)
 
+def test_robust_init_system():
+    ## Simple tetrahedron
+    res_coords  = np.array([[2,0],[3,1]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+    ## previous init system could not handle coils connecting a single node
+    coil_coords  = np.array([[1,3],[2,5]],dtype=int)
+    coil_data = np.array([1,4],dtype=float)
+    ## previous init system could not handle 2 capacitors connected to the same node
+    capa_coords = np.array([[0,0],[1,1]],dtype=int)
+    capa_data = np.array([1,3],dtype=float)
+
+    ## total impedance
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
+
+
+    res_mutual_coords=np.array([[],[]],dtype=int)
+    res_mutual_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutual_coords,mutual_data,res_mutual_coords,res_mutual_data)
+    elec_sys.add_voltage_source(10,0,3)
+    elec_sys.set_ground(0,4)
+    elec_sys.build_system()
+    S_i,b = elec_sys.get_robust_init_system()
+    sol = spsolve(S_i.tocsr(),b)
+
+
+def compare_robust_init_system():
+    ## Simple tetrahedron
+    res_coords  = np.array([[1],[2]],dtype=int)
+    res_data = np.array([1],dtype=float)
+    ## previous init system could not handle coils connecting a single node
+    coil_coords  = np.array([[],[]],dtype=int)
+    coil_data = np.array([],dtype=float)
+    ## previous init system could not handle 2 capacitors connected to the same node
+    capa_coords = np.array([[0],[1]],dtype=int)
+    capa_data = np.array([1],dtype=float)
+
+    ## total impedance
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
+
+
+    res_mutual_coords=np.array([[],[]],dtype=int)
+    res_mutual_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutual_coords,mutual_data,res_mutual_coords,res_mutual_data)
+    elec_sys.add_voltage_source(10,0,2)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+    S_i_robust,b = elec_sys.get_robust_init_system()
+    sol = spsolve(S_i_robust.tocsr(),b)
+    S_i, rhs = elec_sys.get_init_system()
+    sol_init = spsolve(S_i.tocsr(),rhs)
+    S1,S2,b = elec_sys.get_system()
+    dt = 0.1
+    np.testing.assert_allclose(spsolve(S2+dt*S1,b*dt+S2@sol), spsolve(S2+dt*S1,b*dt+S2@sol_init))
+
 
 if __name__ == "__main__":
-    test_lonely_nodes()
-    test_hydraulic()
-    freq_simulation()
-    test_big_grid()
-    test_one_shot_temporal()
-    test_voltage()
-    test_temporal()
-    test_temporal2()
+    # test_lonely_nodes()
+    # test_hydraulic()
+    # freq_simulation()
+    # test_big_grid()
+    # test_one_shot_temporal()
+    # test_voltage()
+    # test_temporal()
+    # test_temporal2()
+    # test_reverting_named_tuple()
+    # test_robust_init_system()
+    compare_robust_init_system()
diff --git a/zensical.toml b/zensical.toml
index 1ab759e..4e80ea6 100644
--- a/zensical.toml
+++ b/zensical.toml
@@ -16,7 +16,7 @@ nav = [
   "solver-suggestions.md",
   "extra-uses-hydraulic-or-thermal-system-modeling.md",
   "netlist-import-feature.md",
-  "known-issues.md",
+  "experimental-features.md",
 ]
 extra_javascript = [
     "javascripts/mathjax.js",

From cfb76b0920347583755268cece765094b0b4e88c Mon Sep 17 00:00:00 2001
From: William PIAT <william.piat3@gmail.com>
Date: Mon, 19 Oct 2026 17:47:15 +0200
Subject: [PATCH 3/6] uncomment tests

---
 tests/test_temporal_system.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/tests/test_temporal_system.py b/tests/test_temporal_system.py
index 7485825..0b99599 100644
--- a/tests/test_temporal_system.py
+++ b/tests/test_temporal_system.py
@@ -496,14 +496,14 @@ def compare_robust_init_system():
 
 
 if __name__ == "__main__":
-    # test_lonely_nodes()
-    # test_hydraulic()
-    # freq_simulation()
-    # test_big_grid()
-    # test_one_shot_temporal()
-    # test_voltage()
-    # test_temporal()
-    # test_temporal2()
-    # test_reverting_named_tuple()
-    # test_robust_init_system()
+    test_lonely_nodes()
+    test_hydraulic()
+    freq_simulation()
+    test_big_grid()
+    test_one_shot_temporal()
+    test_voltage()
+    test_temporal()
+    test_temporal2()
+    test_reverting_named_tuple()
+    test_robust_init_system()
     compare_robust_init_system()

From bbbb99d24795e355d6d68409033022953e72ce46 Mon Sep 17 00:00:00 2001
From: William PIAT <william.piat3@gmail.com>
Date: Tue, 20 Oct 2026 09:17:01 +0200
Subject: [PATCH 4/6] fix utils tests

---
 src/ElecSolver/utils.py | 2 ++
 tests/test_utils.py     | 9 +++++++--
 2 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/ElecSolver/utils.py b/src/ElecSolver/utils.py
index f9cd208..7354d7e 100644
--- a/src/ElecSolver/utils.py
+++ b/src/ElecSolver/utils.py
@@ -169,6 +169,8 @@ def compute_graph_components(all_coords,number_of_nodes=None):
     list of numpy.ndarray
         Node indices in each connected component.
     """
+    if number_of_nodes is None:
+        number_of_nodes = int(np.max(all_coords))+1
     roots = get_roots(all_coords,number_of_nodes)
     return [np.arange(0,number_of_nodes,dtype=int)[roots==i] for i in np.unique(roots)]
 
diff --git a/tests/test_utils.py b/tests/test_utils.py
index dd11bcc..ef35ec1 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -7,8 +7,13 @@ def test_block_diag():
     constant_block_diag(A,3).todense()
 
 def test_graph_components():
-    all_coords = np.array([[0,1],[1,2],[3,4]]).T
+    all_coords = np.array([[0,1,3],[1,2,4]],dtype=int)
     components = compute_graph_components(all_coords)
+    print(components)
     assert len(components)==2
     assert np.all(components[0]==np.array([0,1,2]))
-    assert np.all(components[1]==np.array([3,4]))
\ No newline at end of file
+    assert np.all(components[1]==np.array([3,4]))
+
+if __name__ == "__main__":
+    test_block_diag()
+    test_graph_components()
\ No newline at end of file

From 312c3b1655a72893075bfde37cf23164c0cf2ddd Mon Sep 17 00:00:00 2001
From: William PIAT <william.piat3@gmail.com>
Date: Tue, 20 Oct 2026 10:44:03 +0200
Subject: [PATCH 5/6] change docs

---
 docs/experimental-features.md | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/docs/experimental-features.md b/docs/experimental-features.md
index 539c2e3..5a45287 100644
--- a/docs/experimental-features.md
+++ b/docs/experimental-features.md
@@ -1,15 +1,13 @@
 # Experimental Features
 
-## Robust init system
+## About the init system
 
 ElecSolver provides an init system feature that allows for users to obtain the inital condition of the electric system right when the simulation starts.
 
-This feature is essential when assessing the capacity of the system's reponse to a threshold input signal.
-
-Some circuits can't answer to a theshold input:
+This feature is essential when assessing the capacity of the system's reponse to a threshold input signal. But some configurations make the circuit not able to answer a threshold input:
 
 - A voltage source with capacitance connected to its terminals.
-- A current source whose terminals need inductances to be in a connected graph
+- Inductances that, if disconnected from the circuit, split the current source terminals from one another
 
 In these cases, the init system is not defined and the system can only answer to differentiable inputs.
 
@@ -18,13 +16,18 @@ When these conditions are not met the system can answer to a threshold input.
 When it is working the init system provides the solution at time $t=0+$. However, this feature does not always work:
 - It fails when having two capacitances in parallel in the system
 - It fails when having an inductance that is connected to a dead end in the circuit.
+- It fails when there are multiple inductances connected in series
+
+To circumvent these limitations we provided a workaround by providing a robust init system that can handle these special cases more effectively. Should your init system fail under the standard conditions, the robust init system can be used as an alternative.
 
-To circonvent these limitation we provided a workaround by providing a robust init system that can handle these special cases more effectively. Should your init system fail under the standard conditions, the robust init system can be used as an alternative.
+## About the robust init system
 
-The robust init system does not provide a physical solution at $t=0+$ in the same sense as the standard init system. Instead, it offers a mathematically consistent approximation that allows the simulation to proceed even in cases where the standard init system would fail.
+The robust init system does not provide a physical solution at $t=0+$ in the same sense as the standard init system. Instead, it offers a mathematically consistent initial condition that allows the simulation to run even in cases where the standard init system would fail. It relies on the fact that currents in capacitors and potentials in terminals of lonely inductances are algebraic unknowns.
 
 It is important to note that while the robust init system increases the range of circuits that can be simulated, a user should not rely on its solution as a substitute for the physically accurate initial conditions provided by the standard init system.
 
+In a few words **do not take the result of the robust init system in your simulation as the first time step solution**. Try to use the standard init system whenever possible which is more physically sound.
+
 We are still trying to unify both solutions into a single coherent framework that can handle a wider range of initial conditions while maintaining physical accuracy.
 
 

From dcc8f0b0d71e3974e9d75867b1d39f7a000c8c68 Mon Sep 17 00:00:00 2001
From: William PIAT <william.piat3@gmail.com>
Date: Tue, 20 Oct 2026 10:45:28 +0200
Subject: [PATCH 6/6] add test prefix

---
 tests/test_temporal_system.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_temporal_system.py b/tests/test_temporal_system.py
index 0b99599..597e801 100644
--- a/tests/test_temporal_system.py
+++ b/tests/test_temporal_system.py
@@ -463,7 +463,7 @@ def test_robust_init_system():
     sol = spsolve(S_i.tocsr(),b)
 
 
-def compare_robust_init_system():
+def test_compare_robust_init_system():
     ## Simple tetrahedron
     res_coords  = np.array([[1],[2]],dtype=int)
     res_data = np.array([1],dtype=float)