電漿物理中的不穩定性過程可以使用質點網格法(Particle-In-Cell PIC)進行模擬，由於不穩定性的物理過程將伴隨著較劇烈的能量交換，在數值方法的局限之下，極有可能隨之產生更大的數值偏差，又因不穩定性的過程中必須滿足能量守恆，因此我們將紀錄系統每個時刻的粒子的動能與電場的能量，藉由分析系統的總能量，來評估不穩定性過程所產生的數值誤差是否可以忽略。
本篇論文將以不同數值近似法來進行比較，不穩定性的發生時，其所產生的誤差是否得到更好的控制，我們將動量方程式分別以Euler法、Velocity Verlet法、RK2法、RK4法進行比較，而空間電場的近似分別以二階、四階、六階的中央差分法進行比較。
而模擬的結果顯示，以Velocity Verlet法、RK2法、RK4法近似，在總能量的分析上，都維持不錯的能量守恆，然而透過模擬顯示，我們也確實觀察到，PIC模擬在不穩定性的發生過程中總能量確實存在的較大的數值偏差，雖在此時刻總能量的偏差相較於動能與電能變化，在程度上是可以忽略，但在不穩定性發生時總能量所產生的數值偏差，也會影響著後續能量變化的發展，透過模擬的比較中我們也得知，盡可能地提升PIC每個方程式的準確度，可以獲得較有效的改善，可以避免使用過小的時間步長(time step)，而造成計算時間的增加。

Abstract

Energy Error Analysis of Particle-In-Cell Simulations With Different Numerical Methods

Author：Shao-Wei Hong
Advisor：S.W.Y.Tam
Institute of Space and Plasma Sciences, National Cheng Kung University
SUMMARY
From our simulations, the results show the Velocity Verlet method, RK2 method, and RK4 method all have a good accuracy for energy conservation. However we also observed a larger numerical error during the instability process. Although the error of total energy compared with the variation of energy change between the kinetic energy and the electrical energy is negligible, the numerical deviation would affect the subsequent numerical results. One needs to reduce the error as much as possible. Our simulation results show that using higher order numerical methods to approximate the momentum equation and reducing the time step can decrease the error during the plasma instability process.

Key word：plasma instability, energy conservation, Particle-In-Cell (PIC) method

INTRODUCTION
The Particle-In-Cell (PIC) method can be used to simulate the instability process in plasma physics. The limitations of the numerical method would cause the numerical errors to increase. Since the plasma instability process is associated with intense energy exchange and must satisfy the energy conservation, we record the total energy of the system at every time step, to evaluate the size of the numerical errors.

MATERIALS AND METHODS
We examine different numerical methods for their accuracy in PIC simulations: Euler method, Velocity Verlet method, RK2 method, and RK4 method. Spatially, we compare second-order, fourth-order, and sixth-order central finite-difference method for the Poisson equation.

RESULTS AND DISCUSSION
We use the Particle-In-Cell method to simulate the bump-on-tail instability process, and the results show an intense energy change during the instability process: the particles lose energy, and wave energy grows, and we can see the bump in the velocity distribution being smooth out. These results correctly correspond to plasma physics. Moreover we use the Fourier transform to get the frequency band and the wave number of the growing waves, and use the relation v=ω/k to calculate the corresponding phase velocity in the velocity distribution. We found that the range of phase velocities corresponds to the velocity range of positive slope, so the PIC method can correctly simulate plasma instability. From simulations, we know there is a larger numerical error during the instability process, and the error deviation can be reduced as we improve the accuracy of the approximation for any equation, or reduce the time step.

CONCLUSION
In our study, the Velocity Verlet method, RK2 method, and RK4 method are numerically stable, and we have to reduce all errors as much as possible to avoid the accumulation of errors affecting the subsequent numerical accuracy significantly. Improving the order of numerical method and reducing the time step are both effective ways to reduce the numerical errors. However, practically, to improve the order of numerical method is a more effective way to improve the accuracy of the entire system, because to reduce the time step increase the overall number of steps of the calculation, causing the round-off error to occur more often.

[1] C.K. Birdsall, A.B. Langdon,Plasma physics via computer simulation. New York : Taylor & Francis, 2005
[2]https://www.physics.drexel.edu/~valliere/PHYS305/Diff_Eq_Integrators/Verlet_Methods/Diffrntleqn3.pdf
[3] P. O.J. Scherer, Computational physics : simulation of classical and quantum systems. New York : Springer, c2013.
[4] K. Hussain, F. Ismail, N, Senua., Solving directly special fourth-order ordinary differential equations using Runge-Kutta type method, Journal of Computational and Applied Mathematics 306 (2016) 179-199
[5] http://fourier.eng.hmc.edu/e176/lectures/NM/node42.html
[6] J.C.Butcher.,Coefficients for the study of Runge-Kutta integration processes. J.Austral.Math.Soc.3 (1963), 185–201.
[7] J. C. Butcher, Numerical Methods for Ordinary Differential Equations Second Edition. 2016, John Wiley & Sons, Ltd
[8] M.Tenenbaum,H.Pollard. ,Ordinary Differential Equations. Dover Pubns Chapter 12 pp. 763-765 ,1985
[9]S. Francis, Schaum’s Outlines of Theory and Problems of Numerical Analysis New York, NY: Mcgraw-Hill,Chapter 20 pp. 232-233, 1989
[10] F.F.Chen,Introduction to Plasma Physics and Controlled Fusion.New York and London:Plenum Press, second edition, 1983.
[11]R.L.Burden, ,Numerical Analysis.9th Edition Brooks/Cole 2011