Struct JacobianData 

pub struct JacobianData<'a> {
    pub jacobian: &'a mut DMatrix<f64>,
    pub coefficient_matrix: &'a DMatrix<f64>,
    pub unknowns: &'a DVector<f64>,
    pub excitation: &'a DVector<f64>,
    pub edge_to_system: &'a DMatrix<f64>,
    pub edge_types: &'a [Type],
    pub edge_resistances: &'a DVector<f64>,
    pub edge_resistances_prev: &'a DVector<f64>,
    pub edge_currents: &'a DVector<f64>,
    pub edge_currents_prev: &'a DVector<f64>,
    pub edge_voltages: &'a DVector<f64>,
    pub edge_voltages_prev: &'a DVector<f64>,
}

A struct containing information about the nonlinear equation system (both its general structure and the concrete values of its component). It is provided as an input to a user-supplied Jacobian calculation function,

The network solver tries to solve the general matrix equation A * x = b, with A being the coefficient matrix, x being the unknowns (the mesh currents for MeshAnalysis or the node potential for NodalAnalysis and b being the excitations (current and / or voltage sources).

This is equivalent to searching the root of A * x - b = 0, hence the Newton-Raphson root finding algorithm can be applied. This requires calculating the Jacobian, where every equation of the system is derived by every value of x (see https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant). By default, the network solver uses a finite differences approximation of the Jacobian, but providing an analytical calculation method function return better results with regards to the number of iterations needed.

The following example shows how to write a custom Jacobian calculation function.

use approx; // Needed for result assertion
use network_analysis::*;

// Problem definition
// =====================================================================

/*
This creates the following network with two voltage sources at 0 and 6
 ┌─[1]─┬─[2]─┐
[0]   [6]   [3]
 └─[5]─┴─[4]─┘
*/
let mut edges: Vec<EdgeListEdge> = Vec::new();
edges.push(EdgeListEdge::new(vec![5], vec![1], Type::Voltage));
edges.push(EdgeListEdge::new(vec![0], vec![2, 6], Type::Resistance));
edges.push(EdgeListEdge::new(vec![1, 6], vec![3], Type::Resistance));
edges.push(EdgeListEdge::new(vec![2], vec![4], Type::Resistance));
edges.push(EdgeListEdge::new(vec![3], vec![5, 6], Type::Resistance));
edges.push(EdgeListEdge::new(vec![4, 6], vec![0], Type::Resistance));
edges.push(EdgeListEdge::new(vec![1, 2], vec![4, 5], Type::Voltage));
let network = Network::from_edge_list_edges(&edges).expect("valid network");
let mut mesh_analysis = MeshAnalysis::new(&network);

/*
All resistances are the square of the current going through them + 1
The resulting equation system is therefore nonlinear and requires iterations to be solved.
*/
let resistances = EdgeValueInputs::Function(&|mut args: FunctionArgs<'_>| {
    for (idx, val) in args.edge_value_and_type.iter_mut() {
         *val = 1.0 + args.edge_currents[idx].powi(2);
    }
});

let voltages = EdgeValueInputs::Slice(&[2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]);
let config = SolverConfig::default();

/**
The derivative of `b` with respect to `x` is zero (voltage sources are constant).
The derivative `d(A * x)/dx` can be rewritten as `dA/dx * x + A * dx/dx` (product rule).
`dx/dx` is the unity vector, hence the equation simplifies to `dA/dx * x + A`.

To calculate `dA/dx`, one needs to understand that the elements (sum of resistances in this case) of `A` are a
function of the individual edge currents, not of `x` (`x` is a sum of multiple edge currents).
Therefore, the following formula results for deriving any matrix element by x:
`dA[row, col]/dx[row] = d(edge_current[0])/dx[row] * d(edge_resistance[0])/d(edge_current[0]) + d(edge_current[1])/dx[row] * d(edge_resistance[1])/d(edge_current[1]) + ...`
The value of `d(edge_current[edge])/dx[row]` is equal to `-data.edge_to_system[(row, edge)]` (which describes the coupling between meshes and edges).
*/
fn calc_jacobian(data: JacobianData<'_>) {
    let ncols = data.jacobian.ncols();
    let nrows = data.jacobian.nrows();
    let nedges = data.edge_types.len();

    // DMatrix is column-major
    for col in 0..ncols {
        for row in 0..nrows {
            let mut sum = data.coefficient_matrix[(row, col)]; // this is A[row, col] from the function docstring
            for (edge, edge_type) in data.edge_types.iter().enumerate() {
                // Filter out excitation edges
                if Type::Resistance == *edge_type {
                    // Read out the negative coupling from the system matrix
                    let coupling = -data.edge_to_system[(row, edge)];

                    // Derivative of `edge_currents[idx].powi(2)` is:
                    let derivative = coupling * 2.0 * data.edge_currents[edge];

                    // Add to entry
                    sum += derivative;
                }
            }
            data.jacobian[(row, col)] = sum;
        }
    }
}

// Solving without finite difference approximation
let sol = mesh_analysis.solve(resistances, EdgeValueInputs::None, voltages, None, None, None, &config).expect("can be solved");
let edge_currents_fd = sol.currents().to_vec();
assert_eq!(sol.iter_count(), 119);

// Solving with analytical Jacobian
let sol = mesh_analysis.solve(resistances, EdgeValueInputs::None, voltages, None, None, Some(&mut calc_jacobian), &config).expect("can be solved");
let edge_currents_an = sol.currents().to_vec();
assert_eq!(sol.iter_count(), 34);

for (curr_fd, curr_an) in edge_currents_fd.iter().cloned().zip(edge_currents_an.iter().cloned()) {
    approx::assert_abs_diff_eq!(curr_fd, curr_an, epsilon = 1e-6);
}

Fields

jacobian: &'a mut DMatrix<f64>

Jacobian buffer. The user supplied-function is expected to populate this matrix.

coefficient_matrix: &'a DMatrix<f64>

The A in A * x = b

unknowns: &'a DVector<f64>

The x in A * x = b

excitation: &'a DVector<f64>

The b in A * x = b

edge_to_system: &'a DMatrix<f64>

This is the system structure matrix (MeshAnalysis::edge_to_mesh or NodalAnalysis::edge_to_node). It contains all information regarding the structure of the equation system, please see the documentation of the respective methods for more.

edge_types: &'a [Type]

Edge types of the network.

edge_resistances: &'a DVector<f64>

Edge resistance calculated in the current iteration

edge_resistances_prev: &'a DVector<f64>

Edge resistance calculated in the last iteration.

edge_currents: &'a DVector<f64>

Current through the edges calculated in the current iteration

edge_currents_prev: &'a DVector<f64>

Current through the edges calculated in the last iteration

edge_voltages: &'a DVector<f64>

Voltage drop over the edges calculated in the current iteration.

edge_voltages_prev: &'a DVector<f64>

Voltage drop over the edges calculated in the last iteration.