Simulation Scenarios in One-Dimensional Self-Gravitating System
SSimulation Scenarios in One-Dimensional Self-GravitatingSystem
A guided study for the partial fulfillment of Ph.D. courseworkby
Suman Pramanick
Advisor
Somnath Bharadwaj
Department of PhysicsIndian Institute of Technology KharagpurKharagpur, West Bengal 721302, India a r X i v : . [ phy s i c s . c o m p - ph ] J a n BSTRACT
We study and compare different numerical differential equation solvers on the basisof numerical complexity, energy conservation, and stable solution in phase-space forthe Simple Harmonic Oscillation (SHM) problem. We conclude and show that theLeapfrog method is the best for our problem. We solve Poisson’s equation in gravityon the computation grid by the finite difference method and Fourier method. We solvePoisson’s equation for source not in a grid point by cloud-in cell method. Finally,we simulate one-dimensional self-gravitating system and show the evolution of thesystem via phase-space trajectories. ii ontents
Bibliography 42 iv hapter 1Solutions by Euler, RK-2 andRK-4 Methods To Solve the Simple Harmonic Oscillation (SHM) problem by Euler method, Runge-Kutta second order (RK-2) method and Runge-Kutta fourth order (RK-4) method.
There are very few real systems which can be solve analytically, safely we can say mostof the real life systems we can not solve analytically. However we can formulate thesystem mathematically and can write the differential equations to solve the system.These systems can be solved using numerical methods with obviously some boundaryconditions specified. The solution we get from numerical techniques are not errorfree, because in every numerical technique we assume some simple steps to reducethe complexity of the system. Those assumptions are good approximation of the realsystem but not exactly the real system. That leads to a accumulation of errors in eachcomputational steps (iterations). For a long run these errors sum up and sometimemake the numerical system way deviate from the real system.There are two things we need keep in mind. We have to take each infinitesimalsteps such that they don not deviate much from the real system. The second point is todo such we cannot suggest an infinitesimal process which makes a huge computationalcomplexity, because that will make more run time and sometimes less efficient. So,1asically we have to balance between these two things, and have to choose the mostefficient technique according to our requirements.
In this assignment we will solve the simple harmonic oscillation (SHM) problem. Theequations of the problem are: dxdt = pm and dpdt = − kx (1.1)The parameters for our problem are: k = 1; m = 1; ω = 1 and T = 2 π (1.2)Before solving the equations with different numerical techniques we will brieflydiscuss about the details of some numerical techniques to solve differential equations. Euler method is a first-order method of solving differential equation. The error perstep in Euler method is proportional to the square of the step size. The global errorof Euler method is however proportional to the step size. This is the simplest andbasic method of solving differential equation. (a) fig 1 (b) fig 2
Figure 1.1: Euler method approximate solution2ccording to figure 1.1, we can see that at point ( x , y ) the slope (slope of thetangent at that point) is dydx . lets say dydx = f ( x, y ) then at point ( x , y ) the slopeis f ( x , y ). for a small increment h we can approximate y + hf ( x , y ) as the newvalue of y at point x + h . So, we can write the iterative formula as x n +1 = x n + h (1.3) y n +1 = y n + hf ( x n , y n ) (1.4)So, we need the initial values of x and y to compute the whole function in iterativemethod. If we want other points along the path of the true solution, and yet wedon’t actually have the true solution, then it looks like using the tangent line as anapproximation might be our best bet! After all, at least on this picture, it looks likethe line stays pretty close to the curve if you don’t move too far away from the initialpoint.In this way according to figure 1.1 we can see that we can at least trace thefunction and if the step size h is very small then we can actually very accuratelycalculate the function. Now let us discuss the SHM problem which is our main goal in this chapter.Form figure 1.2 we can see that the solution is diverging for small N values butfor large N , like N = 10000 the solution is very close to the actual elliptical solution.For small N the step size h is very small, so error accumulation is large in each step.The energy of the system should be constant as the system is conservative and hasno damping. From figure 1.4 we can see the energy is diverging for small N . Howeverfor large N the energy is comparably constant of the process. Runge-Kutta methods are among the most popular ODE solvers. They were firststudied by Carle Runge and Martin Kutta around 1900. In contrast to the multistepmethods of the previous section, Runge-Kutta methods are single-step methods —however, with multiple stages per step. They are motivated by the dependence of theTaylor methods on the specific IVP. These new methods do not require derivatives of3
Position (x) -5051015 M o m en t u m ( p ) Phase Space Trajectory For N=10 by Euler Method (a) N = 10 -5 0 5 10 Position (x) -8-6-4-202468 M o m en t u m ( p ) Phase Space Trajectory For N=100 by Euler Method (b) N = 100 -1.5 -1 -0.5 0 0.5 1 1.5 Position (x) -1.5-1-0.500.511.5 M o m en t u m ( p ) Phase Space Trajectory For N=1000 by Euler Method (c) N = 1000 -1.5 -1 -0.5 0 0.5 1 1.5 Position (x) -1.5-1-0.500.511.5 M o m en t u m ( p ) Phase Space Trajectory For N=10000 by Euler Method (d) N = 10000 Figure 1.2: Phase space trajectories for different number of steps ( N ) inside a fulltime period T (For Euler method)the right-hand side function f in the code, and are therefore general-purpose initialvalue problem solvers. Let us consider dydx = f ( x, y ) (1.5)The Taylor expansion of y ( x ) is y ( x + h ) = y ( x ) + hy (cid:48) ( x ) + h y (cid:48)(cid:48) ( x ) + O ( h ) (1.6)The second derivative of y in terms of f is4
10 20 30 40 50 60 70
Time (t) -8-6-4-202468 P o s i t i on ( x ) particle position by Euler Method For N=10 (a) N = 10 Time (t) -8-6-4-202468 P o s i t i on ( x ) particle position by Euler Method For N=100 (b) N = 100 Time (t) -1.5-1-0.500.511.5 P o s i t i on ( x ) particle position by Euler Method For N=1000 (c) N = 1000 Time (t) -1.5-1-0.500.511.5 P o s i t i on ( x ) particle position by Euler Method For N=10000 (d) N = 10000 Figure 1.3: Solution of x ( t ) vs t plot for different number of steps ( N ) inside a fulltime period T (For Euler method) y (cid:48)(cid:48) ( x ) = f x ( x, y ) + f y ( x, y ) f ( x, y ) (1.7) f y is the Jacobian. the Taylor series expansion of y becomes y ( x + h ) = y ( x ) + hf ( x, y ) + h f x ( x, y ) + f y ( x, y ) f ( x, y )] + O ( h ) (1.8) y ( x + h ) = y ( x ) + h f ( x, y ) + h f x ( x, y ) + hf x ( x, y ) + hf y ( x, y ) f ( x, y )] + O ( h ) (1.9)From multivariate Taylor series expansion f ( x + h, y + k ) = f ( x, y ) + hf x ( x, y ) + f y ( x, y ) k + ... (1.10)5
10 20 30 40 50 60 70
Time (t) T o t a l ene r g y ( E ) Euler Method
N=10 (a) N = 10 Time (t) T o t a l ene r g y ( E ) Euler Method
N=100 (b) N = 100 Time (t) T o t a l ene r g y ( E ) Euler Method
N=1000N=10000N=100000 (c)
Figure 1.4: Constant energy value of the solution over time (For Euler method).So, f ( x + h, y + hf ( x, y )) = f ( x, y ) + hf x ( x, y ) + hf y ( x, y ) f ( x, y ) + O ( h ) (1.11)Therefore, we get y ( x + h ) = y ( x ) + h f ( x, y ) + h f ( x + h, y + hf ( x, y ))] + O ( h ) (1.12) y n +1 = y n + h (cid:18) k + 12 k (cid:19) Where, k = f ( x n , y n ) (1.13) k = f ( x n + h, y n + hk ) (1.14)6 .5.2 Explanations of plots From figure 1.5 we can see for small value of N the phase space solution is diverging.However for larger values of N we get very close solution to actual elliptical one.From figure 1.6 we can see that the solution of position x ( t ) is diverging for small N , but for large N we get very stable solution.The constancy of total energy value also can be observed from figure 1.7. Howeverfor small N energy value is not constant after some iteration, but for large N valuesit is well constant. -6 -4 -2 0 2 4 Position (x) -8-6-4-20246 M o m en t u m ( p ) Phase Space Trajectory For N=10 by RK-2 Method (a) N = 10 -1.5 -1 -0.5 0 0.5 1 1.5 Position (x) -1.5-1-0.500.511.5 M o m en t u m ( p ) Phase Space Trajectory For N=100 by RK-2 Method (b) N = 100 -1.5 -1 -0.5 0 0.5 1 1.5 Position (x) -1.5-1-0.500.511.5 M o m en t u m ( p ) Phase Space Trajectory For N=1000 by RK-2 Method (c) N = 1000 -1.5 -1 -0.5 0 0.5 1 1.5 Position (x) -1.5-1-0.500.511.5 M o m en t u m ( p ) Phase Space Trajectory For N=10000 by RK-2 Method (d) N = 10000 Figure 1.5: Phase space trajectories for different number of steps ( N ) inside a fulltime period T (For second-order Runge-Kutta Method)7
10 20 30 40 50 60
Time (t) -8-6-4-20246 P o s i t i on ( x ) Particle position by RK-2 For N=10 (a) N = 10 Time (t) -1.5-1-0.500.511.5 P o s i t i on ( x ) Particle position by RK-2 For N=100 (b) N = 100 Time (t) -1.5-1-0.500.511.5 P o s i t i on ( x ) Particle position by RK-2 For N=1000 (c) N = 1000 Time (t) -1.5-1-0.500.511.5 P o s i t i on ( x ) Particle position by RK-2 For N=10000 (d) N = 10000 Figure 1.6: Solution of x ( t ) vs t plot for different number of steps ( N ) inside a fulltime period T (For second-order Runge-Kutta Method) The classical fourth-order Runge-Kutta Method is given by y n +1 = y n + h (cid:20) k k k k (cid:21) (1.15)where, k = f ( x n , y n ) (1.16) k = f (cid:18) x n + h , y n + h k (cid:19) (1.17)8
10 20 30 40 50 60 70
Time (t) T o t a l ene r g y ( E ) RK-2
N=10 (a)
Time (t) T o t a l ene r g y ( E ) RK-2
N=100 (b)
Time (t) T o t a l ene r g y ( E ) RK-2
N=1000N=10000 (c)
Figure 1.7: Constant energy value of the solution over time (For second-order Runge-Kutta Method) k = f (cid:18) x n + h , y n + h k (cid:19) (1.18) k = f ( x n + h, y n + hk ) (1.19) From figure 1.8 we can see that also for very small value of N ( N = 10) the phasespace trajectory is stable and form a well ellipse.Same thing we can see in x ( t ) plots (figure 1.9). The solution is stable and wellsinusoidal.The constancy of energy value is very good in this method. We can see that a9ery negligible energy change is there (figure 1.10) which is very less compering toother methods. -1 -0.5 0 0.5 1 Position (x) -1-0.8-0.6-0.4-0.200.20.40.60.81 M o m en t u m ( p ) Phase Space Trajectory For N=10 by RK-4 Method (a) N = 10 -1 -0.5 0 0.5 1 Position (x) -1-0.500.51 M o m en t u m ( p ) Phase Space Trajectory For N=100 by RK-4 Method (b) N = 100 -1 -0.5 0 0.5 1 Position (x) -1-0.500.51 M o m en t u m ( p ) Phase Space Trajectory For N=1000 by RK-4 Method (c) N = 1000 -1.5 -1 -0.5 0 0.5 1 1.5 Position (x) -1.5-1-0.500.511.5 M o m en t u m ( p ) Phase Space Trajectory For N=10000 by RK-4 Method (d) N = 10000 Figure 1.8: Phase space trajectories for different number of steps ( N ) inside a fulltime period T (For fourth-order Runge-Kutta Method) From the plots we clearly can see that for a fixed number of steps N the solutionis better in RK-2 method than Euler method. Again RK-4 solution is much betterthan RK-2 solutions. Euler method is first order method, where RK-2 is a secondorder method. so just by increasing an order we see that result improve much. Samething is happening in the case of RK-2 and Rk-4 methods. RK-4 is a Fourth-order10
10 20 30 40 50 60 70
Time (t) -1-0.8-0.6-0.4-0.200.20.40.60.81 P o s i t i on ( x ) Particle position by RK-4 For N=10 (a) N = 10 Time (t) -1-0.500.51 P o s i t i on ( x ) Particle position by RK-4 For N=100 (b) N = 100 Time (t) -1-0.500.51 P o s i t i on ( x ) Particle position by RK-4 For N=1000 (c) N = 1000 Time (t) -1.5-1-0.500.511.5 P o s i t i on ( x ) Particle position by RK-4 For N=10000 (d) N = 10000 Figure 1.9: Solution of x ( t ) vs t plot for different number of steps ( N ) inside a fulltime period T (For Fourth-order Runge-Kutta Method)approximation whereas RK-2 is second-order. So, just by increasing two orders theerror of the solution minimizes a lot. In terms of Taylor series y ( x + h ) = y ( x ) + hy (cid:48) ( x ) + h y (cid:48)(cid:48) ( x ) + ... (1.20)The taylor series is an infinite series. computationally it is impossible to takeinfinite terms. If we just take two terms in the RHS, then that is a first-order approx-imation, which is Euler method. If we take three terms then that is a second-orderapproximation and it has been shown in section 1.5.1 that it is nothing but RK-2method. Now in figure 1.11 we plotted the relative error in energy, i.e. ε = ∆ EE where E is the total energy of the system, which should be constant of the system11
10 20 30 40 50 60 70
Time (t) T o t a l ene r g y ( E ) RK-4
N=10 (a)
Time (t) T o t a l ene r g y ( E ) RK-4
N=100N=1000N=10000 (b)
Figure 1.10: Constant energy value of the solution over time (For Fourth-order Runge-Kutta Method)and ∆ E is the energy change in the process due to accumulation of error. We plottedthis relative error with respect to number of steps N in log-log scale. Theoretically∆ E = 0 so we should get a line parallel to ε axis and passing through 0. But due toerror accumulation in each step, the actual plot that we get for different methods isshown in figure 1.11.From figure 1.11 it is clear for every process the error is less for larger number ofsteps N . The slop of the plot is however less for Euler method, RK-2 is in middleand RK-4 has maximum slop. Larger slope means the error is falling faster with N . From this assignment we conclude that obviously high order methods are less erraticand error per step is also less. But for computational simplicity we cannot choosehigher order methods as much as we want. So, we have choose a simple less ordermethod but again with less error. In next chapter we will talk about such a method,which is first order but still energy conserving.12 log N -15 -10 -5 l og Euler MethodRK-2 MethodRK-4 Method
Figure 1.11: (cid:15) vs N plot in log-log scale for different numerical methods13 hapter 2Solution by Leapfrog Method To solve the SHM problem by Leap frog method and check its energy conserving andreversible properties. To compare the solution of SHM problem by Leapfrog methodwith solutions of the same by previously discussed methods (Euler, RK-2 and RK-4).
Leapfrog method is a numerical differential equation solver of the form d xdt = a ( x ) (2.1)This second order differential equation can be written in form of a coupled firstorder differential equation as dxdt = v ( x ) and dvdt = a ( x ) (2.2)This method is very much useful for dynamical mechanical systems. This is afirst order method, so the computational complexity is less, at the same time it isless erratic process. The accumulation of errors in energy is very less, so this processis energy conserving. For our SHM problem, the force is conservative so the energyshould be constant. In previous Assignment we see that this constancy of energy isno longer retained in first order methods like Euler method. Here we will solve thesame using leapfrog method and will compare our results with previous assignment(Assignment-2). 14 .3 Algorithm of Leapfrog Method We are going to solve a second order differential equation, so to solve it we need twoinitial conditions. The initial conditions are initial position of the particle ( x ) andinitial momentum of the particle ( p ) or velocity ( v ).Figure 2.1: Pictorial representation of iteration steps in Leapfrog methodFor our SHM problem dxdt = p ( x ) /m and dpdt = − kx (2.3)From x and p , the half step of x by Leapfrog method is x = x + p m ∆ t p and x can be calculated for rest of n = 0 to n = N by p ( n + 1) = p n − kx n + ∆ t (2.5)again, x n + = x n + + p n +1 m ∆ t (2.6)There are two primary strengths to leapfrog integration when applied to mechan-ics problems. The first is the time-reversibility of the Leapfrog method. One canintegrate forward n steps, and then reverse the direction of integration and integratebackwards n steps to arrive at the same starting position. The second strength is itssymplectic nature, which implies that it conserves the (slightly modified) energy ofdynamical systems. This is especially useful when computing orbital dynamics, asmany other integration schemes, such as the (order-4) Runge–Kutta method, do notconserve energy and allow the system to drift substantially over time.15 Position (x) -1-0.500.51 M o m en t u m ( p ) Phase Space Trajectory For N=10 by Leapfrog Method (a) N = 10 -1 -0.5 0 0.5 1 Position (x) -1-0.500.51 M o m en t u m ( p ) Phase Space Trajectory For N=100 by Leapfrog Method (b) N = 100 -1.5 -1 -0.5 0 0.5 1 Position (x) -1-0.500.51 M o m en t u m ( p ) Phase Space Trajectory For N=1000 by Leapfrog Method (c) N = 1000 -1.5 -1 -0.5 0 0.5 1 Position (x) -1-0.500.51 M o m en t u m ( p ) Phase Space Trajectory For N=10000 by Leapfrog Method (d) N = 10000 Figure 2.2: Phase space trajectories for different number of steps ( N ) inside a fulltime period T (By Leapfrog method) In figure 2.2 we plot the phase-space trajectories for different number of steps ( N )inside a time period. We can see for N = 10 also we get an ellipse (slightly of axisfrom the analytical solution). The best thing is that like Euler and RK method forsmall N , like the case of N = 10 the solution is not diverging here. This clearly showswhy the Leapfrog method is reversible.In figure 2.3 we plot the x ( t ) versus t . The solution is stable, unlike the solutionsfrom Euler method and RK method, Leapfrog method gives a stable solution overtime. So even for small N the solution is not diverging. That ensures that in leapfrogmethod the error accumulation is very small.16
20 40 60
Time (t) -1-0.500.51 P o s i t i on ( x ) Particle position by Leapfrog For N=10 (a) N = 10 Time (t) -1.5-1-0.500.511.5 P o s i t i on ( x ) Particle position by Leapfrog For N=100 (b) N = 100 Time (t) -1.5-1-0.500.51 P o s i t i on ( x ) Particle position by Leapfrog For N=1000 (c) N = 1000 Time (t) -1.5-1-0.500.51 P o s i t i on ( x ) Particle position by Leapfrog For N=10000 (d) N = 10000 Figure 2.3: Solution of x ( t ) vs t plot for different number of steps ( N ) inside a fulltime period T (By Leapfrog method)Figure 2.4 shows the constancy of total energy with time. Unlike the Euler methodand RK method, here we are not getting a ever growing or ever falling energy pattern.The total energy here oscillates around a constant energy value. The average line ofthis oscillation is the constant energy value of the system. As we discussed the error in energy in Leapfrog method is negligible. In figure 2.4 wehave plotted relative error ( ε = ∆ EE ) with number of steps N in log scale. In figure 2.4the first plot (plot (a)) shows the negligible amount of energy change for a change in N . The change is so small that when we plot the relative error in Leapfrog method17
10 20 30 40 50 60 70
Time (t) T o t a l ene r g y ( E ) Leapfrog Method
N=10 (a) N = 10 Time (t) T o t a l ene r g y ( E ) Leapfrog Method
N=100 (b) N = 100 Time (t) T o t a l ene r g y ( E ) Leapfrog Method
N=1000 (c) N = 1000 Time (t) T o t a l ene r g y ( E ) Leapfrog Method
N=10000 (d) N = 10000 Figure 2.4: Constant energy value of the solution over time (By Leapfrog method).We can see the energy is oscillating around a constant value with very small amplitude.with the relative errors from other methods like Euler, RK-2 and RK-4 (figure 2.4),the Leap frog method form a horizontal line close to ε = 0 in ε vs. N plot. This plotconfirms that Leapfrog method is best energy conserving method.So, whenever in a system we need energy conserving results, Leapfrog method isthe best way.1) It is very simple to implement.2) It is first order so computational complexity is less.3) It gives stable solution.4) The error accumulation in each step is negligible.5) The solutions are reversible. 18 log N l og Leapfrog Method (a) For Leapfrog only log N -15 -10 -5 l og Euler MethodRK-2 MethodRK-4 MethodLeapfrog Method (b) Comparison with other methods
Figure 2.5: (cid:15) vs N plot in log-log scale for different numerical methods6) The solutions are energy conserving.7) For small number of steps N the solution is not diverging.19 hapter 3Solution of Poisson’s Equationusing Finite Difference Method To discuss the Finite Difference Method (FDM) and solve Poisson’s equation usingFDM.
Our goal is to solve Poisson equation. In previous assignments we learnt how to solvePoisson equation using Fourier method. In this assignment we will solve Poissonequation by more direct approach, Finite difference method.
We want to solve ∇ φ = 4 πGρ (3.1)where ρ is given and we want to solve for φ . In 1D the Poisson equation is d φdx = 4 πGρ (3.2)Now, by Taylor series expansion we have φ ( x + h ) = φ ( x ) + h φ (cid:48) ( x ) + h φ (cid:48)(cid:48) ( x ) + ... (3.3)20 ( x − h ) = φ ( x ) − h φ (cid:48) ( x ) + h φ (cid:48)(cid:48) ( x ) − ... (3.4)By eq (3 . − eq (3 .
4) we have φ (cid:48) ( x ) = φ ( x + h ) − φ ( x − h )2 h + O ( h ) h (3.5)and by eq (3 .
3) + eq (3 .
4) we have φ (cid:48)(cid:48) ( x ) = φ ( x + h ) + φ ( x − h ) − φ ( x ) h + O ( h ) h (3.6)Eq. 3.5 and eq. 3.6 shows the discrete version of first and second derivative of afunction respectively. In the discrete form eq 3.2 can be written as φ ( x + h ) + φ ( x − h ) − φ ( x ) h = 4 πGρ ( x ) + O ( h ) h (3.7)In discrete space (computational space) x is replaced by discrete integer a . φ a +1 + φ a − − φ a L = 4 πGρ a (3.8)neglecting higher order terms, L is the grid spacing. We done the simulation for a 2D grid, the parameters are:1) square grid side length = 1002) grid spacing L = 13) The value of potential φ at boundaries of the square grid is assumed to be zero(As with distance the potential value falls down, and for infinite distance it is zero.However we can not work with infinity, so it is a good approximation to set potentialvalues zero at boundaries of the finite square grid over that we are simulating). Figure 3.1 to figure 3.4 shows potential solutions for different source mass configura-tions. (For getting a clear plot I just reverse the sign of the potential before plotting).21 a) 3D surface plot x y Color-plot for 2D potential solution (b) Color plot
Figure 3.1: Potential ( φ ) plot for a point mass at the centre (50 ,
50) of 2D grid(100 × (a) 3D surface plot x y Color-plot for 2D potential solution (b) Color plot
Figure 3.2: Potential ( φ ) plot for two point masses at (33 ,
66) and (66 ,
33) of a(100 × a) 3D surface plot x y Color-plot for 2D potential solution (b) Color plot
Figure 3.3: Potential ( φ ) plot for four point masses at (33 , , ,
33) and(66 ,
66) of a (100 × (a) 3D surface plot x y Color-plot for 2D potential solution (b) Color plot
Figure 3.4: Potential ( φ ) plot for a square 2D uniform bulk mass distributed between x = 33, x = 66, y = 33 and y = 66 of a (100 × hapter 4Fourier Method Solution ofPoisson’s Equation To find the potential-field and acceleration-field created by a known mass densitydistribution ( ρ ) for 1D system by Fourier transform and inverse Fourier transformmethod. We can write any function in Fourier series expansion, as the terms in Fourier seriesforms a complete basis set. Same is true for Fourier transform. Fourier transformforms a complete continuous basis to expand any function.We want to solve the potential and finally acceleration field created by a massdensity distribution ( ρ ). So, the source mass distribution is ρ , which is known asfunction of distance.We cannot deal with continuous variables in computation. The distance x is justa number a in discrete system. So we have ρ a . By discrete Fourier transform of thiswe have ˜ ρ n = 1 N N − (cid:88) a =0 ρ a exp (cid:20) − πinaN (cid:21) (4.1)From this mass density in Fourier space we can calculate the potential in Fourierspace via the Formula derived in the class:24 φ n = − (4 πG ) ˜ ρ n L πnN for n > A n = − i sin (cid:18) πnN (cid:19) ˜ φ n L (4.3)Now, we have ˜ φ n and ˜ A n , so, we can just take inverse Fourier transform to calculateour requre potential field φ a and acceleration field A a in real space: φ a = N − (cid:88) n =0 ˜ φ n exp (cid:20) πinaN (cid:21) (4.4)and, A a = N − (cid:88) n =0 ˜ A n exp (cid:20) πinaN (cid:21) (4.5) In this assignment we take a 1D grid of grid size N = 128. We work with three typesof mass density ρ a . In figure 4.1 and figure 4.2 we have a point mass (source) at a = 64. For thismass configuration these figure shows the plots for all different parameters in realand Fourier space. We plotted real and imaginary parts separately for all complexparameters.In figure 4.3 and figure 4.4 the source is two point masses situated at a = 40 and a = 80.In figure 4.5 and figure 4.6 the source is three point masses situated at a = 20, a = 40 and a = 60. 25
50 100 150 a a=64 (a) n -1-0.500.51 a=64 (b) a -10-505 R e () -9 a=64 (c) a -505 R e ( A ) -10 a=64 (d) n -1-0.500.511.522.533.5 10 -7 a=64 (e) n -1-0.500.51 10 -9 a=64 (f) Figure 4.1: Real parts of all different complex parameters (a point mass at a = 64 isthe source mass) 26
20 40 60 80 100 120 140 n -1-0.8-0.6-0.4-0.200.20.40.60.81 a=64 (a) a -101 I m () -23 a=64 (b) a -6-4-20246 I m ( A ) -25 a=64 (c) n -2-1.5-1-0.500.511.52 10 -8 a=64 (d) n -2-1.5-1-0.500.511.52 10 -8 a=64 (e) Figure 4.2: Imaginary parts of all different complex parameters (a point mass at a = 64 is the source mass) 27
20 40 60 80 100 120 140 a a=40 and a=80 (a) n -2-1.5-1-0.500.511.52 a=40 and a=80 (b) a -8-6-4-202468 R e () -9 a=40 and a=80 (c) a -6-4-20246 R e ( A ) -10 a=40 and a=80 (d) n -1-0.500.511.522.533.54 10 -7 a=40 and a=80 (e) n -8-6-4-20246 10 -9 a=40 and a=80 (f) Figure 4.3: Real parts of all different complex parameters (two point masses at a = 40and a = 80 are the source distribution) 28
20 40 60 80 100 120 140 n -2-1.5-1-0.500.511.52 a=40,80 (a) a -2-1012 I m () -23 a=40,80 (b) a -6-4-202468 I m ( A ) -25 a=40,80 (c) n -1-0.500.51 10 -7 a=40,80 (d) n -2-1.5-1-0.500.511.52 10 -8 a=40,80 (e) Figure 4.4: Imaginary parts of all different complex parameters (two point masses at a = 40 and a = 80 are the source distribution)29
20 40 60 80 100 120 140 a a=20,40,60 (a) n -3-2-10123 a=20,40,60 (b) a -1.5-1-0.500.511.5 R e () -8 a=20,40,60 (c) a -1-0.500.51 R e ( A ) -9 a=20,40,60 (d) n -0.500.511.522.5 10 -7 a=20,40,60 (e) n -1-0.500.511.522.533.5 10 -8 a=20,40,60 (f) Figure 4.5: Real parts of all different complex parameters (three point masses at a = 20, a = 40 and a = 60 are the source distribution)30
50 100 150 n -3-2-10123 a=20,40,60 (a) a -3-2-10123 I m () -23 a=20,40,60 (b) a -6-4-20246 I m ( A ) -25 a=20,40,60 (c) n -8-6-4-202468 10 -7 a=20,40,60 (d) n -1.5-1-0.500.511.5 10 -8 a=20,40,60 (e) Figure 4.6: Imaginary parts of all different complex parameters (three point massesat a = 20, a = 40 and a = 60 are the source distribution)31 hapter 5Solution of Poisson’s Equation byCIC Method To find the acceleration value at test particle position, where in generally test particlemay not be situated on a grid position.
In Assignment-4 we solved the acceleration and potential field for a given mass density.However the acceleration and potential field that we calculated there is on the grid.The density of particles also was distributed on the grid points only. In real physicalproblem the source particle and test particle positions may not be always on a gridpoint. In generally it can be anywhere between the grid points.In that case we have to use Cloud in Cell (CIC) method.
If the position for any source particle between two grid points is x then according toCIC method the mass of the particle will distribute among the grid points. The massdistribution will be weighted average. The closer a mass particle to a grid point, moreweight will be given for that grid point.For a 1D case if the source particle position is x and it is situated between thegrid points a = f loor and a = ceil . If the mass of the particle is m then32 = 1 − | x − f loor | and w = 1 − | ceil − x | (5.1)and, ρ floor = w × m and ρ ceil = w × m (5.2)Now, as we have density of mass on the grid the calculation of potential field andacceleration on the grid is as the previous assignment (Assignment-4).If the test particle is at x test and situated between the grid points a = f t and a = ct then acceleration at x test is A ( x test ) = w × A ft + w × A ct (5.3)where, w = 1 − | x − f t | and w = 1 − | ct − x | (5.4) We simulate for a 1D grid of L = 128.The source is an unit point mass at x = 64 . x test = 10 . The value of the acceleration at the test particle position ( x test = 10 .
8) is A ( x test ) =6 . × − N/kg
For this system we plot the real and imaginary solutions of different real andFourier space parameters in the figure 5.1 (real parts of the solutions) and figure 5.2(imaginary parts of the solutions). 33
50 100 150 a a Source =64.4; a tp =10.8 (a) n -1-0.500.51 a Source =64.4; a tp =10.8 (b) a -10-505 R e () -9 a Source =64.4; a tp =10.8 (c) a -505 R e ( A ) -10 a Source =64.4; a tp =10.8 (d) n -101234 10 -7 a Source =64.4; a tp =10.8 (e) n -6-4-20246 10 -10 a Source =64.4; a tp =10.8 (f) Figure 5.1: Real parts of all different complex parameters (The source particle is at x = 64 . x test = 10 .
20 40 60 80 100 120 140 n -0.6-0.4-0.200.20.40.6 a Source =64.4; a tp =10.8 (a) a -101 I m () -23 a Source =64.4; a tp =10.8 (b) a -4-3-2-1012345 I m ( A ) -25 a Source =64.4; a tp =10.8 (c) n -1.5-1-0.500.511.5 10 -8 a Source =64.4; a tp =10.8 (d) n -2-1.5-1-0.500.511.52 10 -8 a Source =64.4; a tp =10.8 (e) Figure 5.2: Imaginary parts of all different complex parameters (The source particleis at x = 64 . x test = 10 . hapter 6Simulation of 1D Self-GravitatingSystem To simulate a 1D distribution of particle under the influence of its own gravity. Todraw the phase space snapshots over time.
After doing all the previous assignments we now know how to solve differential equa-tion, we know that by Leapfrog method the solution is energy conserving. Fromassignment-4 we know how to solve Poisson’s equation using Fourier method. Fromassignment-5 we know how to even solve for acceleration and potential when the testand source particles are not on grid (by CIC method). Now finally we are going toaddress our goal, to solve a self gravitating system. In this assignment we solve the1D self gravtating system.
We normalize the time and space using appropriate scaling parameters. We normalizetime using t g = a √ Gρ . This t g is the gravitational time scale, which is the time takenby a uniformly distributed spherical body to collapse. The length x is normalizedusing the grid spacing L as a scaling length.So, dimensionless X = x/L and dimensionless T = t/t g dXdT = t g PLm and dPdT = t g fLm = F (dimensionless force) (6.1)The form of dimensionless force F can be calculated as F = 4 π (cid:20) −∇ x . ( ∇ − x ) ρ ¯ ρ (cid:21) (6.2)¯ ρ = 1 for our simulation, but ρ will evolve as the system evolve. Simulation parameters are: L = 100; total space in 1D dx = 1; length of each segment N = S/dx ; number of grid points T = 100; total run time dt = 0 .
1; time stepAll the variables are normalized in this simulation.
All the plots in figure 6.1 to figure 6.4 are the phase space snapshots for time t = 10to t = 800 with an interval of t = 10. We start with particles uniformly distributedover a 1D line with zero momentum. As time pass particles will get momentum dueto gravitational attraction and try to collapse (cluster) at the centre. We use periodicboundary condition here. A video of the simulation done in this assignment can be found in this link https://youtu.be/tzBqdSGEpTQ . However with this assignment PDF file the video is alsoattached. 37
20 40 60 80 100 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (a) t = 10 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (b) t = 20 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (c) t = 30 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (d) t = 40 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (e) t = 50 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (f) t = 60 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (g) t = 70 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (h) t = 80 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (i) t = 90 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (j) t = 100 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (k) t = 110 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (l) t = 120 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (m) t = 130 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (n) t = 140 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (o) t = 150 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (p) t = 160 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (q) t = 170 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (r) t = 180 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (s) t = 190 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (t) t = 200 Figure 6.1: Phase Space snapshots for t = 10 to t = 200 with time separation of 1038
20 40 60 80 100 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (a) t = 210 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (b) t = 220 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (c) t = 230 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (d) t = 240 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (e) t = 250 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (f) t = 260 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (g) t = 270 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (h) t = 280 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (i) t = 290 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (j) t = 300 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (k) t = 310 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (l) t = 320 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (m) t = 330 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (n) t = 340 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (o) t = 350 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (p) t = 360 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (q) t = 370 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (r) t = 380 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (s) t = 390 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (t) t = 400 Figure 6.2: Phase Space snapshots for t = 210 to t = 400 with time separation of 1039
20 40 60 80 100 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (a) t = 410 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (b) t = 420 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (c) t = 430 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (d) t = 440 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (e) t = 450 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (f) t = 460 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (g) t = 470 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (h) t = 480 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (i) t = 490 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (j) t = 500 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (k) t = 510 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (l) t = 520 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (m) t = 530 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (n) t = 540 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (o) t = 550 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (p) t = 560 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (q) t = 570 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (r) t = 580 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (s) t = 590 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (t) t = 600 Figure 6.3: Phase Space snapshots for t = 410 to t = 600 with time separation of 1040
20 40 60 80 100 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (a) t = 610 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (b) t = 620 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (c) t = 630 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (d) t = 640 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (e) t = 650 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (f) t = 660 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (g) t = 670 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (h) t = 680 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (i) t = 690 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (j) t = 700 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (k) t = 710 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (l) t = 720 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (m) t = 730 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (n) t = 740 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (o) t = 750 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (p) t = 760 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (q) t = 770 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (r) t = 780 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (s) t = 790 position -25-20-15-10-50510152025 m o m en t u m Phase-space-plot (t) t = 800 Figure 6.4: Phase Space snapshots for t = 610 to t = 800 with time separation of 1041 ibliography Anastassi, Z. & Simos, T. (2005), ‘An optimized runge–kutta method for the solutionof orbital problems’,
Journal of Computational and Applied Mathematics (1), 1–9.Atkinson, K. (1985), ‘The numerical evaluation of particular solutions for poisson’sequation’.Bhat, P., Curless, B., Cohen, M. & Zitnick, C. L. (2008), Fourier analysis of the 2dscreened poisson equation for gradient domain problems, in ‘European Conferenceon Computer Vision’, Springer, pp. 114–128.Butcher, J. C. (1976), ‘On the implementation of implicit runge-kutta methods’, BITNumerical Mathematics (3), 237–240.Cash, J. R. & Karp, A. H. (1990), ‘A variable order runge-kutta method for ini-tial value problems with rapidly varying right-hand sides’, ACM Transactions onMathematical Software (TOMS) (3), 201–222.Gel’fand, I. M. & Dikii, L. A. (1979), ‘Integrable nonlinear equations and the liouvilletheorem’, Functional Analysis and its Applications (1), 6–15.Ghosh, T., Pramanick, S., Sarkar, S., Dey, A. & Chandra, S. (2020), ‘Dynamicalproperties and effects of quantum diffraction on the propagation of e-a-solitarywaves in three-component fermi plasma’, arXiv:2012.13616v1 p. 8.Greiner, A., Strittmatter, W. & Honerkamp, J. (1988), ‘Numerical integration ofstochastic differential equations’, Journal of Statistical Physics (1-2), 95–108.Hahn, G. (1991), ‘A modified euler method for dynamic analyses’, International Jour-nal for Numerical Methods in Engineering (5), 943–955.42serles, A. (1986), ‘Generalized leapfrog methods’, IMA Journal of Numerical Analysis (4), 381–392.Jordan, T. F. (2004), ‘Steppingstones in hamiltonian dynamics’, American journal ofphysics (8), 1095–1099.Liu, T., Zhao, C., Li, Q. & Zhang, L. (2012), ‘An efficient backward euler time-integration method for nonlinear dynamic analysis of structures’, Computers &structures , 20–28.Lyashko, A. & Solov’Yev, S. (1991), ‘Fourier method of solution of fe systems withhermite elements for poisson equation’,
Russian Journal of Numerical Analysis andMathematical Modelling (2), 121–130.Matthews, A. P. (1994), ‘Current advance method and cyclic leapfrog for 2d multi-species hybrid plasma simulations’, Journal of Computational Physics (1), 102–116.Merle, F. & Zaag, H. (2000), ‘A liouville theorem for vector-valued nonlinear heatequations and applications’,
Mathematische Annalen (1), 103–137.Mikkola, S. & Aarseth, S. (2002), ‘A time-transformed leapfrog scheme’,
CelestialMechanics and Dynamical Astronomy (4), 343–354.Pramanick, S. (2020), ‘Onsets and outflow distributions in abelian and stochastic btwmodels’, arXiv preprint arXiv:2009.04516 .Pramanick, S., Dey, A. & Chandra, S. (2020), ‘Electron-acoustic solitary waves infermi plasma with two-temperature electrons’, DOI: 10.2139/ssrn.3711910 .Rice, J. R. (1960), ‘Split runge-kutta method for simultaneous’,
Journal of Researchof the National Bureau of Standards: Mathematics and mathematical physics. B , 151.Sahana, S., Mitra, S., Chandra, S. & Pramanick, S. (2021), ‘Linear properties of elec-trostatic wave in two-component fermi plasma’, arXiv preprint arXiv:2101.00510 . 43hampine, L. F. (2009), ‘Stability of the leapfrog/midpoint method’, Applied mathe-matics and computation (1), 293–298.Skeel, R. D. (1993), ‘Variable step size destabilizes the st¨ormer/leapfrog/verletmethod’,
BIT Numerical Mathematics (1), 172–175.Sk¨ollermo, G. (1975), ‘A fourier method for the numerical solution of poisson’s equa-tion’, Mathematics of Computation (131), 697–711.Xiong, X., Chen, W., Jin, S. & Kamal, S. (2019), ‘Discrete-time implementation ofcontinuous terminal algorithm with implicit-euler method’, IEEE Access7