Fundamentals of Plasma Simulation (I)

1. Fundamentals of Plasma Simulation (I) 核融合基礎学（プラズマ・核融合基礎学） /京都大学 李継全（准教授）/岸本泰明（教授）/今寺賢志（D1） 2007.4.9 — 2007.7.13 Part two: Fundamentals of simulation methods ➣ Introduction Importance of simulation & limitations of simulation Developments of computer system & parallel computation ➣ Time advance Initial value problem Numerical stability & CFL condition Explicit methods Semi-implicit & fully implicit methods Predictor/Corrector methods References: T. Tajima, computational plasma physics

2. Importance of numerical simulation Analytical solutions exist only in limit cases. We have established several theoretical approaches to describe plasmas—plasma equations. Generally speaking, we can investigate physical phenomenon by applying these theories if we might have the solutions of these equations. However, plasma physics, like all sciences, is full of high nonlinearity and complexity, analytical solutions exist only in extremely rare cases. Experimental tools are also limited. On the other hand, traditional experiments is limited to measure only a small fraction of the interested quantities with a limited degree of accuracy. Hence, one is usually faced to interpret limited observations with incomplete analytical solutions. Numerical modeling and simulation can bridge the gap between experiments and the ability of theorists. Numerical experiment is to simulate the physical behavior of complicated plasma by solving an appropriate set of equations based on an accepted model.

3. Importance of numerical simulation (cont.) Numerical simulation is popular and powerful in almost all scientific research fields. On the one hand, it includes many advantages of numerical modeling to helpfully solve some traditional problems or reproduce complex physical phenomenon. On the other hand, simulation is essential in studying plasma turbulence and related transport, nonlinear dynamics like chaos and structure formation. Simulation can try innovative idea Most importantly, simulation can solve some problems which is out of analytical solutions or not realistic due to scale in experiments. For example, it is difficult to reproduce space plasma process in laboratory but it is possible to have a numerical observatory. It can perform some innovative ideas to predict new physics.

4. Descriptions and simulations of plasmas Mathematic modeling Description Plasma theories; et al Numerical simulation Discrete algebraic appro.; Numerical algorithms; Codes; et al. Real physical Processes Experiments; et al. Computer simulation may be regarded as the theoretical exercise of numerically solving an initial-value-boundary-value problem. Simulation consists of following the temporal evolution of the configuration. Central task in simulation is the calculation of partial derivatives: ∂nu(t,x,v)/ ∂xn;∂nu(t,x;v)/ ∂tn ; ∂nu(t,x;v)/ ∂vnwith high accuracy. In the following lectures, we will introduce some fundamental methods (schemes) to do space discretization and time advance in numerically solving equation systems.

5. Plasma physics theory Kinetic description (Kinetic/gyro-kinetic theory) fluid description (MHD; two-fluid; …) Gyro-fluid description (with Landau and FLR damping) Fluid ions Kinetic electrons Particle codes Vlasov & Fokker-Planck codes Kinetic/fluid Hybrid code Fluid/MHD code/ Gyrofluid code Targets of plasma simulations Target equations in simulation Vlasov & F-P simulation: (Vlasov code) ● kinetic equations (Vlasov; drift-kinetic; gyrokinetic; F-P) in (x, v) phase space ●gyrokinetic Poisson-Ampere eqs. Particle simulation: (particle code) ● kinetic-Maxwell system (gyrokinetic;drift-kinetic eq.) ● particle motion equation ● poisson-Ampere eqs. Fluid simulation:(fluid code) ● fluid moment equations (MHD; fluid; gyro-fluid)

