Suppose I have some time dependent PDE, which can be written in the strong form as

$$ \frac{\partial u}{\partial t} + \mathcal{L}(u) = f $$

Where $\mathcal{L}$ is some differential operator. If I write the weak form of the problem without discretizing the time term, I end up with the standard weak form problem of finding $u_h \in X_h$ such that

$$ m\left(\frac{\partial u_h}{\partial t},\, v\right) + a(u_h,\,v) = f(v), \qquad \forall v \in X_h $$

At this point, I usually discretize in time in the following way (using backward-euler here) as

$$ m\left(\frac{u_h^{k+1} - u_h^k}{\Delta t},\, v\right) + a(u_h^{k+1 },\,v) = f(v), \qquad \forall v \in X_h $$

at which point, I know how to solve the problem by putting the $m(u_h^k, v)$ term to the right hand side, choosing a test, trial space, discretizing and solving.

But what if I wanted to start from the second equation and structure my code to use a black box integrator? Is there a standard way of doing this? I know that the right hand side terms can become much more complicated with higher-order time discretizations.


It seems that the answer is to first discretize in space, and express the discrete system as a matrix system $B\dot{U_h} + AU_h = F$. However, if the black box time integrator library does not support DAE systems (allowing user to input $B$) matrix, you're out of luck (this seems to be the case with scipy.integrate.ode). If $B$ can be inverted cheaply, then the black box integrator can handle the equivalent problem $\dot{U_h} = -B^{-1}(AU_h - F)$. Alternatively, mass lumping can be used, but this seems ridiculous if you were interested in implicit time-marching schemes to begin with.


1 Answer 1


Do you have a particular black box ode integrator in mind to solve? $$ m\left(\frac{\partial u_h}{\partial t},\, v\right) + a(u_h,\,v) - f(v) = 0, \qquad \forall v \in X_h $$ and initial conditions, those I will skip also in the rest of these notes. In general, this is just a differential-algebraic system of equations (DAE)$$ F(t, u, u_{t}) = 0. $$ that could also be solved with a solver for DAEs, such as ode15i in MATLAB.

  • $\begingroup$ I was planning on using the scipy.integrate.ode libraries for some type of implicit time-stepping. While I agree with everything you wrote, it still isn't clear to me how to do it; the time derivative term is being hit with a mass operator term. It seems what you've written is $g(y'(t)) = f(t, y)$ instead of the form $y'(t) = f(y,t)$. $\endgroup$ May 6, 2018 at 20:10
  • $\begingroup$ Yes, but there are solvers to deal with the $g(y'(t)) = f(t,y)$ type ODE (actually called DAE), also implicit ones. On the other hand, you could multiply the whole system of equations with the inverse mass operator - whether that is feasible or not depends on the size of the finite element space. But then you get $g(y'(t)) = f(t,y) \Leftrightarrow y'(t) = g^{-1}(f(t,y))$. $\endgroup$
    – Jonas
    May 6, 2018 at 21:57
  • $\begingroup$ I see your point. It seems that my issue stems from the fact that scipy.integrate.ode doesn't seem to support this type of DAE feature. $\endgroup$ May 7, 2018 at 18:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.