The hyperdiffusion method is a technique for the numerical solution of partial differential equations (PDEs) in a way which is numerically stable. It is known that physically interesting problems involve shocked and unstable systems, obtaining stable solutions for such systems may be numerically challenging. We review numerical methods for solving PDEs and attempt to understand how numerical instabilities arise. This is aided with the presentation of an example program with matlab and scilab to understand numerical instability and the different stabilisation methods.

The method known as finite differencing is used to solve PDEs when analytical methods are not available.The approach adopted is to map out a problem onto a mesh of finite elements. For example, consider the first order differential equation

A number of finite differencing techniques are available for solving this.

The first order Euler method solves this as

y_(n+1) and y_(n) are values of y at the n th and (n+1)th mesh point, h is the difference between the two points

f(x_n, y_n) is the function evaluated at the nth point. To second order the Euler method is written

Here, we have used the notation defined above, a third order Euler method also available. One of the most powerful soultion techniques for first order PDEs is the fourth order Runge-Kutta method, using this method the solution may be obtained using the following expressions.

Where,

Now consider more general PDEs of the form

Equation 10: General PDE |

Consider a Taylor expansion around the point x, so we may write

Manipulating the first expansion we obtain the derivative of F in terms of the forward difference

the derivative can also be obtained using the backward difference
Taking the difference between the two expansions, gives the derivative in terms of the central difference

Notice that we have neglected the higher order terms in the Taylor expansion. Central differencing formula for the derivatives to order h^2 are as follows:

The notation used above is as follows

and for example

The central differencing expressions to order h^4 for the first three derivatives is as follows

Now, solve equation 10 using first the forward difference

and then a backward difference

Combining the forward and backward difference we can solve the PDE using the following

This method is known as the Crank-Nicolson method. It is a so called implicit method in that knowledge of the solution at time t_{i} gives a relationship with the function at time t_{i+1}, which is not given explicitly. We now apply these methods to the solution of the following advection problem.

Equation 29: Advection problem |

For the above equation, when c=u, then we have an equation called the Burgers equation, this so called inviscid equation can generate solutions which can develop discontinuities and shock waves. In our matlab and scilab script we attempt to solve this numerically for a propagating square wave and a sine wave. Solutions for the central differencing method at two different times are shown below.

Advection of a square wave at time t1 |

Advection of a square wave at time t2 |

The first image shows the initial square wave and how it is advected after a small number of time steps, it can be seen clearly that numerical disturbances arise at the leading and trailing edges of the wave form. At the even later time t2, the noise has grown and it is apparent that these solutions are numerically unstable. Why is this?

To obtain stabilised solutions a number of solution approaches are possible including the hyperdiffusion approach, but first we consider another approach to obtaining numerically stable solutions by considering the Lax-Wendroff method. Consider the following PDE.

To obtain numerically stable solutions of such a problem we may use two stages. In the first stage the solution is advanced by half a time step and by half a mesh point in the forward/backward direction as follows

We advance the solution using the results from the first step as follows

Using this approach our advection problem, equation 29 may be solved using the following finite difference expression.

Using VonNeumann stability analysis the stability of the solutions can be checked, with solutions of the form

For the lax-Wendroff method the amplification factor A is found to be 1. If A<1 the amplitude for numerical diffusion is found to decrease. If A>1 the amplitude for numerical diffusion is found to increase and the solution is numerically unstable. The result of the application of this solution method is shown below.

Advection of a square wave at time t1 using the Lax-Wendroff method |

Advection of a square wave at time t2 using the Lax-Wendroff method |

Using the matlab script file comment out the lines for the central differencing method and the hyperdiffusion corrections, Also comment out the correct lines before the statement, for i = starti:finishi. It is necessary to set the correct values for starti and finishi, for the Lax-Wendroff method we use the second order differencing values. When this script is run we should observe the advection of the square wave as illustrated by the two figures shown above. Once again, observe the numerical instability at the leading/trailing edges of the square wave. However at the later time t2 although there is still numerical noise this has not grown significantly and we can successfully compute solutions at later times.

