Wednesday, 4 February 2015

Stable Numerical Solution of PDEs Using Hyperdiffusion

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
The pictures above illustrate the advection of a square wave using our matlab/scilab code with the implementation of the hyperdiffusion method described here. We can see that the correction improves on the Lax-Wendroff approach and is effective at preserving the wave shape.

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

No comments:

Post a Comment