exponential integrators

Using exponential integrators is a way to eliminate stiffness from linear dynamical systems. Many pathsim blocks are pure linear odes such as the StateSpace blocks and its derivates, as well as the Differentiator and the PID.

They are more or less of the following form:

\dot{\vec{x}} = \mathbf{A} \vec{x} + \mathbf{B} \vec{u} 

Stiffness occurs when the eigenvalues of A are on vastly different scales, i.e. there are fast and slow time constants in the system.

Typpically the ODE is solved from one timestep to the next like this:

\vec{x}_{n+1}  = \vec{x}_{n} + \int_{t_n}^{t_{n+1}}   \mathbf{A} \vec{x}(\tau) + \mathbf{B} \vec{u}(\tau)  \ d\tau

But since the exact solution can be expressed for linea systems, the problem can be reduced to:

\vec{x}_{n+1}  = e^{\mathbf{A}h} \vec{x}_{n} + \int_{t_n}^{t_{n+1}}   e^{\mathbf{A} (t_{n+1} - \tau)} \mathbf{B} \vec{u}(\tau)  \ d\tau

It removes the state dependency from the integration. Intuitively it can be explained by superposition due to the linearity. The solution is an exponentially dacaying state added to the input convolved with the impulse response.

The resulting convolution integral can then be solved by some ode solver.

Implementing this as an option for linear ode blocks will help reduce stiffness induced by these kinds of blocks.

Start by locating the StateSpace blocks, their derivatives, the Differentiator, and the PID implementations that handle linear ODEs. Review how these blocks currently advance state, then define and test an optional exponential-integrator path that reduces stiffness for the stated linear form.