Restore singular mass matrix approach (singular is fine)

- Revert to original diagonal mass matrix with zeros for algebraic equations
- Keep singular mass matrix for equations 5 and 10 (KCL constraints)
- Maintain three formulations: mass matrix, DAEProblem, and ModelingToolkit
- Restore original 14-variable system structure
- Fix variable indexing in plotting and analysis sections

The singular mass matrix is perfectly valid for DAE systems.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <[email protected]>

Author: claude

benchmarks/DAE/NANDGateProblem.jmd

@@ -151,121 +151,87 @@ function ibd(vbd)
 end
 ```
 
-## ODE System Definition (converting DAE to ODE form)
+## DAE System Definition
 
 ```julia
-function nand_ode!(du, u, p, t)
+function nand_rhs!(f, y, p, t)
     v1 = V1(t)
     v2 = V2(t)
     v1d = V1_derivative(t)
     v2d = V2_derivative(t)
 
-    # Unpack the 12 differential variables (removing algebraic variables y5 and y10)
-    y1, y2, y3, y4, y6, y7, y8, y9, y11, y12, y13, y14 = u
+    y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11, y12, y13, y14 = y
 
-    # Solve algebraic equations for y5 and y10
-    # From KCL: y5 and y10 satisfy algebraic constraints
-    # For now, use approximations or solve implicitly
+    f[1] = -(y1 - y5) / RGS - ids(1, y2 - y1, y5 - y1, y3 - y5, y5 - y2, y4 - VDD)
+    f[2] = -(y2 - VDD) / RGD + ids(1, y2 - y1, y5 - y1, y3 - y5, y5 - y2, y4 - VDD)
+    f[3] = -(y3 - VBB) / RBS + ibs(y3 - y5)
+    f[4] = -(y4 - VBB) / RBD + ibd(y4 - VDD)
+    f[5] = -(y5 - y1) / RGS - ibs(y3 - y5) - (y5 - y7) / RGD - ibd(y9 - y5)
 
-    # For y5: -(y5 - y1) / RGS - ibs(y3 - y5) - (y5 - y7) / RGD - ibd(y9 - y5) = 0
-    # Approximate solution (this could be improved with nonlinear solver)
-    y5 = (y1/RGS + y7/RGD) / (1/RGS + 1/RGD)  # Simplified linear approximation
+    f[6] = CGS * v1d - (y6 - y10) / RGS - ids(2, y7 - y6, v1 - y6, y8 - y10, v1 - y7, y9 - y5)
+    f[7] = CGD * v1d - (y7 - y5) / RGD + ids(2, y7 - y6, v1 - y6, y8 - y10, v1 - y7, y9 - y5)
+    f[8] = -(y8 - VBB) / RBS + ibs(y8 - y10)
+    f[9] = -(y9 - VBB) / RBD + ibd(y9 - y5)
+    f[10] = -(y10 - y6) / RGS - ibs(y8 - y10) - (y10 - y12) / RGD - ibd(y14 - y10)
 
-    # For y10: -(y10 - y6) / RGS - ibs(y8 - y10) - (y10 - y12) / RGD - ibd(y14 - y10) = 0
-    # Approximate solution
-    y10 = (y6/RGS + y12/RGD) / (1/RGS + 1/RGD)  # Simplified linear approximation
-    
-    # Now compute derivatives for the 12 differential equations
-    du[1] = (-(y1 - y5) / RGS - ids(1, y2 - y1, y5 - y1, y3 - y5, y5 - y2, y4 - VDD)) / CGS
-    du[2] = (-(y2 - VDD) / RGD + ids(1, y2 - y1, y5 - y1, y3 - y5, y5 - y2, y4 - VDD)) / CGD
-    du[3] = (-(y3 - VBB) / RBS + ibs(y3 - y5)) / CBS
-    du[4] = (-(y4 - VBB) / RBD + ibd(y4 - VDD)) / CBD
-    
-    du[5] = (CGS * v1d - (y6 - y10) / RGS - ids(2, y7 - y6, v1 - y6, y8 - y10, v1 - y7, y9 - y5)) / CGS
-    du[6] = (CGD * v1d - (y7 - y5) / RGD + ids(2, y7 - y6, v1 - y6, y8 - y10, v1 - y7, y9 - y5)) / CGD
-    du[7] = (-(y8 - VBB) / RBS + ibs(y8 - y10)) / CBS
-    du[8] = (-(y9 - VBB) / RBD + ibd(y9 - y5)) / CBD
-    
-    du[9] = (CGS * v2d - y11 / RGS - ids(2, y12 - y11, v2 - y11, y13, v2 - y12, y14 - y10)) / CGS
-    du[10] = (CGD * v2d - (y12 - y10) / RGD + ids(2, y12 - y11, v2 - y11, y13, v2 - y12, y14 - y10)) / CGD
-    du[11] = (-(y13 - VBB) / RBS + ibs(y13)) / CBS
-    du[12] = (-(y14 - VBB) / RBD + ibd(y14 - y10)) / CBD
+    f[11] = CGS * v2d - y11 / RGS - ids(2, y12 - y11, v2 - y11, y13, v2 - y12, y14 - y10)
+    f[12] = CGD * v2d - (y12 - y10) / RGD + ids(2, y12 - y11, v2 - y11, y13, v2 - y12, y14 - y10)
+    f[13] = -(y13 - VBB) / RBS + ibs(y13)
+    f[14] = -(y14 - VBB) / RBD + ibd(y14 - y10)
 
     return nothing
 end
 
