# time discretization via FEM to solve ODEs

It seems that FEMs are never seen to implement time discretization. Is time discretization via FEMs to solve ODEs ever possible?

• I don't think there is anything stopping anyone from using FEM for time discretization, but finite different does the job well for 1-D problems. I am not a mathematician but, I think that the mathematical foundations for time discritaziation are also quite mature. There are space-time FEM techniques where both time and space are discritized with FEM. Adding an extra dimension for the problem never justified the use of something so expensive so they never became popular. Maybe in the future. Dec 19, 2022 at 8:56
• In fact, it was. Search for space-time finite element discretisation methods. Quite an active research area. Dec 19, 2022 at 12:58
• People do derive discontinuous Galerkin schemes for time discretization. Both the explicit and implicit Euler methods arise from this family. Dec 21, 2022 at 3:00

Let's consider the Newton's law of cooling $$u' = k(u_1 - u)$$ with the initial condition $$u(0)=u_0$$. Following the standard process of deriving a weak formulation, we multiply by a test function and integrate over the domain $$[0, T]$$ where $$T$$ is the final time: $$\int_0^T u'v\,dx + k \int_0^T uv\,dx = k \int_0^T u_1 v\,dx.$$ To me it looks like a completely reasonable weak formulation for finite elements. Let's try implementing it using piecewise-linear elements and comparing to the analytical solution $$u(t) = u_1 + (u_0 - u_1)e^{-kt}.$$ In the following snippet we have $$u_0 = 10$$, $$u_1=5$$, $$T=10$$.

# Dependencies: pip install scikit-fem[all]==8.0.0
from skfem import *
from skfem.helpers import *
import numpy as np
import matplotlib.pyplot as plt

m = MeshLine(np.linspace(0., 10., 10))
basis = Basis(m, ElementLineP1())
k = .5  # conductivity
u1 = 5.  # target temperature

@BilinearForm
def bilinf(u, v, w):
return d(u)[0] * v + k * u * v

@LinearForm
def linf(v, w):
return k * u1 * v

u = basis.zeros()
u[0] = 10.  # initial temperature

A = bilinf.assemble(basis)
b = linf.assemble(basis)

u = solve(*condense(A, b, x=u, D=np.array([0])))

ax = basis.plot(u, nrefs=0)
ax.plot(np.linspace(0, 10, 50), u1 + (10. - u1) * np.exp(-k * np.linspace(0, 10, 50)))
plt.legend(['discrete', 'analytical'])


Running the above snippet leads to the following graph:

So it is definitely possible!

• Assuming exact integration, this is equivalent to a central difference approximation in the interior of the domain and a one-sided difference for the boundary DOFs. Dec 20, 2022 at 8:30