B. Approximating PDEs by System of ODEs
Deriving the System of ODEs
A parabolic PDE can be approximated by a system of ODEs. This is known as the method of lines. In this section we use the simplest heat equation as an example to show how to derive the system of ODEs. Consider the initial boundary value problem
where the boundary conditions are set to zero to further simplify the argument below. First discretize the PDE in space. Let \(\{x_j\}_{j=0}^{m}\) be a uniform grid in \([-1, 1]\) with \(x_j = -1+2j/m\). We denote the solution at \(x_j\) by
They are functions of time. In particular, we have \begin{align*} U_0(t) &= u(t, x_0) = -1, \\ U_m(t) &= u(t, x_m) = 1. \end{align*} Next, denote by \(\vec U(t)\) the vector \begin{align*} \vec U(t) = \begin{pmatrix} U_1(t)\\ U_2(t)\\ U_3(t)\\ \vdots\\ U_{m-1}(t)\\ \end{pmatrix}. \end{align*} We will derive the system of ODEs for \(\vec U(t)\). Define the “shift-to-the-right” operator \(S\): \begin{align*} SU_j(t) = U_{j+1}(t) \end{align*} and the “shift-to-the-left” operator \(S^{-1}\): \begin{align*} SU_j(t) = U_{j-1}(t). \end{align*} Then define the finite difference operators: \begin{align*} \delta_x^+ &= \frac{S-I}{\Delta x},\\ \delta_x^- &= \frac{I-S^{-1}}{\Delta x},\\ \delta_x^+\delta_x^- &= \frac{S-2I+S^{-1}}{(\Delta x)^2}. \end{align*} We have \begin{align*} \delta_x^+\delta_x^- U_j(t) = \frac{U_{j+1}(t) -2 U_j(t) + U_{j-1}(t)}{(\Delta x)^2}\qquad\forall j=1, 2, \ldots, m-1. \end{align*} For a fixed \(t\), this is the familiar finite difference approximation of \(u_{xx}(t, x)\). We can write
But the PDE gives us \(u_{xx}(t, x_j) = u_t(t, x_j) = dU_j(t)/dt\). Thus we have
If we drop the error term and write in vector form, we obtain a system of ODEs that approximates the original heat equation: \begin{align*} \frac{d{\vec U(t)}}{dt} = \delta_x^+\delta_x^- \vec U(t). \end{align*}
Next we write the system in matrix form. The shift operators are \begin{align*} S\vec U(t) &= \begin{pmatrix} U_2(t)\\ U_3(t)\\ U_4(t)\\ \vdots\\ U_{m}(t)\\ \end{pmatrix} = \begin{pmatrix} 0 & 1 &&&\\ & 0 & 1 &&\\ && 0 & 1 &\\ &&& \ddots \\ &&&& 0\\ \end{pmatrix}\vec U(t) = \begin{pmatrix} 0 & 1 &&&\\ & 0 & 1 &&\\ && 0 & 1 &\\ &&& \ddots \\ &&&& 0\\ \end{pmatrix}\begin{pmatrix} U_1(t)\\ U_2(t)\\ U_3(t)\\ \vdots\\ U_{m-1}(t)\\ \end{pmatrix},\\ S^{-1}\vec U(t) &= \begin{pmatrix} U_0(t)\\ U_1(t)\\ U_2(t)\\ \vdots\\ U_{m-2}(t)\\ \end{pmatrix} = \begin{pmatrix} 0 &&&&\\ 1 & 0 &&&\\ & 1 & 0 &&\\ &&& \ddots \\ &&& 1 & 0\\ \end{pmatrix}\vec U(t) = \begin{pmatrix} 0 &&&&\\ 1 & 0 &&&\\ & 1 & 0 &&\\ &&& \ddots \\ &&& 1 & 0\\ \end{pmatrix} \begin{pmatrix} U_1(t)\\ U_2(t)\\ U_3(t)\\ \vdots\\ U_{m-1}(t)\\ \end{pmatrix}, \end{align*} where we have used the zero boundary conditions \(U_0(t) = U_m(t) = 0\). The finite difference operator \(\delta_x^+\delta_x^-\) in matrix form is \begin{align*} \delta_x^+\delta_x^-\vec U(t) = \frac{S-2I+S^{-1}}{(\Delta x)^2}\vec U(t) = \frac{1}{(\Delta x)^2} \begin{pmatrix} -2 & 1 &&&\\ 1 & -2 & 1 &&\\ & 1 & -2 & 1 &\\ &&& \ddots \\ &&& 1 & -2\\ \end{pmatrix}\vec U(t). \end{align*} In summary, the system of ODEs is
with the initial conditions \(U_j(0) = u_0(x_j)\) for all \(j=1, 2, \ldots, m-1\).
Why ODEs?
In most applications, we still numerically solve the system of ODEs obtained from discretizing a PDE in space only. In other words, we still need to discretize in time. It is natural to wonder why method of lines is needed?
There are benefits turning a PDE to ODEs. On the theory side, numerical schemes solving ODEs, e.g. Runge Kutta, and their behaviors have been well studied, including stability and convergence analysis, all directly applicable. On the engineering side, once the system of ODEs are derived, one can apply the existing numerical ODEs solver softwares. As an example, SciPy does not have PDE solvers at the time of this writing, only scipy.integrate.OdeSolver for system of ODEs, which can be utilized through method of lines for quick prototyping. Furthermore, space and time need not be discretized the same way. If we discretize the weak form of the PDE in space, we actually get a different system of ODEs. This leads to the spectral method and exponential convergence in space which finite difference method cannot achieve due to its local nature.