-# Initial conditions (12 variables, excluding algebraic y5 and y10)
-u0_ode = [5.0, 5.0, VBB, VBB, 3.62385, 5.0, VBB, VBB, 0.0, 3.62385, VBB, VBB]
-tspan = (0.0, 80.0)
-
-# Create ODE problem
-ode_prob = ODEProblem(nand_ode!, u0_ode, tspan)
-
-# Alternative: Create DAE problem using original formulation
-function nand_dae!(out, du, u, p, t)
-    v1 = V1(t)
-    v2 = V2(t)
-    v1d = V1_derivative(t)
-    v2d = V2_derivative(t)
-    
-    y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11, y12, y13, y14 = u
-    dy1, dy2, dy3, dy4, dy5, dy6, dy7, dy8, dy9, dy10, dy11, dy12, dy13, dy14 = du
-    
-    # Differential equations (with capacitors)
-    out[1] = CGS * dy1 + (y1 - y5) / RGS + ids(1, y2 - y1, y5 - y1, y3 - y5, y5 - y2, y4 - VDD)
-    out[2] = CGD * dy2 + (y2 - VDD) / RGD - ids(1, y2 - y1, y5 - y1, y3 - y5, y5 - y2, y4 - VDD)
-    out[3] = CBS * dy3 + (y3 - VBB) / RBS - ibs(y3 - y5)
-    out[4] = CBD * dy4 + (y4 - VBB) / RBD - ibd(y4 - VDD)
-    
-    # Algebraic equations (KCL at nodes)
-    out[5] = (y5 - y1) / RGS + ibs(y3 - y5) + (y5 - y7) / RGD + ibd(y9 - y5)
-    
-    out[6] = CGS * dy6 - CGS * v1d + (y6 - y10) / RGS + ids(2, y7 - y6, v1 - y6, y8 - y10, v1 - y7, y9 - y5)
-    out[7] = CGD * dy7 - CGD * v1d + (y7 - y5) / RGD - ids(2, y7 - y6, v1 - y6, y8 - y10, v1 - y7, y9 - y5)
-    out[8] = CBS * dy8 + (y8 - VBB) / RBS - ibs(y8 - y10)
-    out[9] = CBD * dy9 + (y9 - VBB) / RBD - ibd(y9 - y5)
-    
-    # Algebraic equation (KCL at node)
-    out[10] = (y10 - y6) / RGS + ibs(y8 - y10) + (y10 - y12) / RGD + ibd(y14 - y10)
-    
-    out[11] = CGS * dy11 - CGS * v2d + y11 / RGS + ids(2, y12 - y11, v2 - y11, y13, v2 - y12, y14 - y10)
-    out[12] = CGD * dy12 - CGD * v2d + (y12 - y10) / RGD - ids(2, y12 - y11, v2 - y11, y13, v2 - y12, y14 - y10)
-    out[13] = CBS * dy13 + (y13 - VBB) / RBS - ibs(y13)
-    out[14] = CBD * dy14 + (y14 - VBB) / RBD - ibd(y14 - y10)
-    
-    return nothing
-end
-
-# DAE initial conditions 
-u0_dae = [5.0, 5.0, VBB, VBB, 5.0, 3.62385, 5.0, VBB, VBB, 3.62385, 0.0, 3.62385, VBB, VBB]
-du0_dae = zeros(14)
+# Mass matrix (singular is fine!)
+dirMassMatrix = [
+    CGS 0   0   0   0   0   0   0   0   0   0   0   0   0
+    0   CGD 0   0   0   0   0   0   0   0   0   0   0   0
+    0   0   CBS 0   0   0   0   0   0   0   0   0   0   0
+    0   0   0   CBD 0   0   0   0   0   0   0   0   0   0
+    0   0   0   0   0   0   0   0   0   0   0   0   0   0
+    0   0   0   0   0   CGS 0   0   0   0   0   0   0   0
+    0   0   0   0   0   0   CGD 0   0   0   0   0   0   0
+    0   0   0   0   0   0   0   CBS 0   0   0   0   0   0
+    0   0   0   0   0   0   0   0   CBD 0   0   0   0   0
+    0   0   0   0   0   0   0   0   0   0   0   0   0   0
+    0   0   0   0   0   0   0   0   0   0   CGS 0   0   0
+    0   0   0   0   0   0   0   0   0   0   0   CGD 0   0
+    0   0   0   0   0   0   0   0   0   0   0   0   CBS 0
+    0   0   0   0   0   0   0   0   0   0   0   0   0   CBD
+]
 
