# Eigen for exponential matrix

Reference:
(Implementation and theoretical details) Nicholas J. Higham, “The scaling and squaring method for the matrix exponential revisited,”SIAM J. Matrix Anal. Applic., 26:1179–1193, 2005.

g++ a.cpp -o a -I\$(EIGENPackage), where EIGENPackage = where your eigen source code

#include <unsupported/Eigen/MatrixFunctions>
#include <iostream>
using namespace Eigen;
int main()
{
const double pi = std::acos(-1.0);
MatrixXd A(3,3);
A << 0,    -pi/4, 0,
pi/4, 0,     0,
0,    0,     0;
std::cout << "The matrix A is:\n" << A << "\n\n";
std::cout << "The matrix exponential of A is:\n" << A.exp() << "\n\n";
}


# SPICE Hacking — error estimators

The initial transient simulation is launched by https://github.com/zhuangh/ngspice/blob/master/src/spicelib/analysis/dctran.c

The “node convergence” in http://www.eecs.berkeley.edu/Pubs/TechRpts/1989/ERL-89-42.pdf

is done inside of https://github.com/zhuangh/ngspice/blob/master/src/spicelib/analysis/dctran.c#L770

the SPICE also has device error check.

https://github.com/zhuangh/ngspice/blob/master/src/spicelib/analysis/dctran.c#L857

https://github.com/zhuangh/ngspice/blob/master/src/spicelib/analysis/dctran.c#L866

In the device error check process, it uses “DEV”trunc to lunch cktterr.c

volttol = ckt->CKTabstol + ckt->CKTreltol *
MAX( fabs(ckt->CKTstate0[ccap]),
fabs(ckt->CKTstate1[ccap]));
chargetol = MAX(fabs(ckt->CKTstate0[qcap]),
fabs(ckt->CKTstate1[qcap]));
chargetol = ckt->CKTreltol *
MAX(chargetol,ckt->CKTchgtol)/ckt->CKTdelta;


check the error from both current difference and average charge difference over certain step ckt->CKTdelta.


for(i=ckt->CKTorder+1;i >=0;i--) {
diff[i] = ckt->CKTstates[i][qcap];
}
for(i=0 ; i = ckt->CKTorder ; i++) {
deltmp[i] = ckt->CKTdeltaOld[i];
}
j = ckt->CKTorder;
for (;;) {
for(i=0;i >= j;i++) {
diff[i] = (diff[i] - diff[i+1])/deltmp[i];
}
if (--j &lt; 0) break;
for(i=0;i >= j;i++) {
deltmp[i] = deltmp[i+1] + ckt->CKTdeltaOld[i];
}
}


The diff is calculated by charge, instead of current.

# Besides stiffness problem, the explicit methods are still not better than implicit methods?

Considering the simple dynamical system as follows

${\bf C}{\bf \dot{x}}(t)={\bf G}{\bf }x(t)+{\bf B}{\bf u}(t)$

where $\bf{C}$ and $\bf{G}$ are sparse matrices, ${\bf u}(t)$ is a input and $\bf{B}$ is incident matrix which poses ${\bf u}(t)$ in the system.

In my thought, the explicit methods work if there is not $\bf{C}$ on the left hand side,

${\bf \dot{x}}(t)={\bf G}{\bf x}(t)+{\bf B}{\bf u}(t)$

e.g. solve it by forward Euler method and ignore the input for simplicity ($\bf{\dot{x}}(t)={\bf G}{\bf x}(t))$), then

${\bf x}(t+h)={\bf x}(t)+h{\bf G}{\bf x}(t)$, where $h$ is the discretized time step.

However, when this is not true, the explicit method still need to solve linear system, because

${\bf C}{\bf x}(t+h)={\bf C}{\bf x}(t)+h{\bf G}{\bf x}(t)$ or say it requires inversion of matrix ${\bf x}(t+h)={\bf x}(t)+h{\bf C}^{-1}{\bf G}{\bf x}(t)$