Runge-Kutta method for PDE

Question

I consider certain partial differential equation (PDE). For example, let it be heat equation

$$u_t = u_{xx}$$

I want to apply numerical Runge-Kutta method for solving it. As a first step I approximate $u_{xx}$ with difference scheme of several order. Let it be

$$ u_{xx} \approx \frac{u_{i+2} - 2 u_{i+1} + u_i}{h^2} $$ Then I take $t$ as a continious value and want to solve ODE with Runge-Kutta method:

$$u_t = \frac{u_{i+2} - 2 u_{i+1} + u_i}{h^2}$$

Here we may know only $u_0$ and $u_{N}$ values (boundary conditions).

This way I should now replace $u_t$ with difference scheme, for example

$$u_t \approx \frac{u_{\alpha + 1} - u_\alpha}{\tau}$$ (greek indices stand for time-lattice points; all values are taken in given x-lattice point $x_j$). And now I may write the second equation:

$$u_{xx} = \frac{u_{\alpha + 1} - u_\alpha}{\tau} $$ where $x$ is treated as continious value.

So I have a system of ODEs:

$$u_t = \frac{u_{i+2} - 2 u_{i+1} + u_i}{h^2} $$ $$u_{xx} = \frac{u_{\alpha + 1} - u_\alpha}{\tau}$$ with some initial-boundary conditions.

And now the problem is finding a way to solve it numerically (Runge-Kutta method). How should I obtain values that are necessary for solving (e.g. $u_{\alpha+1}$ or $u_{i+2}$)? It's not clear for me because for finding $u_{\alpha+1}$ we should know $u_{i+2}$ and equations are engaging.

Answer 1

Assuming you're using method of lines.

Let the original initial-boundary problem be $$ u_t = u_{xx}\\ u(0, x) = f(x)\\ u(t, 0) = a(t)\\ u(t, 1) = b(t). $$ Introduce a set of points $x_j = jh,\; j = 0,1,\dots, N,\;Nh = 1$. Let $u_j(t) = u(t, x_j)$. Note that $u_0(t) = a(t),\; u_N(t) = b(t)$ are already known. Unknown are $u_j(t),\; j = 1, 2, \dots, N-1$. Then $u_{xx}(t, x_j)$ can be approximated as $$ u_{xx}(t, x_j) \approx \frac{u_{j-1}(t) - 2u_j(t) + u_{j+1}(t)}{h^2}. $$ Plugging that into PDE gives $$ u'_j(t) = \frac{u_{j-1}(t) - 2u_j(t) + u_{j+1}(t)}{h^2} $$ a system of $N-1$ ODEs with initial conditions $u_j(0) = f(x_j)$. That could be solved using any RK method (provided that method is stable). For explicit Euler method that would be $$ \frac{u_j(t_{n+1}) - u_j(t_n)}{\tau} = \frac{u_{j-1}(t_n) - 2u_j(t_n) + u_{j+1}(t_n)}{h^2}, \; j = 1, 2, \dots, N-1\\ u_j(0) = f(t_j)\\ u_0(t) = a(t), \; u_N(t) = b(t). $$ This well known explicit scheme is stable when $\frac{\tau}{h^2} \leqslant \frac{1}{2}$.

Edit. For those who ask how to use this method with higher order RK methods. Let's take for example an RK2 method (explicit midpoint) with the following Butcher's tableau: $$ \begin{array}{c|cc} 0 & 0 & 0\\ 1/2 & 1/2 & 0\\ \hline & 0 & 1 \end{array} $$ Applied to ODE in form $\mathbf u' = \mathbf F(t, \mathbf u)$ this method expands to $$ \mathbf r = \mathbf F(t_n, \mathbf u_n)\\ \mathbf s = \mathbf F\left(t_n + \frac{\tau}{2}, \mathbf u_n + \frac{\tau}{2} \mathbf r\right)\\ \frac{\mathbf u_{n+1} - \mathbf u_n}{\tau} = \mathbf s $$ I've used $\mathbf r$ and $\mathbf s$ instead of common $\mathbf k_{1,2}$ to reduce the number of indices involved. Here $\mathbf r$ and $\mathbf s$ are intermediate values that depend solely on values of $u_j$ at $t = t_n$.

For our case the right hand side of the ODE is $$ F_j(t, u_1, u_2, \dots, u_{N-1}) = \frac{1}{h^2}\begin{cases} u_0(t) - 2 u_1 + u_2, &j = 1\\ u_{j-1} - 2 u_j + u_{j+1}, &1 < j < N-1\\ u_{N-2} - 2 u_{N-1} + u_N(t), &j = N-1 \end{cases} $$ Note that $u_0(t)$ and $u_N(t)$ are given explicitly as $a(t)$ and $b(t)$. This is why I have separated cases $j=1$ and $j = N-1$ in the definition of $F_j$.