6. Limitations of numerical simulation Numerical modeling: Numerical simulation results depend on the numerical modeling which should include most essential factors due to the complexity of plasma physics. That is, the physical description of the system by the chosen equations must be accurate for the considered problem (e.g., length and time scales) Computer resource (speed; memory……) is always limited so that the physical size of the problem in time and in space has to be shortened. Simulation uses different statistical average from the strict ensemble averaging. Computer produces only one ensemble per run in general, however, for most statistical arguments, the quantities are statistically averaged ones. It is not practical to perform many runs to make ensemble averaging. If the system is in a quasi-steady state, the average and spatial average are used to replace the ensemble averaging. Numerical accuracy due to the time-spatial discretization. How close can we approach the true description of plasma processes (continuous functions in equation system) by discretizing all the physical processes (discrete grid or particle quantities) for the computer. This is related the numerical discretization methods and also capacity of computer if the accuracy increases linearly with the number of grids or particles.

7. Single scalar processor (CPU) Vector processor parallel supercomputer Developments of computer system & parallel computation S – speed of each processor F – efficiency of the parallelized code Np – number of processors P – proportion (can be parallelized) τC – time for inter-processor communication τS – S /Np Parallel execution time for problem N The speed of computer increases exponentially through the years (from Top500 data). Speed-up is:

8. Summary of plasma equations Particle simulations The microscopic sources are determined by with Or Boltzmann (FP, Vlasov) equation plus equation of motion and Maxwell equations Fluid simulations Continuity equation Equation of motion Energy equation MHD simulations

