闪电代写 -代写CS作业_CS代写_Finance代写_Economic代写_Statistics代写_代码代做_IT代写_加急帮助

Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit

Physics 305 – Computational Physics, Fall 2021

Term Project

The program in your term project can be either submitted as a python program or ipython notebook, where the latter is preferred. The program, an explanation of what the program does, along with answers to all questions asked should be uploaded to d2l.

You are expected to write a term paper (in word or Latex) on your project that discusses the problem you are trying to solve, the basic equations that govern the problem, includes plots that show the solutions, and describes the solution and the numerical method involved. In addition, you must demonstrate that your solution is correct by showing that the code converges at the expected order. If your code does not converge at the expected order you should try to identify potential reasons for why this is the case. You are expected to work on your term project by yourself. .

Your term project will receive full credit only if: (a) the program runs successfully without errors using python 3, (b) the programs have explanatory comments and variable names that identify with the problem equations you are trying to solve, (c) give the correct output, and (d) demonstrate the validity of the solution through convergence plots. No credit will be given to late term projects.

The term paper is as important as the code (50% of the term project credit will go to the code and the other 50% to the paper). Answers to the questions and analysis requested below should be elaborated in the report. Plots should be clearly labeled and be properly described in the report, and not just shown. You will need to explain what each and every plot demonstrates. A polished paper written in word or LaTex (preferred, e.g. please try overleaf) is expected to get full credit.

Note: Before you present results from numerical integrations that answer the questions in the project, it is critical to *ﬁrst* perform the convergence tests for one case, and to estimate errors. This will tell you how small a step size is necessary for accurate solutions. Only after errors are estimated, does it make sense to run your code for producing results that help you learn more about the system you study.

I. EQUATORIAL KERR GEODESICS: ORBITS AROUND SPINNING BLACK HOLES

This project is about understanding numerically orbits around a spinning black hole for massive and massless particles. For simplicity we will focus on orbits restricted to the orbital plane of the black hole.

According to the theory of general relativity, the equations of motion (EOM) on the equatorial plane of a single spinning black hole in polar coordinates are

dγ ,[E(γ2 + a2 ) - La]2 - ∆[u2 γ 2 + (L - aE)2]

d≠ -a(aE - L) + (γ2 + a2 )(E(γ2 + a2 ) - La)/∆

d。 -(aE - L) + a[E(γ2 + a2 ) - La]/∆

d≠ -a(aE - L) + (γ2 + a2 )(E(γ2 + a2 ) - La)/∆ (1)

where E , L are the conserved energy and angular momentum of the particle, u = 1 for a massive particle and u = 0 for a massless particle, e.g., a photon. In addition, a is a parameter that controls the angular momentum of the black hole. In particular the angular momentum of the black hole is J = aM , where M is the mass of the black hole. The ratio J/M2 is dimensionless (we are adopting units in which G = c = 1), and for a black hole it must be -1 ≤ J/M2 ≤ 1. Thus, -1 ≤ a/M ≤ 1. We also have

∆ = γ 2 - 2M γ + a2 (2)

Note that the black hole has a horizon located at

γ五 = M +,M2 - a2 , (3)

which becomes singular in these coordinates. In addition, another singularity appears at the outer edge of the black hole ergosphere γ = 2M , that we mention below. When integrating the equations always terminate the integration if (γ - γ五 )/γ五 ≤ 0.02. Essentially if you reach distances that are so close to the horizon the particle will fall into the BH.

You will be considering circular orbits, and parabolic orbits. Parabolic ortbis are marginally unbound and cor- respond to zero binding energy. The energy E contains the rest-mass energy of the particle, thus parabolic orbits correspond to E = 1. Bound orbits (circular and elliptic) correspond to energy E < 1. So, the energy and angular momenta for the orbits are as follows:

$ Parabolic orbits: The energy is E = 1. The angular momentum L for parabolic orbits can be determined by the pericenter distance (closest approach) γ夕 at which dγ/d≠ = 0, which from the EOM implies

[(γ夕(2) + a2 ) - La]2 - ∆[u2 γ夕(2) + (L - a)2] = 0 (4)

Thus, by choosing γ夕 one can solve this last equation to determine L

2aM 土 ′γ夕 ╱a2 - 2M γ夕 + γ夕(2)、(2u2 M - (u2 - 1) γ夕 ) (5)

2M - γ夕