Putting this altogether gives $$ r_j = \frac{u_{j-1}(t_n) - 2 u_j(t_n) + u_{j-1}(t_n)}{h^2}\\ s_j = \frac{1}{h^2}\begin{cases} u_0\left(t + \frac{\tau}{2}\right) - 2 \left(u_{1}(t_n) + \frac{\tau}{2}r_{1}\right) + \left(u_2(t_n) + \frac{\tau}{2}r_2\right), &j = 1\\ \left(u_{j-1}(t_n) + \frac{\tau}{2}r_{j-1}\right) - 2 \left(u_j(t_n) + \frac{\tau}{2}r_j\right) + \left(u_{j+1}(t_n) + \frac{\tau}{2}r_{j+1}\right), &1 < j < N-1\\ \left(u_{N-2}(t_n) + \frac{\tau}{2}r_{N-2}\right) - 2 \left(u_{N-1}(t_n) + \frac{\tau}{2}r_{N-1}\right) + u_N\left(t + \frac{\tau}{2}\right), &j = N-1 \end{cases}\\ \frac{u_j(t_{n+1}) - u_j(t_n)}{\tau} = s_j $$ Note that $r_j$ and $s_j$ are helper values to step from $u_j(t_n)$ to $u_j(t_{n+1})$ and are different for each time step. One may want to attribute values $r_j$ and $s_j$ to some moment of time. A reasonable choice would be attributing all each value $\mathbf k_i$ with moment $t_n + \tau c_i$. Here $c_i$ is the first column of the Butcher's tableau.

$$ r_j(t_n) = \frac{u_{j-1}(t_n) - 2 u_j(t_n) + u_{j-1}(t_n)}{h^2}\\ s_j\left(t_n + \frac{\tau}{2}\right) = \frac{1}{h^2} \times \\ \times \begin{cases} u_0\left(t + \frac{\tau}{2}\right) - 2 \left(u_{1}(t_n) + \frac{\tau}{2}r_{1}(t_n)\right) + \left(u_2(t_n) + \frac{\tau}{2}r_2(t_n)\right), &j = 1\\ \left(u_{j-1}(t_n) + \frac{\tau}{2}r_{j-1}(t_n)\right) - 2 \left(u_j(t_n) + \frac{\tau}{2}r_j(t_n)\right) + \left(u_{j+1}(t_n) + \frac{\tau}{2}r_{j+1}(t_n)\right), &1 < j < N-1\\ \left(u_{N-2}(t_n) + \frac{\tau}{2}r_{N-2}(t_n)\right) - 2 \left(u_{N-1}(t_n) + \frac{\tau}{2}r_{N-1}(t_n)\right) + u_N\left(t + \frac{\tau}{2}\right), &j = N-1 \end{cases}\\ \frac{u_j(t_{n+1}) - u_j(t_n)}{\tau} = s_j\left(t_n + \frac{\tau}{2}\right) $$

While this is the answer to the question "How to apply RK2 to this ODE" I really don't like the final form. Instead let's write the same method in a slightly different form: $$ \frac{\mathbf u_{n+1/2} - \mathbf u_n}{\tau / 2} = \mathbf F(t_n, \mathbf u_n)\\ \frac{\mathbf u_{n+1} - \mathbf u_n}{\tau} = \mathbf F(t_n + \frac{\tau}{2}, \mathbf u_{n+1/2}). $$ One can check that the methods are the same by plugging $\mathbf u_{n+1/2} = \mathbf u_n + \frac{\tau}{2} \mathbf r$. Just like $\mathbf r$ and $\mathbf s$ the $\mathbf u_{n+1/2}$ is a helper value used to perform a step from $\mathbf u_n$ to $\mathbf u_{n+1}$.

Applied to our ODE this method gives $$ \frac{u_j(t_{n+1/2}) - u_j(t_n)}{\tau / 2} = \frac{u_{j-1}(t_n) - 2 u_j(t_n) + u_{j-1}(t_n)}{h^2}, \quad j = 1, 2, \dots, N-1\\ u_0(t_{n+1/2}) = a\left(t_n + \frac{\tau}{2}\right), \quad u_N(t_{n+1/2}) = b\left(t_n + \frac{\tau}{2}\right),\\ \frac{u_j(t_{n+1}) - u_j(t_n)}{\tau / 2} = \frac{u_{j-1}(t_{n+1/2}) - 2 u_j(t_{n+1/2}) + u_{j-1}(t_{n+1/2})}{h^2}, \quad j = 1, 2, \dots, N-1\\ u_0(t_{n+1}) = a\left(t_n + \tau\right), \quad u_N(t_{n+1}) = b\left(t_n + \tau\right). $$ Though not strictly necessary I have defined values $u_0(t_{n+1/2})$ and $u_N(t_{n+1/2})$ to get rid of treating $j=1$ and $j=N-1$ as separate cases.