Eigen for exponential matrix

Useful link
http://eigen.tuxfamily.org/dox/unsupported/group__MatrixFunctions__Module.html

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)

Any comments?