17 Numerical methods for linear DAE
Consider the DAE
\begin{aligned} x(t) &=\sin(t) \\ x'(t)+y(t) &= 0 \end{aligned}
with exact solution
\begin{aligned} x(t) &= \sin(t) \\ y(t) &= -\cos(t) \end{aligned}
Is a linear DAE \bm{F}\bm{x}'(t)+\bm{G}\bm{x}(t)=\bm{h}(t) where
\begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} x' \\ y' \end{pmatrix}+ \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix}= \begin{pmatrix} \sin(t) \\ 0 \end{pmatrix}
Using Kronecker/Weierstrass normal form with
\bm{U}=\bm{V}=\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \qquad\textrm{and}\qquad \bm{N}=\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}
the DAE is written as
\bm{N} \begin{pmatrix} y' \\ x' \end{pmatrix}+ \bm{I} \begin{pmatrix} y \\ x \end{pmatrix}= \begin{pmatrix} 0 \\ \sin(t) \end{pmatrix}
because \bm{N}\neq\bm{0} and \bm{N}^2=\bm{0} the DAE is of index 2. No initial conditions are required because block \bm{C} is absent.
17.1 Consistent initial values
Reconsider
\begin{aligned} x(t) &=\sin(t) \\ x'(t)+y(t) &= 0 \end{aligned}
I initial values x_0 and y_0 are consistent if
x(t) =\sin(t)
the hidden constraint is obtained by deriving the algebraic part
-y(t) = x'(t) = \cos(t)
and thus consistent initial condition must satisfy:
\begin{aligned} x_0 &=\sin(0) \\ y_0 &=-\cos(0) \end{aligned}
17.2 Explicit Euler
Consider the ODE \bm{x}'=\bm{f}(\bm{x},t), explicit Euler take the form
\bm{x}_{k+1} = \bm{x}_k + h \bm{f}(\bm{x}_k,t_k)
applied to the DAE become (t_k=kh)
\begin{aligned} 0 &= h(\sin(t_k) -x_k) \\ x_{k+1} &= x_k - h y_k \end{aligned}
Thus, the method cant be applied (y_{k+1} do not appear in the equations).
17.3 Implicit Euler
For \bm{x}'=\bm{f}(\bm{x},t), implicit Euler take the form
\bm{x}_{k+1} - h \bm{f}(\bm{x}_{k+1},t_{k+1}) = \bm{x}_k
Applied to the DAE result in (t_k=k\,h)
\begin{aligned} h(\sin(t_{k+1}) -x_{k+1}) &= 0 \\ x_{k+1} + h y_{k+1} &= x_k \end{aligned}
so that
\begin{aligned} x_{k+1} &= \sin(t_{k+1}) \\ y_{k+1} &= (x_k - \sin(t_{k+1})/h \\ &= (\sin(t_k) - \sin(t_{k+1}))/h =\color{red} -\cos(t_{k+1}) + \mathcal{O}(h) \end{aligned}
17.4 Numerical result
Crank-Nicolson
\frac{\bm{x}_{k+1}-\bm{x}_k}{h} = \frac{\bm{f}(\bm{x}_{k+1},t_{k+1})+\bm{f}(\bm{x}_k,t_k)}{2}
Applied to the DAE result in
\begin{aligned} \sin(t_{k+1}) -x_{k+1} &= x_k - \sin(t_k) \\ x_{k+1} - \frac{h}{2} y_{k+1} &= x_k + \frac{h}{2} y_{k} \end{aligned}
so that
\begin{aligned} x_{k+1} &= \sin(t_{k+1})+\sin(t_k)-x_k \\ y_{k+1} &= \frac{2}{h}(x_{k} - x_{k+1}) - y_k \end{aligned}
solving for x_k
x_k = \sin(t_k)+(-1)^k x_0
17.5 Damped stabilization (Baumgarte)
Re-consider the DAE
\begin{aligned} x(t) &= \sin(t) \\ x'(t)+y(t) &= 0 \end{aligned}
stabilize the constraint x(t)-\sin(t)=0 using Baumgarte
\begin{aligned} \tau (x(t)-\sin(t))' + (x(t)-\sin(t)) &= 0 \\ x'(t)+y(t) &= 0 \end{aligned}
The ODE \tau z'(t) + z(t) = 0 has solution
z(t) = z(0)\exp\bigg(-\frac{t}{\tau}\bigg) and \lim_{t\to\infty}z(t)=0
Using Crank-Nicolson
\begin{aligned} \frac{\tau}{h}(x_{k+1}-x_k)- \frac{\tau}{h}(\sin(t_{k+1})-\sin(t_k)) &+ \\[1em] \frac{1}{2}(x_{k+1}+x_k)- \frac{1}{2}(\sin(t_{k+1})+\sin(t_k)) &= 0 \\[1em] \frac{x_{k+1}-x_k}{h}+\frac{y_{k+1}+y_k}{2} &=& 0 \end{aligned}
solving for x_{k+1} and y_{k+1}
\begin{aligned} x_{k+1}&=\frac{2\,\tau-h}{2\,\tau+h}\big(x_k-\sin(t_k)\big)+\sin(t_{k+1}) \\[1em] y_{k+1}&=\frac{ \sin(t_k)(4\tau-2h)-\sin(t_{k+1})(4\tau+2h)+4hx_k}{h(2\tau+h)} -y_k \end{aligned}
17.6 Numerical result
Implicit Euler appears to work also with inconsistent data.
Crank–Nicolson does not work with inconsistent data.
Baumgarte stabilization works well, but there is no general recipe for determining stabilization parameters (\tau in this case).