Can we solve a matrix differential equation $\textbf{X}'(t)=\textbf{A}\textbf{X}(t)+\textbf{B}(t)$ by Finite Element Method?

I will be very happy if you give some special example or suggest some books about the method for solving matrix differential eauations.

In here, \begin{equation} \begin{aligned} &\textbf{X}(t)=\left( \begin{array}{cccc} x_{1}(t) \\ x_{2}(t) \\ x_{3}(t)\\ x_{4}(t) \end{array} \right), \textbf{A} = \left( \begin{array}{cccc} a_{11} & a_{12} &a_{13} & a_{14}\\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32}& a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{array} \right),\\ &\textbf{B}(t)=\left( \begin{array}{cccc} b_{1}(t)\\ b_{2}(t)\\ b_{3}(t)\\ b_{4}(t) \end{array} \right). \end{aligned} \end{equation} initial conditions: $$x_1(0)=m_1$$ $$x_2(0)=m_2$$ $$x_3(0)=m_3$$ $$x_4(0)=m_4.$$

  • 2
    $\begingroup$ Finite elements serve to discretize a solid body or some area. What is the continuum that you can discretize here? Time? Then see the trapezoidal method or other symmetric implicit Runge-Kutta methods. $\endgroup$ – LutzL Jul 7 '17 at 16:39

The differential system ${\bf X}' = {\bf A} {\bf X} + {\bf B}$ where $t \mapsto {\bf X}(t)$ is a vector-valued function of size $n=4$, ${\bf A}$ is a $n$-by-$n$ matrix and $t \mapsto {\bf B}(t)$ is a vector-valued function of size $n$, can be solved by hand in many particular cases. Indeed, this system is an order-1 linear system of ordinary differential equations (ODE). Solutions write as $$ {\bf X}(t) = \exp(t{\bf A})\, {\bf X}(0) + \int_{0}^t \exp((t-\tau) {\bf A}) \, {\bf B}(\tau)\, \mathrm{d}\tau \, , $$ where $\exp(s{\bf A})$ is a matrix exponential. Therefore, if the matrix exponential and the integral can be computed, then the solution is given above.

(1) FEM. To implement FEM, the first step is to write a weak formulation of the ODE. We write the scalar product of the ODE with the test vector field ${\bf Y}(t)$, and we integrate over $t\in \Theta$. Hence, using integration by parts, $$ \int_\Theta ({\bf Y}' + {\bf A}^\top{\bf Y})^\top\, {\bf X} = -\int_\Theta {\bf Y}^\top {\bf B}, $$ for all ${\bf Y}$ which trace equals zero on the boundaries of $\Theta$. One notes that the above weak formulation writes as $a({\bf X},{\bf Y}) = f({\bf Y})$, where $a(\cdot,\cdot)$ is a bilinear form over $H^1(\Theta)$ and $f(\cdot )$ is a linear form. It remains to discretize the above equations by considering that ${\bf X}$, ${\bf Y}$ belong to finite-dimension spaces, e.g. by using appropriate bases of Lagrange polynomials of the variable $t$.

(2) FD. The system ${\bf X}' = {\bf A} {\bf X} + {\bf B}$ can be integrated numerically by finite-difference methods. The simplest method is the first-order accurate Euler method. Here, the forward Euler method writes $$ \left\lbrace \begin{aligned} &{\bf X}_{n+1} = {\bf X}_{n} + \Delta t \left({\bf A}{\bf X}_{n} + {\bf B}(t_n)\right) , \\ &{\bf X}_{0} = {\bf X}(0) \, , \end{aligned} \right. $$ where ${\bf X}_{n} \simeq {\bf X}(t_{n})$, and $\Delta t = t_{n+1} - t_n$ denotes the time step. The method is stable provided that $ \Delta t \leq {2}/{\varrho({\bf A})} $, where $\varrho({\bf A})$ is the spectral radius of the matrix $\bf A$. The backward Euler method $$ \left\lbrace \begin{aligned} &{\bf X}_{n+1} = \left({\bf I} - \Delta t {\bf A}\right)^{-1}\, \left({\bf X}_{n} + \Delta t\, {\bf B}(t_{n+1})\right) , \\ &{\bf X}_{0} = {\bf X}(0) \, , \end{aligned} \right. $$ where ${\bf I}$ is the identity matrix, is unconditionally stable. Thus, it bypasses stability issues, but requires the computation of a matrix inverse. Many more accurate numerical methods exist, such as Runge-Kutta methods. If the system is stiff, then the time step $\Delta t$ must be chosen very small to get good accuracy, which is costly from a computational point of view. Hence, dedicated adaptive methods may be preferred. For instance, the reader is referred to the book Numerical Recipes in C: The Art of Scientific Computing for details.


This is perfectly possible. First we need a weak formulation: take smooth $v$ then $$ \int_0^T v^T(x' - Ax) dt = \int_0^T v^T Bdt $$ is a nice weak formulation. It remains to choose the spaces $X_h$ and $V_h$ to take $x$ and $v$ from.

If $X_h=V_h$=piecewise polynomials, then the resulting system completely couples over all time steps, which is not nice.

A sensible choice is to take one of $X_h$ and $V_h$ piecewise polynomial and continuous, the other space piecewise polynomial but discontinuous. You end up with a Petrov-Galerkin or discontinuous Galerkin scheme. Due to the discontinuity in one of the variables, the time steps are decoupled.

Choosing piecewise linear and continuous functions for $x$, but discontinuous and piecewise constant functions for $v$ results in a variant of the Crank-Nicolson scheme.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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