Numerical ODE, DAE & Co-Simulation Solvers

Implementation: packages/simulate/src/solvers/ and packages/runtime/src/solvers/

---

1. Tsitouras 5(4) Adaptive Runge-Kutta (Tsit5)

Academic Citations

• Tsitouras, Ch. (2011). "Runge-Kutta pairs of order 5 (4) satisfying only the first column simplifying assumption."
  Computers & Mathematics with Applications, 62(2), pp. 770–775.
  DOI: 10.1016/j.camwa.2011.06.002
• Rackauckas, C., & Nie, Q. (2017). "DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia."
  Journal of Open Research Software, 5(1), 15. DOI: 10.5334/jors.151

ModelScript Rationale & Modifications

For non-stiff systems (e.g. mechanical multi-body dynamics, kinematics, flight control), high-order adaptive explicit Runge-Kutta pairs offer peak throughput without Newton-Raphson linear algebra costs. Tsit5 demonstrates significantly lower error coefficients and wider stability regions than classical Dormand-Prince (DOPRI5).

• Modifications: Employs the First-Same-As-Last (FSAL) optimization eliminating stage 1 evaluation in successive steps; embeds a 5th-order continuous dense output polynomial for sub-interval event zero-crossing localization.

---

2. TR-BDF2 Composite Stiff Solver

Academic Citations

• Bank, R. E., Coughran, W. M., et al. (1985). "Transient simulation of silicon devices and circuits."
  IEEE Transactions on CAD of Integrated Circuits, 4(4), pp. 436–446.
  DOI: 10.1109/TCAD.1985.1270142
• Hosea, M. E., & Shampine, L. F. (1996). "Analysis and implementation of TR-BDF2."
  Applied Numerical Mathematics, 20(1-2), pp. 21–37.
  DOI: 10.1016/0168-9274(95)00115-8

ModelScript Rationale & Modifications

Cyber-physical models featuring frequent hybrid discrete events (switches, impacts, hysteresis) severely penalize multistep BDF solvers (which must drop to order 1 after every event). TR-BDF2 is an L-stable one-step composite method pairing a trapezoidal half-step with a BDF2 completion step:

• Modifications: Shares the diagonal Jacobian coefficient $d=1-1/\sqrt{2}$ across both internal stages, requiring only one sparse LU factorization per time step; connects directly with ModelScript's in-WASM sparse matrix kernel (sparse_lu.ts).

---

3. RODAS-4P Linearly Implicit Rosenbrock Solver

Academic Citations

• Steinebach, G. (1995). "Order-reduction of ROW-methods for DAEs and method of lines applications."
  TH Darmstadt Preprint 1742.
• Hairer, E., & Wanner, G. (1996). Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems.
  Springer-Verlag. Section IV.8.

ModelScript Rationale & Modifications

Solves stiff index-1 DAEs and method-of-lines continuum discretizations without iterative non-linear Newton solves. It solves only a sequence of linear systems sharing the same system matrix $I-\gamma hJ$.

• Modifications: Optimized for in-WASM sparse Jacobian caching; utilizes automatic differentiation coloring (coloring.ts) for fast Jacobian updates.

---

4. Itô Stochastic Differential Equations (SDE) Solver

Academic Citations

• Rößler, A. (2010). "Runge-Kutta methods for Itô stochastic differential equations with additive noise."
  SIAM Journal on Numerical Analysis, 48(3), pp. 922–952. DOI: 10.1137/09076636X
• Kloeden, P. E., & Platen, E. (1992). Numerical Solution of Stochastic Differential Equations.
  Springer. DOI: 10.1007/978-3-662-12616-5

ModelScript Rationale & Modifications

Physical systems operate under thermal noise, turbulence, and sensor fluctuations:

$$dx=f(t,x)dt+g(t,x)d{W}_t$$

• Modifications: Employs an in-WASM Xoshiro256pp PRNG with Box-Muller transformation; provides both fixed-step Euler-Maruyama (strong order 0.5) and adaptive Rößler SRIW1 (strong order 1.5) integrators for Monte Carlo uncertainty rollouts.

---

5. Co-Simulation Master Orchestrator

Academic Citations

• Bastian, J., et al. (2011). "Master for Co-Simulation Using FMI."
  Proceedings of the 8th International Modelica Conference, pp. 115–120. DOI: 10.3384/ecp11063115
• Kübler, R., & Schiehlen, W. (2000). "Two methods of simulator coupling."
  Mathematical and Computer Modelling of Dynamical Systems, 6(2), pp. 93–113. DOI: 10.1076/1387-3954(200006)6:2;1-M;FT093

ModelScript Rationale & Modifications