9. Initial value problem The equations above can be written as the following general form with and boundary conditions Typically, this is an initial value problem. In kinetic modeling (Particle or Vlasov simulations, it includes the dependence of velocity and its derivatives. In plasma fluid (MHD) simulations, it is simplified as with Initial value problem is common in computational physics, also in plasma simulations, for example, particle motion in electric and magnetic field; plasma turbulence and transport simulation; wave-particle interaction; wave propagation; …... We will separate the problem into two issues: time integration and spatial discretization. Before that, we first introduce preliminary computational techniques in plasma simulations.

10. General procedures in simulations To perform a simulation for a given physical problem, we should consider these steps 1. Mathematical formula (equations with initial value and boundary conditions) 2. Numerical discretization (discrete continuous physical quantities in time and space, calculate derivatives). In this step, we obtain a system of algebraic equations. For the spatial discretization, some well-established methods can be applied such as finite-differences; finite element; finite volume; and spectral methods. (Numerical error is produced in this step during discretization!) 3. Time advance. For the time integration, various explicit or implicit schemes have been constructed. If implicit scheme is employed, we need to have a solver. 4. Numerical diagnoses and result analyses to approximate solution.

11. t→i (i+1,j) (i,j-1) (i,j) (i,j+1) (i-1,j) x→ j Initial value problem—1D case Considering an initial value problem represented by partial differential equation with In order to solve (numerically integrate ) this equation by computer, we should discretize it in time and in space, i.e., we must approximate the differentiation by finite difference. For example, simply we can have These expressions reduce to the original differentiation when ∆t tends to zero. Many schemes representing the differential operators can be chosen, which depends on the accuracy and economy in simulation.

12. - - - - - - - - - + + + + + + + + + ions electrons n0 d Numerical stability—Courant-Friederich-Levy(CFL) condition When we choose the finite difference scheme, more important consideration should be its numerical stability in time. We will take the harmonic oscillator as an example to analyze it. This harmonic oscillator may originate from plasma oscillation. Consider an infinite slab of electrons and ion with a width of d (in x ) and particle density of n0. Assume that the electrons are displaced by a small distance x in the x direction. This creates two regions of nonzero charge density. The electric field is produced as The electrons move under the electric field acceleration with an initial position x0(t=0). This initial value problem is described as It has an analytical solution Now we can choose different finite difference scheme to discretize these equations. It will show different numerical stability!

13. Numerical stability—CFL condition (cont.) Firstly, a so-called forward differencing is applied We look for linear stability with an amplification factor g over he adjacent time step to analyze the numerical properties of this equation system Assuming The finite difference equations become The amplification factor g can be solved by setting the determinant of the amplification matrix as zero It can be seen that if |g|2=1+ω2∆t2>1, the solution blows up in time. The amplitude of oscillator is amplified numerically. So, this forward differencing scheme is always numerically unstable for any small ∆t! (This means we can not use it.)

14. t υ i+1 i i+2 t x i+1/2 i+3/2 Numerical stability—CFL condition (cont.) Now we choose another finite difference scheme, the so-called leapfrog scheme Using the same amplification factor g defined in We have the amplification matrix The determinant of amplification matrix equals to zero, we have

15. Numerical stability—CFL condition (cont.) Then If In this case This scheme is stable. If , it is unstable. The stability condition is Time step should be enough small to achieve desirable accuracy. This is the Courant –Friederich—Levy (CFL) condition. Leapfrog scheme is common in particle simulation since two first-order differential equations can be separately integrated for each particle. In simulations, especially in particle simulation with enough more particles, we should consider two aspects when we choose the scheme: — The scheme is as fast a method as possible and still retain acceptable accuracy; — Minimum information needs to be stored when the equation is integrated. On the one hand, it will take more calculation, on the other hand, it will cost much time for data communication among processors in modern parallel simulations. Leapfrog scheme is of advantages in these two aspects compared to high-order scheme like Runge-Kutta method.

16. When τ=0, it reduces to the forward differencing scheme which is numerical stable; when τ=1/2, it is the leapfrog scheme, which is numerical stable for . Biasing scheme For the same differential equation system of plasma oscillator We generalize the leapfrog scheme to introduce a non-integer time step, (i+τ)th time step to push the above equation, similarly discretize the equations to get Similar numerical stability analysis for this scheme. The amplification matrix is When τ=1, the scheme becomes implicit, we have Hence, this scheme is always stable.

17. t→i (i+1,j) (i,j-1) (i,j) (i,j+1) (i-1,j) x→ j Discretization of derivatives in numerical calculation In numerical calculations, there are many different ways to discrete the derivatives. For example for the first order derivative in time and space, we can have And also in the initial value problem with spatial derivatives, we can calculate the spatial derivatives in different time. For example in a diffusion problem We can have two different ways of descretization

18. t υ i+1 i i+2 t x i+1/2 i+3/2 t→i (i+1,j) (i,j-1) (i,j) (i,j+1) (i-1,j) x→ j Old values Old values Explicit and implicit schemes Based on the discretization above, we can see two different schemes: In the first case, the new future value [t=(i+1)∆t] can be given entirely by the old values at t=i ∆t. This is called as explicit scheme. In this case, the calculation is straightforward and easy. In the second case, the new future value [t=(i+1)∆t] should be calculated by the old values and some undetermined new values. This is called as implicit scheme. Since we do not know the undetermined new values, we cannot advance the differential in an explicit way. In this case, we should have a solver to calculate the matrix. So let’s come back to the Bias scheme when τ=0, the forward differencing scheme is explicit method (unstable); when τ=1/2, the leapfrog scheme is also explicit method (but conditionally stable)

19. Old values Old values t t υ υ i+1 i+1 i i+2 i i+2 t t x x i+1/2 i+1/2 i+3/2 i+3/2 Explicit and implicit schemes (cont.) when τ≤1/2, we call this is the backward bias scheme. It is an explicit method; when τ≥1/2, it is called forward bias scheme, it becomes implicit method. Here the reason we call the forward bias scheme with τ≥1/2 as implicit method is because we have to express partially the terms at t=(i+ τ)∆t by using old values at t=i∆t and new values t=(i+ 1)∆t.

20. Low order or high order schemes The accuracy of a finite difference approximation depends on a number of aspects. Higher accuracy can be achieved by the inclusion of more terms in the Taylor expansion and by using preferably centered difference schemes rather then one sided differencing such as backward or forward differencing like follows. Usually, higher order schemes can get a better approximation to the actual result, however, how high order an approximation is necessary depends on numerical purposes. The characteristics of leapfrog scheme are: two-order accuracy; explicit and simple; numerically stable. Taking the plasma oscillator as an example, equations become

21. Low order or high order schemes (cont.) 4th Runge-Kutta scheme is a typical high-order method, which is characterized by: 4th order accuracy; explicit (or implicit); numerically stable. It needs 3~4 times storage of leapfrog scheme. Usual 4th order RK scheme for equation dx/dt=f (t,x) is Comparison between low-order leapfrog scheme and high-order RK scheme: the former gives numerical stability and energy conservation, but can not follow sudden changes in the exact solution (discontinuity). However, RK scheme may cause an increase of energy at a rate determined by the allowed error along trajectories in a non-dissipative system. High-order RK scheme can calculate sharp change of the exact solution much more closely than the leapfrog method.

22. Low order or high order schemes (cont.) Periodic orbit Non-periodic orbit Example: Phase space orbits of pendulum Leapfrog Runge-Kutta Leapfrog scheme correctly gives a closed orbit; RK scheme shows the orbit open up due to error. Leapfrog Runge-Kutta Example: Phase space trajectories of Rayleigh oscillator RK scheme shows true trajectory with change corners; Leapfrog one overshoots sharp change and damps the error.

23. Low order scheme? high order ? Some important aspects should be considered when we choose suitable scheme: 1. Higher order scheme is of high accuracy but also require more computational effort. The choice of the approximation order is one of computational efficiency. 2. A second consideration is how well a chosen grid resolves the actual solution. It requires a good resolution to realize the advantage of a high order scheme. However, a lower order scheme may already provide reasonable results for a good grid resolution. 3. The first order schemes should always be avoided. In many cases it is sufficient to use a scheme of second order accuracy. In terms of numerical efficiency one often has the choice between a simple second order scheme and larger number of grid points or a more complex higher order approximation with less grid points (and thus larger grid spacing). Advantages and disadvantages should be carefully considered, i.e., numerical efficiency to achieve a particular accuracy. The complexity of higher order schemes usually requires more effort in terms of programming and programs are of higher complexity. 4. High order schemes have only limited advantage for discontinuities.

24. Explicit scheme Explicit scheme is traditionally of simplicity in programming code and high speed to push the time integration. In addition, it possesses high parallelization efficiency in modern parallel computation. We will introduce some typical explicit scheme in the application of plasma simulations.

25. Explicit scheme— Particle motion— Time centered scheme Let us consider the particle motion in a static magnetic filed due to Lorentz force In the code, we should reproduce the gyromotion with constant velocity. The equation of motion is discretized in the time-centered scheme as Now we must approximate by the average of vicinal grids Then we have finite difference equations

26. Explicit scheme— Time centered scheme (cont.) If let g=υi+1/υi be the amplification factor, stability analysis shows the amplification matrix Stability equation is It gives the stability condition |g|2=1 when This is similar to the leapfrog stability for the plasma oscillation.

27. t→i (i+1,j) (i,j-1) (i,j) (i,j+1) x→ j Explicit scheme—transport problem (—FTCS scheme) Transport phenomenon such as heat conductivity or magnetic field diffusion are very popular in plasma physics. It is described by a parabolic differential equation Corresponding difference equation for 1D case by using FTCS scheme is FTCS— forward in time (FT) and centered in space (CS) With truncation error The general approach in considering the numerical stability of a scheme is to derive an amplification factor for a Fourier mode and the amplification factor is defined as . It can be derived as CFL condition is |g|2≤1 for all wave-numbers. That is

28. Explicit scheme—transport problem (—2D FTCS scheme) In 2D case, the difference equation becomes CFL condition is An improved scheme can be derived if ∆x=∆y,

29. t→i (i+1,j) (i,j-1) (i,j) (i,j+1) (i-1,j) x→ j Explicit scheme—transport problem (—Leapfrog scheme) —Richardson scheme Difference equation of diffusion problem by using a leapfrog scheme is This is a scheme with centered difference in time and in space. It is called as Richardson scheme. Stability analysis yields the equation of the amplification factor The amplification factor is It can be seen that |g|2>1. This scheme is unconditionally unstable

30. t→i (i+1,j) (i,j-1) (i,j+1) (i-1,j) x→ j Explicit scheme—transport problem (—DuFort-Frankel scheme) The unconditionally unstable Richardson scheme can be modified by splitting into two time levels. It leads to the so-called Dufort-Frankel scheme We can explicitly solve this difference equation as The equation of amplification factor can be derived as The amplification factor is It shows this scheme is stable for all wave-numbers. However, it may give oscillatory solution so that this method in this range is not accurate.

31. means the rate of change at fixed point (Eulerian) means the rate of change at moving point (Lagrangian) the change due to motion Explicit scheme—fluid (MHD) model— convective problem Fluid (MHD) equations are hyperbolic and they are of convection type. This equation system is not straightforward to be integrated in time since they include the convective derivative terms (also called advective term) In Eulerian frame, there is difficulty to integrate in spatially-fixed time due to this nonlinear term. (still current research topics!) On the other hand, one attempts to come to Lagrangian frame to solve some difficulty appeared in Eulerian approach, it also bring in some own new problems (such as grid entanglements as the fluid churns; complex Lagrangian algorithm).

32. Explicit scheme—fluid (MHD) model— Linear convection The convective difficulty associated with the Eulerian algorithm may be looked at in a simple 1D example. All fluid (MHD) equations include such a structure. For example, in density equation, υ is ρ; in momentum and pressure equations, υ is u. At first, we consider the simplest case with constant υ. The analytical solution at time t for an initial value condition f(t=0, x)=F(x) is given by the transport of initial profile, Now we will take numerical approximation. The simplest way is direct discretization by using FTCS scheme Letting and then, Stability analysis shows that the amplification factor is So this is always numerically unstable scheme. This indicates that the algorithm using FTCS scheme cannot be utilized.

33. t→i (i+1,j) (i,j-1) (i,j) x→ j Explicit scheme—fluid (MHD) model— Upwind scheme If we replace the spatial differencing by using upwind scheme is last FTCS method, we get Stability analysis shows that the amplification factor is The stability require (CFL condition) Expanding the upwind scheme in Taylor series This scheme is the first order accurate with leading error term represented by second derivative of x, therefore it represents a numerical diffusion.

34. t→i (i+1,j) (i,j-1) (i,j) (i,j+1) (i-1,j) x→ j Explicit scheme—fluid (MHD) model— leapfrog scheme A rather simple and second order accurate scheme is the so-called Leapfrog scheme Stability analysis shows that the amplification factor is We obtain If , |g|=1. In this case this scheme is numerically stable. On the other hand, , |g|>1, the algorithm is unstable. Thus the linear stability condition (CFL condition) for all wave-numbers is Leapfrog scheme only requires either the oven or the odd gridpoints/timelevels. The differencing decouples odd and even grid points at any given time step such that a solution can develop independently on the interlaced odd/even grid. Thus it may lead to strong oscillations on the grid scale if the two grids are combined. This problem in connection to the issues of diffusion and dispersion.

35. Explicit scheme—fluid (MHD) model— Lax scheme If we replace by a spatially averaged term in the FTCS scheme, we can have the so-called Lax scheme Similarly, we can do stability analysis to get the amplification factor Therefore the stability (CFL) condition is Although this scheme is stable, the accuracy is the first order in space. When the CFL condition is not satisfied, an oscillating solution grows exponentially in time.

36. Explicit scheme—fluid (MHD) model— Lax scheme (cont.) If we rewrite the difference equation The corresponding differential form is An effective diffusion is incurred in this equation. With , the spatial numerical diffusion incurred by this algorithm When , the numerical diffusion coefficient “D” is positive, the scheme is stable. On the contrary, it is unstable. It is very difficult to satisfy the exact condition since the velocity υ is a variable. In addition, the algorithm also causes distortion in the dispersion relation of the eigenmodes.

37. t→i (i+1,j) (i+1/2,j-1) (i+1/2,j+1) (i,j-2) (i,j) (i,j+2) x→ j Explicit scheme—fluid (MHD) model— Lax-Wendroff scheme A scheme closely related to the Leapfrog is the Lax-Wendroff method, which can avoid excessive numerical diffusion and distortion of the dispersion relation in the LAX scheme. In this scheme, we have twice discretization by using Lax scheme. The first step is over a half time step The second step is a full time step Combining two steps, we have (Actually, this is a predictor-corrector scheme, we will explain the details later.)

38. Explicit scheme—fluid (MHD)— Lax-Wendroff scheme (cont.) Stability analysis shows Stability condition |g|2≤1 gives Hence the stability condition gives CFL condition Lax-Wendroff scheme is a second order method in time and in space. The second term in the brace of right side is a diffusion term called the Lax-Wendroff diffusion.

39. Explicit scheme—fluid (MHD)— Lax-Wendroff scheme (cont.) Lax-Wendroff scheme can be derived by using the forward time discretization with a correction that eliminates the lowest order error of the forward time differencing as follows. Expanding the convective equation in time using a Taylor series Substituting the time derivatives from the differential equation Using centered differences for the spatial derivatives Hence this scheme can be regarded as a correction of the FTCS scheme to eliminate the lowest order error of the forward time differencing due to numerical diffusion. From the stability analysis, in the small k∆x limit, Hence the numerical diffusion is fourth order in k∆x , it is small for long wavelength.

40. Implicit scheme The stability of explicit schemes is governed by the CFL condition about the allowable time step. These methods could become inefficient when the equation system becomes high stiff or when the spatial mesh is largely non-uniform. In such cases it is desirable to use implicit scheme. The advantage of this method lies in the numerical stability and large time step.

41. t→i (i+1,j-1) (i+1,j) (i+1,j+1) (i,j) x→ j Implicit scheme – transport problem – FTCS scheme Transport phenomenon is described by a parabolic differential equation FTCS— forward in time (FT) and centered in space (CS) Corresponding difference equation for 1D case by using implicit FTCS scheme is From here, we can see that we must calculate the inversion of a matrix at each time step. This is the general nature of implicit scheme. This method is expensive to use in the code. With truncation error For the amplification factor . Stability analysis shows This indicates this scheme is unconditionally stable for all wave-numbers

42. t→i (i+1,j-1) (i+1,j) (i+1,j+1) (i,j-1) (i,j) (i,j+1) x→ j Implicit scheme – transport problem – Crank-Nicholson scheme In explicit FTCS scheme, if we average the spatial diffusion term in time, the difference equation becomes implicit This scheme produces a second order error in time and in space. The scaling of second order error in time is due to the centered time derivative. The equation of the amplification factor is The amplification factor is This indicates this scheme is unconditionally stable for all wave-numbers.

43. Implicit particle code – first order accurate method In introducing the bias scheme for plasma oscillation, we have known the feature of implicit scheme in solving the equation of motion. Here we consider a 1D unmagnetized homogeneous plasma of finite size particle with only electrostatic interaction. The equation can be discretized by using a fully implicit leapfrog scheme This scheme is restricted by the following time step

44. Implicit particle code – Implicit time filtering The primary goal of the implicit method is to model low frequency phenomena accurately. We will improve the first-order-scheme above so that damping of unwanted high frequency modes increases and the damping of low frequency mode decreases. One way to do so is to replace the fully implicit acceleration in the equation of motion with a time averaged acceleration. This means it will approximately filter out the high frequency responses. One of scheme is It can be seen that these two difference equation are virtually time centered, then preserving the second order accuracy in the leapfrog scheme. The same techniques are used again to solve recursively for the time advanced electric field. This scheme requires more memory and more computation, but has improved frequency filtering over the first order accurate methods. We can do stability analysis to confirm the improvement.

45. Semi-implicit scheme Advantages of implicit schemes: In order to describe the long-time scales of MHD behavior such as the resistive processes, implicit scheme is very useful to remove the time step constraint, by selectively damping out the high frequency components of the field response in the system. However, larger time step causes the narrower spectral width of damping modes. Disadvantages of implicit schemes: On the other hand, implicit scheme requires a large block matrix inversion. Semi-implicit scheme: This method allows a large time step and avoids a large matrix inversion in solving the full implicit equation. The main idea of semi-implicit methods is to add new terms in the time-difference scheme of the MHD equations. These new terms do not affect the solution as the time step tends to zero, and yet the method is absolutely stable relative to fast modes. The step restriction characteristic of ideal MHD is also eliminated. In the semi-implicit method, the terms approximating linear behavior of fast modes are taken implicitly to achieve numerical stability. The method is similar to implicit methods in 1D case. In the multidimensional anisotropic case, the simplicity of the additional terms makes it possible to obtain matrices that consist of three diagonal blocks and can be easily inverted. Because of efficient matrix inversion, the time complexity of the implicit step is the same as for the explicit step. The form of the semi-implicit terms added in various schemes depends on the nature of the fast modes that must be suppressed.

46. Semi-implicit scheme – Alfven wave equation D S Harned & W Kerner, J Comput. Phys. 60, 62(1985) For example, for the fast compressible Alfven wave equation Here b may be magnetic field; υis fluid velocity, a is the Alfven velocity. Considering discretizing in time, wave equation becomes This is an explicit scheme with time step CFL restriction for fast mode. To make the scheme become unconditional stable, we can move the second spatial derivative term to the time step (i+1), but it becomes an implicit scheme. So we add a new term in both sides to make them implicit on the LHS implicit and explicit on the RHS to avoid the calculation of the inversion of matrix.

47. Semi-implicit scheme – Alfven wave equation (cont.) Where a0 is a constant, which is determined the stability condition. The subtraction of these new terms can have the same effect as making the fast modes implicit, even if it may be very different from the Alfven velocity a. This method is unconditionally stable when a0≥a. This is equivalent to adding a linear term with a coefficient proportional to ∆t in original wave equation Stability analysis shows it unconditional stable provided that

48. Predictor-corrector scheme The previous schemes like leapfrog and Runge-Kutta scheme are called single-step methods because they use only the information from one previous point to compute the successive point. that is, only the initial point (i=0, j) is used to compute (i+1, j), and in general, (I, j) is needed to compute (i+1, j). After several points have been found, it is feasible to use several prior points in the calculation. We can develop multiple step method. A desirable feature of a multi-step method is that the local truncation error can be determined and a correction term can be included, which improves the accuracy of the answer at each step. Also, it is possible to determine if the step size is small enough to obtain an accurate value for (i+1, j), yet large enough so that unnecessary and time-consuming calculations are eliminated. Actually, we have met this treatment (predictor-corrector) in Lax-Wendroff scheme by a two-step procedure. The first is over a half time step (predictor): The second is a full time step (corrector):

49. Predictor-corrector scheme – Alfven wave equations In applying the semi-implicit scheme for the Alfven wave problem, a practical way to advance the differential equation in a code is a two-step predictor-corrector scheme. From the wave equations, by writing the time partial derivative as , the predicted values are The corrected values are Or they can be written as Usual stability analysis gives the equation of amplification factor The explicit system of this predictor-corrector equations is stable if For θ=1/2, the scheme is unstable for all ∆t, for θ=1, is stable when ∆t <1/a.

50. Predictor-corrector scheme – Alfven wave equations (cont.) D S Harned & W Kerner, J Comput. Phys. 65, 57(1986) If we try to modify this predictor-corrector scheme according to the semi-implicit technique, The predicted values are The corrected values are a0 is some constant. For the velocity equation, we have With 0.5≤θ≤1. This method is unconditionally stable if