# Runge-Kutta method for PDE

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.

• You're using a biased approximation for second derivative. The scheme would be unstable whatever time integrator you are using. Your numerical solution at $x_j$ depends only on the right part of the domain $x > x_j$, while the true solution depends on the whole domain. Jun 27, 2015 at 12:49

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.

• Oh, thanks a lot.
– newt
Jun 27, 2015 at 16:05
• Works well with the Euler time stepping, but other Runge-Kutta methods are not applicable very well... Nov 16, 2015 at 23:27
• The RK2 method I have been introduced to uses a $k_1$ and $k_2$ "Predictor-Corrector" formula. This just looks a lot like finite difference to me. How does it relate to RK? I have a similar question, perhaps even the same, but I can't understand your answer. Mar 2, 2018 at 1:16
• I'm wondering the same thing. @uranix if you want to apply a RK method of higher order, you'd need to evaluate it at e.g. $t+h/2$ but you don't know this value and cannot evaluate $u$. How can one handle this? Jun 12, 2020 at 6:18
• @vogs Please see my expanded answer Jun 14, 2020 at 13:23