How do I determine what $f(y, t)$ is when using Runge-Kutta to solve a PDE?

Question

I've been given the formula for RK2 and shown how to use it on an ODE to solve $y'(t) = f(y, t)$, where $f(y,t)$ is explicitly known.

Equipped only with this knowledge I have been asked to use RK2 on a PDE, where neither side of the equation is known, except that I see that I must substitute the spacial derivative with a finite difference operator. I am required, for the purpose of this assignment to use a center space finite difference, as shown below.

My work so far:

Given $\partial_tu(x,t) = -c\partial_xu(x,t)$ is analogous to $y'(t) = f(y,t)$ when we consider all variables besides t as constants:

$$ k_1 = \Delta t f(y_n, t_n) \\k_2 = \Delta t f(y_n + k_1, t_n + \Delta t)\\ y_{n+1} = y_n + \frac{1}{2}(k_1 +k_2) + \mathcal{O}(\Delta t^3) $$

Becomes:

$$ k_1 = \Delta t f(u_{j+1}, u_{j-1}) = -c\Delta t \frac{u^n_{j+1} - u^n_{j-1}}{2\Delta x} \\k_2 = \Delta t f(u_n +k_1, t_n + \Delta t) = ????? \\ u^{n+1}_j = u^n_j + \frac{1}{2}(k_1 +k_2) + \mathcal{O}(\Delta t^3) $$ $$ \text{where }u^n_j \text{ means } u(x_0 + j\Delta x, n\Delta t) $$

I have no idea how to add $k_1$ to my new "function" of $u_{j+1}$ and $u_{j-1}$ in order to find $k_2$. My textbook (Numerical Recipes [in fortran 77!]) does not review RK2 in combination with PDEs, and neither do my class notes. I have no idea how to functionally put RK2 to work on a PDE without knowing this step. I am not far enough along in the course to make use of linear algebra, though I have seen suggestions to do so.

This question seems quite similar, but I do not understand the accepted answer:

Runge-Kutta method for PDE

The answer seems to just put finite difference to use. I don't see Runge Kutta used at all. It looks like a finite difference operator is used and that's it.

Answer 1

One approach for performing such time-integration is the method of lines (MOL). Let us introduce the nodal values $u_j(t) = u(x_j,t) = u(x_0 + j\Delta x,t)$. Using centered differencing to approximate the spatial derivative, the PDE gives \begin{aligned} \frac{\text d}{\text{d} t} u_j(t) &= -c\, \partial_x u(x_j,t) \\ & \simeq -c\, \frac{u_{j+1}(t) - u_{j-1}(t)}{2\, \Delta x} \, . \end{aligned} The system of ordinary differential equations \begin{aligned} \frac{\text d}{\text{d} t} {\bf u}(t) &= -\frac{c}{2\, \Delta x} \left( \begin{array}{ccccc} 0 & 1 & &\\ -1 & \ddots & \ddots &\\ & \ddots & \ddots & 1\\ & & -1 & 0 \end{array} \right) \, {\bf u}(t) - \frac{c}{2\, \Delta x} \left( \begin{array}{c} -u_0(t)\\ 0\\ \vdots\\ 0\\ u_N(t) \end{array} \right) ,\\ &= {\bf f}({\bf u}(t)) \, , \end{aligned} made by the vector ${\bf u} = (u_1,\dots ,u_{N-1})^\top$ of nodal values can then be integrated in time explicitly using Runge-Kutta methods.

In particular, if the forward Euler (RK1) method ${\bf u}^{n+1} = {\bf u}^{n} + \Delta t\, {\bf f}({\bf u}^n)$ is used for time-integration, then an unstable scheme $$ u_j^{n+1} = u_j^n -c\, \frac{\Delta t}{2\,\Delta x} (u_{j+1}^n - u_{j-1}^n)\, , \qquad 1 \leq j \leq N-1 $$ is obtained. A slight modification of this method gives the stable Lax-Friedrichs scheme.

If the proposed improved Euler (RK2) method is used for time integration, one has \begin{aligned} \tilde u_j^{n+1} &= u_j^n + \Delta t\, k_1 \\ &= u_j^n -c\, \frac{\Delta t}{2\, \Delta x} (u_{j+1}^n - u_{j-1}^n) \, ,\\ u_j^{n+1} &= u_j^n + \frac{1}{2} (k_1+k_2)\\ &= u_j^n -c\, \frac{\Delta t}{4\, \Delta x} (u_{j+1}^n - u_{j-1}^n + \tilde u_{j+1}^{n+1} - \tilde u_{j-1}^{n+1}) \, . \end{aligned} Finally, for $2 \leq j \leq N-2$, $$ u_j^{n+1} = u_j^n -c\, \frac{\Delta t}{2\, \Delta x} (u_{j+1}^n - u_{j-1}^n) + c^2\, \frac{\Delta t^2}{8\, \Delta x^2}(u_{j+2}^n - 2u_{j}^n + u_{j-2}^n) \, , $$ which looks quite similar to the Lax-Wendroff scheme, but may have different stability properties.