# How to implement Galerkin Method of Lines / FEM with black box integrators in scipy

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.

Edit:

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.

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.
• 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)$. – user3482876 May 6 '18 at 20:10
• 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))$. – Jonas May 6 '18 at 21:57
• 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. – user3482876 May 7 '18 at 18:16