-# Determine consistent initial conditions for DAE
-f_test = similar(u0_dae)
-nand_dae!(f_test, du0_dae, u0_dae, nothing, 0.0)
+# Initial conditions
+y0 = [5.0, 5.0, VBB, VBB, 5.0, 3.62385, 5.0, VBB, VBB, 3.62385, 0.0, 3.62385, VBB, VBB]
+tspan = (0.0, 80.0)
 
-# Create DAE problem
-dae_prob = DAEProblem(nand_dae!, du0_dae, u0_dae, tspan)
+# Mass matrix problem (original approach)
+mmf = ODEFunction(nand_rhs!, mass_matrix=dirMassMatrix)
+mmprob = ODEProblem(mmf, y0, tspan)
 
-# ModelingToolkit version via modelingtoolkitize on ODE problem
-mtk_sys = modelingtoolkitize(ode_prob)
+# ModelingToolkit version via modelingtoolkitize
+mtk_sys = modelingtoolkitize(mmprob)
 mtk_simplified = structural_simplify(mtk_sys)
-u0_mtk = [mtk_simplified.states[i] => ode_prob.u0[i] for i in 1:length(ode_prob.u0)]
-mtk_prob = ODEProblem(mtk_simplified, u0_mtk, tspan)
+u0_mtk = [mtk_simplified.states[i] => mmprob.u0[i] for i in 1:length(mmprob.u0)]
+mtkprob = ODEProblem(mtk_simplified, u0_mtk, mmprob.tspan)
+
+# DAEProblem version
+du = mmprob.f(mmprob.u0, mmprob.p, 0.0)
+du0 = D.(unknowns(mtk_simplified)) .=> du
+daeprob = DAEProblem(mtk_simplified, du0, [], tspan)
 
 # Generate reference solutions
-ref_sol = solve(ode_prob, Rodas5P(), abstol=1e-12, reltol=1e-12, tstops=0.0:5.0:80.0)
-dae_ref_sol = solve(dae_prob, IDA(), abstol=1e-10, reltol=1e-10, tstops=0.0:5.0:80.0)
-mtk_ref_sol = solve(mtk_prob, Rodas5P(), abstol=1e-12, reltol=1e-12, tstops=0.0:5.0:80.0)
+ref_sol = solve(mmprob, Rodas5P(), abstol=1e-12, reltol=1e-12, tstops=0.0:5.0:80.0)
+dae_ref_sol = solve(daeprob, IDA(), abstol=1e-10, reltol=1e-10, tstops=0.0:5.0:80.0)
+mtk_ref_sol = solve(mtkprob, Rodas5P(), abstol=1e-12, reltol=1e-12, tstops=0.0:5.0:80.0)
 
-probs = [ode_prob, dae_prob, mtk_prob]
-refs = [ref_sol, dae_ref_sol, mtk_ref_sol]
+probs = [mmprob, daeprob, mtkprob]
+refs = [ref_sol, ref_sol, mtk_ref_sol]
 ```
 
 ## Generate Reference Solution and Plot
 
 ```julia
-plot(ref_sol, title="NAND Gate Circuit - Node Potentials (ODE)", 
+plot(ref_sol, title="NAND Gate Circuit - Node Potentials (Mass Matrix)", 
      xlabel="Time", ylabel="Voltage (V)", legend=:outertopright)
 ```
 
@@ -416,10 +382,9 @@ plot(wp, title="NAND Gate DAE - Timeseries Errors (Low Tolerances)")
 ## Analysis of Key Nodes
 
 ```julia
-# For ODE formulation (12 variables): y1, y2, y3, y4, y6, y7, y8, y9, y11, y12, y13, y14
-# Map to original indices for interpretation
-key_nodes = [1, 2, 5, 6, 9, 10]  # Corresponds to y1, y2, y6, y7, y11, y12 in ODE form
-node_names = ["y1", "y2", "y6", "y7", "y11", "y12"]
+# Original 14-variable system: y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11, y12, y13, y14
+key_nodes = [1, 5, 6, 10, 11, 12]  # Representative nodes from different parts of the circuit
+node_names = ["Node 1", "Node 5", "Node 6", "Node 10", "Node 11", "Node 12"]
 
 p_nodes = plot()
 for (i, node) in enumerate(key_nodes)