Solving differential algebraic equations with SDC¶
This project contains the sweepers, hooks, example problems, plotting and simulation scripts for a Master’s thesis investigating the usage of SDC methods to solve differential algebraic equations.
To run the scripts contained in this project a standard installation of pySDC should suffice.
Project overview¶
- misc
HookClass_DAE.py
Simple hook classes to read out the approximate solution and error after each time step.
- plotting
linear_plot.py
Reads a previously generated data file in .npy format and generates a plot on linear axis of the specified parameters.loglog_plot.py
Reads a previously generated data file in .npy format and generates a plot on logarithmic axis of the specified parameters. Commonly used to generate convergence plots.semilogy_plot.py
Reads a previously generated data file in .npy format and generates a plot on logarithmic y-axis and linear x-axis.
- problems
simple_DAE.py
A number of simple examples of differential algebraic equations are implemented here: a linear index-2 system with known analytical solution, the 2D pendulum as an index-3 system and a very simple fully implicit index-2 system that is not solvable by most numerical methods for certain values of a parameter.synchronous_machine.py
Synchronous machine model attached to an infinite bus undergoing torque disturbance test. Results in an index-1 system.transistor_amplifier.py
A two transistor amplifier model that results in an index-1 differential algebraic system. A nice example of a system resulting from a common real world situation.
- run
run_convergence_test.py
Script to generate convergence data of applying SDC to the simple linear index-2 differential algebraic system mentioned above.run_iteration_test.py
Script to generate data describing behaviour of error and residual of applying SDC to the simple linear index-2 differential algebraic system mentioned above.fully_implicit_dae_playground.py
Testing arena for the fully implicit sweeper.synchronous_machine_playground.py
Testing arena for the synchronous machine model.
- sweepers
fully_implicit_DAE.py
Sweeper that accepts a fully implicit formulation of a system of differential equations and applies to it a modified version of spectral deferred correction
Theoretical details¶
A fully implicit representation of a system of differential equations takes the form
A special case of such an implicit differential equation arises when the Jacobian \(\partial_{u'}F\) is singular. This implies that the derivative of some of the components of \(u(t)\) do not appear in the system of equations. The system is thus denoted a differential algebraic system of equations.
Since the derivative \(u'(t)\) cannot be isolated the Picard formulation used in ordinary SDC cannot be used here. Instead the derivative, henceforth denoted by \(U(t)\), is cast as the new unknown solution and the implicit system of equations is written as
Based on this equation and an initial approximate solution \(tilde{U}\), the following error equation is formed
This results directly in
from which the following time marching discretisation becomes obvious
The spectral integration matrix \(\mathbf{Q}\) is used to approximate the integral of the current approximation \(\tilde{U}\) and a low order approximation, in this case implicit Euler, is used for the unknown error \(\delta(t)\). Combining each step in the time marching scheme into a vector results in the following matrix formulation
with the integration matrix of the implicit Euler method
Finally, the iterative nature of the method is made clear by considering that the approximate solution can be updated repeatedly with a \(\tilde{\mathbf{\delta}}\) that is recalculated after each iteration and using the previously updated solution as the initial condition for the next iteration. In this way, reformulation of the previous equation as
results in the following iterative scheme
In practice each iteration is carried out line by line and the resulting implicit equation for \(U_{m+1}^{k+1}\) is solved using the familiar scipy.optimize.root()
function.