Coordinates multi-rate, multi-vendor FMI Functional Mock-up Units (FMUs) and SSP topologies. Implements Richardson extrapolation error estimation to dynamically adapt the inter-simulator communication step size $H_{\text{comm}}$, preventing divergence while maximizing simulation speed.

---

6. Dommel Companion-Circuit EMT & Bergeron Traveling-Wave Solver

Academic Citations

• Dommel, H. W. (1969). "Digital Computer Solution of Electromagnetic Transients in Single- and Multiphase Networks."
  IEEE Transactions on Power Apparatus and Systems, PAS-88(4), pp. 388–399.
  DOI: 10.1109/TPAS.1969.292459
• Marti, J. R. (1982). "Accurate Modelling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Simulations."
  IEEE Transactions on Power Apparatus and Systems, PAS-101(1), pp. 147–157.
  DOI: 10.1109/TPAS.1982.317332

ModelScript Rationale & Modifications

Implemented in EmtSolver (packages/simulate/src/solvers/emt/emt-solver.ts), this solver targets microsecond-scale electromagnetic transients, power electronics switching (IGBTs/breakers), and distributed transmission lines imported from @modelscript/pscad:

• Trapezoidal Companion Models & Constant Nodal Matrix: Inductors ($G_L=\frac{\mathrm{\Delta }t}{2L}$) and capacitors ($G_C=\frac{2C}{\mathrm{\Delta }t}$) are discretized into Norton companion equivalents $\mathbf{G}\mathbf{v}(t)={\mathbf{i}}_{\text{hist}}(t-\mathrm{\Delta }t)+{\mathbf{i}}_{\text{src}}(t)$, allowing $\mathbf{G}$ to remain constant between switching events.
• Critical Damping Adjustment (CDA): Automatically transitions from Trapezoidal to Backward Euler integration ($G_L=\frac{\mathrm{\Delta }t}{L},G_C=\frac{C}{\mathrm{\Delta }t}$) for two half-steps ($\mathrm{\Delta }t/2$) immediately after any breaker or semiconductor switching discontinuity to suppress artificial $dv/dt$ numerical oscillations.
• Bergeron Traveling-Wave Line Decoupling & Norton FMU Interface: Distributed transmission lines with propagation delay $\tau =\ell /v\ge \mathrm{\Delta }t$ decouple the nodal admittance matrix into independent sub-grids using interpolated history ring buffers, and expose getNortonEquivalent(portNode) ($I_{\text{eq}},G_{\text{eq}}$) for zero-delay FMI Co-Simulation coupling.

---

7. High-Order Taylor Model Flowpipe Integration

Academic Citations

• Makino, K., & Berz, M. (2003). "Taylor Models and Other Validated Functional Inclusion Methods."
  International Journal of Pure and Applied Mathematics, 4(4), pp. 379–456.
• Chen, X., Ábrahám, E., & Sankaranarayanan, S. (2013). "Flow*: An Analyzer for Non-Linear Hybrid Systems."
  Computer Aided Verification (CAV 2013), LNCS vol. 8044, pp. 258–263.
  DOI: 10.1007/978-3-642-39799-8_18

ModelScript Rationale & Modifications

Implemented in TaylorModel (packages/runtime/src/analysis/wasm_taylor_model.ts), validated Taylor Model arithmetic bounds all trajectories of a non-linear ODE/DAE system subject to uncertain initial sets $x(0)\in X_0$ and parameter intervals without wrapping-effect blowup:

• Polynomial + Interval Remainder Representation: Represents each state reachable tube over normalized domain $\mathbf{z}\in [-1,1{]}^d$ as a multivariate polynomial $P_k(\mathbf{z})$ of order $k$ plus a sound interval remainder $I=[\underset{―}{r},\bar{r}]$, where truncated higher-order terms ($>k$) are conservatively swept into $I$.
• Picard-Lindelöf Validated Operator: Evaluates the contractive Picard integral operator $\mathcal{P}(X)(t)=X_0+{\int }_0^{t}f(s,X(s))ds$ directly on Taylor models to construct guaranteed continuous-time flowpipes consumed by ContinuousSafetyOracle and SysML v2 hybrid state-machine reachability (hybrid-flowpipe-bridge.ts).

---

Upstream & Downstream Pipeline Connections

• Upstream Inputs:
  • DAE blocks partitioned and torn by BltEngine and TornBlock.
  • Sparse Jacobians computed via Dual Numbers & Graph Coloring.
• Downstream Consumers:
  • Powers msx simulate, msx optimize, and the Web IDE charting canvas.
  • Feeds state trajectories and validated flowpipes into Signal Temporal Logic (STL) Monitor & Theory Coordinator and formal verification oracles.