In other words if you choose the pericenter distance γ夕 , you immediately know what the constant E and L should go into determining the EOM for parabolic orbits. Notice that at the edge of the ergoshere L blows up. Thus, do not set a parabolic orbit with pericenter distance of γ夕 = 2M .

$ A somewhat involved algebraic calculation yields that for circular orbits the energy and angular momentum are given by

γ 3/2 - 2M γ1/2 土 aM1/2

γ 3/4(γ3/2 - 3M γ1/2 土 2aM1/2)1/2

土M1/2(γ2 ￥ 2aM1/2γ 1/2 + a2 )

γ 3/4(γ3/2 - 3M γ1/2 土 2aM1/2)1/2 (6)

In other words if you choose the radius of the circular orbit γ, you immediately know what the constant E and L should go into determining the EOM.

In these equations the upper signs apply to prograde orbits (i.e., particle angular momentum is aligned with the black hole angular momentum), while the lower signs apply to retrograde motion.

For all calculations that follow we will assume that the black hole has unit mass M = 1 and spin parameter a = 0.46. Strictly speaking the mass of the black hole sets the length and time scales in the problem, and when solving the EOM for M = 1 we can always ﬁnd the solution for any M by rescaling the length and time scales. Show this for yourself: introduce a new dimensionless time coordinate = ≠M and a new dimensionless radial coordinate = γM , and rewrite the above equations in dimensionless form. Show, that the form of the dimensionless equations is the same as those with dimension but with M = 1.

The energy of the

where

particle

is also equal to

E = -夕 (g芒芒芒 - g芒φ ),

(7)

芒 -a(aE - L) + (γ2 + a2 )[E(γ2 + a2 ) - La]/∆

g芒芒 = - ╱ 1 - 、 (8)

2Ma

The expression (7) can serve as a useful diagnostic for the validity of the numerical integration. You should ensure that at ≠ = 0 Eq. (7) yields the same value as the energy E that you input for the numerical integration. You will need to plot the orbits to determine at what the integration should be terminated. For example for circular orbits stop after you have completed a few periods or the particle reaches the ergoregion and falls into the black hole. For parabolic orbits stop when the particle goes back to the large distance it starts out or falls into the black hole. For the purpose of ploting the orbits, after having γ and 。the cartesian coordinates can be obtained by

α = ,(γ2 + a2 ) cos(。), 夕 = ,(γ2 + a2 ) sin(。). (9)

1. Use RK4, to integrate numerically the dimensionless form of system of ODEs (1) with initial conditions = 0 and 。= 0, i.e., initially on the x-axis. .

$ Treating only massive particles u = 1, for parabolic orbits consider (≠ = 0) = 100. For prograde orbits consider a range of pericenter distances 夕 = 3.5, 3.1, 2.9, 2.5. What do you notice as the pericenter distance becomes smaller than 3.0? Make a plot of the orbits. For retrograde orbits consider a range of pericenter distances 夕 = 5, 4.9, 4.7, 4.5. What do you notice as the pericenter distance becomes smaller than 4.9? Make a plot of the orbits.

$ For circular orbits consider both massive (u = 1 and massless (u = 0) particles. For u = 1 and prograde orbits consider a range of radii = 5.0, 4.5, 4.3, 4.0. What do you notice as the radius becomes smaller than 4.5? For u = 1 and retrograde orbits consider a range of radii = 8.0, 7.5, 7.0. What do you notice as the radius becomes smaller than 7.5? Make a plot of the orbits. For massless particles consider prograde orbits only. Consider a range of radii 夕 = 3.5, 3.1, 2.5, 2.0. What do you notice is diﬀerent with respect to the program massive particle orbits? Make a plot of the orbits.

Use your judgement as to how small a step size you need to solve this system accurately. If you cannot ﬁgure this out from pure thought, experiment with diﬀerent step sizes. Use 6E = |(Etot - Etot (≠ = 0))/Etot (≠ = 0)| to determine the accuracy. If 6E is smaller than 10 −3 for all integration times, then you have decent accuracy.

2. Convergence: Demonstrate for one stable circular orbit that as you decrease the step size 左, 6E at the ﬁnal time you integrate converges to 0. Make a plot to determine the order of convergence of 6E, and discuss why or why not this matches your expectation.

3. Self-convergence: Use a number of step sizes and make a plot to demonstrate that the code solution for α and 夕 at the ﬁnal time self-converges.

4. Using the order of convergence you determined, employ Richardson extrapolation to determine an error for the solution for α(≠), 夕(≠) at a time ≠ of your choosing.