# A question about Finite Element Method.

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.$$

• 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. – LutzL Jul 7 '17 at 16:39

## 2 Answers

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.