We now consider the source of the numerical noise and how the hyperdiffusion method may be used to control that noise.

To help understand the source of the numerical noise consider the solution of the following PDE

let u be cosntant and solve the problem in 1D, so that we write

Using a Taylor expansion we can write the follwoing Taylor expansions

Taking the difference of these two expressions our finite difference version of the PDE becomes

If we choose a solution of the form

then we can obtain an expression for the frequency as follows

It can be seen that the second term in the series causes a reduction of the frequency and the solution will change shape. The phase speed of the signal is

Shock waves and square waves may become dominated by high frequencies. Therefore, f we can use a method, such as central differencing to truncate the solution and apply this near the discontinuity it may turn out that the solution is actually dominated by the higher frequencies. An nth order term may be written as

For the 1D problem the second order correction may be written as

The introduction of the terms such as this leads to a damping term which can be used to damp out some of the higher frequency terms of a given wavelength. Solving the above equation solutions of the following form are obtained

The coefficient nu_2 damps out the spatial wavelengths.

Now consider the 1D problem

In this problem we have numerical discontinuities arising from shock, advection and sound propagation. If the introduced term is adjusted to balance the jumps in the solution arising from the higher order terms then a solution may be stabilised. Following the definition of our damping term earlier a general expression for the viscosity is as follows

The problem now is to set the coefficients so that the discontinuities are successfully cancelled. So consider our solution of the form

If a wave of length 3delta x dies out quickly then after 4 lengths it should be cancelled.

The time t_e is the time such that the signal will decay to 1/e of its original what must nu_2 be in order to satisfy this.

If we define

Now, define

Using these definitions we obtain

i.e.

Thus our differential equation with numerical correction becomes

We have computed the value of v_2 to enable the correct cancellation of the higher order terms. We can compute the wave speed and we can compute delta x from the Courant condition. Now consider the equation of momentum conservation from our perturbed system of equations for ideal MHD. Below rho_b is the fixed background density and rho is the perturbed density. For full details see reference 1 below.

We will now consider how we compute the hyperdiffusion corrections for this system of equations

where

The above expression takes the same form as the corrective term we described earlier, the second term is a correction for the shock viscosity.

delta^{3} is a third order difference and delta^{1} is a first order difference, i.e. using central differencing,

These expressions will be used to pick out the numerically unstable parts of the solution for our computational mesh. For our MHD code we split these contributions into a left handed and right handed differencing contributions.

We express this numerically as

The maximum is taken over the left and right differences over the computational mesh. We have implemented this term in the example matlab and scilab codes, here each of the terms are named d3l,d1l,, d3r and d1r.

The shock viscosity term is calculated for negative compression, the compression term is

Thus

Only regions undergoing compression are corrected using the shock viscosity. The hyperdiffusive operator

where

c_s is the sound speed, v_a is the Alfven speed.

is the hyperdiffusion function and

We may write

These terms have been implemented in our matlab/scilab code through the terms labelled nushk(i) for each of the mesh points.

Advection of a square wave at time t1 using the hyperdiffusion method |

Advection of a square wave at time t2 using the hyperdiffusion method |

We now consider the hyperdiffusion terms for the equation of continuity and the evolution of energy.

For the energy term

where

Where the thermal energy perturbation is used and given by,

In the expression above, epsilon_klm is the levi-civita symbol

### Codes

matlab example script runanmstringhyper1.mscilab example script runanmstringhyper1.sci

### References

- Magnetohydrodynamic code for gravitationally-stratified media, S. Shelyag, V. Fedun and R. Erdelyi, A&A 486, 655–662 (2008)
- Numerical Computation with Shocks - Sondergaard
- A 3D MHD model of astrophysical flows: Algorithms, testsand parallelisation - S. E. Caunt and M. J. Korpi - A&A 369, 706-728 (2001)
- SIMULATIONS OF SOLAR GRANULATION. I. GENERAL PROPERTIES, R. F. STEIN ANDÃ“. NORDLUND, ApJ. 499:914-933, 1998 June 1

## No comments:

## Post a Comment