Random walk

Increasing the order of a numerical scheme is important as soon as non-linearities are contained into an Ordinary Differential Equation (ODE) or a Partial Differential Equation (PDE).

In this example, I give the recipe to build an eulerian upwind numerical scheme, useful to compute an integration in time, for an evolution equation of the type \[ \frac{\partial T}{\partial t} = f(T) \].

Of course, it could be of interest to use an old good 4th order Runge Kutta scheme, however, as soon as you have a partial differential equation, then the RK4 scheme does not stand properly anymore. In that case, a high order integration scheme can be of interest.

The important questions are :

• How to find the expression of the numerical approximation of the derivative at the timestep number $n$?
• Which precision order can I reach if I use N variables to calculate it ?

In the multiscale simulations codes, it is rather typical to change the integration time step as a function of time, then we will directly generalize the procedure to a set of different timesteps.

As an example, let's decide to use the variables $T(t)$, $T(t-h_1)$, and $T(t-h_2)$, where $h_i$ are the time between $t$ and the previous timesteps. The variable timesteps will be calculated at the end of the process.

1. We start by developing the expression around $T(t)$ we are interested in. \[ T(t-h_1) = T(t) - h_1 \times \frac{\partial T}{\partial t} + \frac{h_1^2}{2!} \times \frac{\partial^2 T}{\partial t^2} \] \[ T(t-h_2) = T(t) - h_2 \times \frac{\partial T}{\partial t} + \frac{h_2^2}{2!} \times \frac{\partial^2 T}{\partial t^2} \]
2. We want to build a linear combination of $T(t)$, $T(t-h_1)$, and $T(t-h_2)$ to build an expression equal to $\frac{\partial T}{\partial t}$, since we want to numerically express the derivative via three different points in time. \[ \alpha_1 T(t) + \alpha_2 T(t-h_1) + \alpha_3 T(t-h_2) = \frac{\partial T}{\partial t} \] where $\alpha_1$, $\alpha_2$ and $\alpha_3$ belong the real-valued space $\mathcal{R}$.
3. Then, we count by identification the number of $\alpha_i$ needed to obtain one $\frac{\partial T}{\partial t}$, zero $\frac{\partial^2 T}{\partial t^2}$ and zero $T(t)$. This operations leads to write the following matrix \[ A = \left[ \left( 1, 1, 1 \right), \left(0, -h_1, -h_2 \right), \left(0, h_1^2, h_2^2 \right) \right] \].
4. In the basis $\left(T(t), \frac{\partial T}{\partial t}, \frac{1}{2!} \frac{\partial^2 T}{\partial t^2} \right)$, to define the expression that we want to build, we define the vector \[ B = ( 0, 1, 0) \] and then solve the system \[ A X = B \] where \[ X=\left( \alpha_1, \alpha_2, \alpha_3 \right) \]
5. This leads to the generalized solution \[ X = \left( \frac{h_1 + h_2}{h_1 h_2}, \frac{h_2}{h_1 \left( h_1 - h_2 \right)}, \frac{h_1}{h_2 \left( h_2 - h_1 \right)} \right) \].
6. We have obtained the values of the $\alpha_i$ allowing to express the time derivative using the ingredients that we wanted to put inside. For convenience into the numerical code, it just remains to shape the final expression as a function of temporal interval between the numerical steps $dt_1$ and $dt_2$. We thus define and replace in the previous expression \[ h_1 = dt_1 \\ h_2 = dt_1 + dt_2 \].

We obtain the expression of the numerical derivative $\frac{\partial T}{\partial t}$ at the time step number $n$ as a function of the values of the solution $T(t)$ at the previous timesteps $T_{n-1}$ and $T_{n-2}$.

To draw a conclusion, we have seen that, by considering N different points to build an expression of \frac{\partial T}{\partial t}, we have build a value of a precision order N-1.

As the last stage, we express the numerical solution $T(t)$ at the time step number $n+1$ via the following formula \[ T_{n+1} = \frac{dt_1 \left( dt_1 + dt_2 \right) }{2 dt_1 + dt_2} \left( \frac{dt_1 + dt_2}{dt_1 dt_2} T_{n} - \frac{dt_1}{\left( dt_1 + dt_2 \right) dt_2 } T_{n-1} + f(T) \right) \]

An eulerian upwind scheme of the second order precision has been successfully built and can be directly implemented in a code. This procedure can be generalized to use N coefficients, which leads to a precision of the N-1 th order in the case of a first order derivative.