I am reading the article "First few overtones probe the event horizon geometry". I want to replicate the results of Table 1 in the above article using the analytical continuation method as mentioned in the article. However, my code doesn't run out of any results. I don't know what's wrong with it, because there is no warning message.
Any suggestions would be greatly appreciated. Thanks very much.
The code is as follows:
M = 1/2; x = 1 - r0/r; \[Epsilon] = (2 M - r0)/r0; a0 = b0 = 0;
A[r_] = 1 - \[Epsilon] (1 - x) + (a0 - \[Epsilon]) (1 - x)^2 +
Al[r] (1 - x)^3;
B[r_] = 1 + b0 (1 - x) + Bl[r] (1 - x)^2;
Nm[r_] = Sqrt[x A[r]];
Al[r_] = 1/x ContinuedFractionK[a[n] x, 1, {n, 3}];
Bl[r_] = 1/x ContinuedFractionK[b[n] x, 1, {n, 3}];
f[r_] = Nm[r]^2/B[r];
V[r_] = Nm[r]^2 l (l + 1)/r^2 + (1 - s)/2/r D[f[r]^2, r];
eq = f[r] D[f[r] D[R[r], r], r] + (\[Omega]^2 - V[r]) R[r];
R[r_] = Exp[\[Sigma] r] x^\[Rho] (x - 1)^\[Lambda] y[r];
a[1] = 1/2; a[2] = 100; a[3] = 0; b[1] = 0; \[Epsilon] = 0; r0 = 1;
fctor = Coefficient[R[r], y[r]];
eqr = eq/fctor // Simplify;
eqrr = Collect[
eqr, { y[r], Derivative[1][y][r], (y^\[Prime]\[Prime])[r]}];
epeq = Series[eqrr, {r, Infinity, 3}];
epf = Solve[
Coefficient[Coefficient[epeq, y[r]], r, 0] == 0 &&
Coefficient[Coefficient[epeq, y[r]], r, -1] ==
0, {\[Sigma], \[Lambda]}];
(*solve the parameter {\[Sigma],\[Lambda]} using the boundary \
condition at infinity*)
eqnr = Collect[eqr, (-1 + r)] // Together // Simplify;
eqrho = eqnr /. r -> 1;
rho = Solve[Coefficient[eqrho, y[1]] == 0, \[Rho]];
(*solve the parameter \[Rho] using the boundary condition at the \
event horizon*)
Equation = Simplify[eqr /. rho[[1]] /. epf[[2]]];
(*the radial function with rho[[1]] and epf[[2]] will satisfy the \
quasinormal boundary conditions*)
r[z_] = First[r /. Solve[z == 1 - 1/r, r]];
CompactCoordinateEquation =
Equation /. (y^\[Prime]\[Prime])[r] ->
D[D[y[z], z]/D[r[z], z], z]/D[r[z], z] /.
Derivative[1][y][r] -> D[y[z], z]/D[r[z], z] /. y[r] -> y[z] /.
r -> r[z] // Expand;
y2Coeff = Coefficient[CompactCoordinateEquation, y''[z]] // Together;
y1Coeff = Coefficient[CompactCoordinateEquation, y'[z]] // Together;
yCoeff = Coefficient[CompactCoordinateEquation, y[z]] // Together;
fr = PolynomialGCD[y2Coeff, y1Coeff, yCoeff];
y2RCoeff = Simplify[y2Coeff/fr]; y1RCoeff =
Simplify[y1Coeff/fr]; yRCoeff = Simplify[yCoeff/fr];
PolynomialQ [y2RCoeff, z] && PolynomialQ [y1RCoeff, z] &&
PolynomialQ [yRCoeff, z];
CoefficientListMatrix =
CoefficientList[{y2RCoeff/z, y1RCoeff, z yRCoeff}, z];
z1 = 15/1000; z2 = 30/1000; z3 = 60/1000; z4 = 120/1000; z5 =
240/1000; z6 = 480/1000; z7 = 960/1000;
(*midpoints for continuation*)
y2RCoeffst = y2RCoeff /. z -> t + z1;
y1RCoeffst = y1RCoeff /. z -> t + z1;
yRCoeffst = yRCoeff /. z -> t + z1;
CoefficientListMatrix1 =
CoefficientList[{y2RCoeffst, t y1RCoeffst, t^2 yRCoeffst}, t];
y2RCoeffnd = y2RCoeff /. z -> t + z2;
y1RCoeffnd = y1RCoeff /. z -> t + z2;
yRCoeffnd = yRCoeff /. z -> t + z2;
CoefficientListMatrix2 =
CoefficientList[{y2RCoeffnd, t y1RCoeffnd, t^2 yRCoeffnd}, t];
y2RCoeffrd = y2RCoeff /. z -> t + z3;
y1RCoeffrd = y1RCoeff /. z -> t + z3;
yRCoeffrd = yRCoeff /. z -> t + z3;
CoefficientListMatrix3 =
CoefficientList[{y2RCoeffrd, t y1RCoeffrd, t^2 yRCoeffrd}, t];
y2RCoefffr = y2RCoeff /. z -> t + z4;
y1RCoefffr = y1RCoeff /. z -> t + z4;
yRCoefffr = yRCoeff /. z -> t + z4;
CoefficientListMatrix4 =
CoefficientList[{y2RCoefffr, t y1RCoefffr, t^2 yRCoefffr}, t];
y2RCoefffv = y2RCoeff /. z -> t + z5;
y1RCoefffv = y1RCoeff /. z -> t + z5;
yRCoefffv = yRCoeff /. z -> t + z5;
CoefficientListMatrix5 =
CoefficientList[{y2RCoefffv, t y1RCoefffv, t^2 yRCoefffv}, t];
y2RCoeffsx = y2RCoeff /. z -> t + z6;
y1RCoeffsx = y1RCoeff /. z -> t + z6;
yRCoeffsx = yRCoeff /. z -> t + z6;
CoefficientListMatrix6 =
CoefficientList[{y2RCoeffsx, t y1RCoeffsx, t^2 yRCoeffsx}, t];
y2RCoeffsh = y2RCoeff /. z -> t + z7;
y1RCoeffsh = y1RCoeff /. z -> t + z7;
yRCoeffsh = yRCoeff /. z -> t + z7;
CoefficientListMatrix7 =
CoefficientList[{y2RCoeffsh, t y1RCoeffsh, t^2 yRCoeffsh}, t];
terms = CoefficientListMatrix // Dimensions // Last; rcf =
Function[{i,
n}, (n + 2 - i) (n + 1 - i) CoefficientListMatrix[[1,
i]] + (n + 2 - i) CoefficientListMatrix[[2, i]] +
CoefficientListMatrix[[3, i]]];
terms1 = CoefficientListMatrix1 // Dimensions // Last; rcf1 =
Function[{i,
n}, (n + 3 - i) (n + 2 - i) CoefficientListMatrix1[[1,
i]] + (n + 3 - i) CoefficientListMatrix1[[2, i]] +
CoefficientListMatrix1[[3, i]]];
terms2 = CoefficientListMatrix2 // Dimensions // Last; rcf2 =
Function[{i,
n}, (n + 3 - i) (n + 2 - i) CoefficientListMatrix2[[1,
i]] + (n + 3 - i) CoefficientListMatrix2[[2, i]] +
CoefficientListMatrix2[[3, i]]];
terms3 = CoefficientListMatrix3 // Dimensions // Last; rcf3 =
Function[{i,
n}, (n + 3 - i) (n + 2 - i) CoefficientListMatrix3[[1,
i]] + (n + 3 - i) CoefficientListMatrix3[[2, i]] +
CoefficientListMatrix3[[3, i]]];
terms4 = CoefficientListMatrix4 // Dimensions // Last; rcf4 =
Function[{i,
n}, (n + 3 - i) (n + 2 - i) CoefficientListMatrix4[[1,
i]] + (n + 3 - i) CoefficientListMatrix4[[2, i]] +
CoefficientListMatrix4[[3, i]]];
terms5 = CoefficientListMatrix5 // Dimensions // Last; rcf5 =
Function[{i,
n}, (n + 3 - i) (n + 2 - i) CoefficientListMatrix5[[1,
i]] + (n + 3 - i) CoefficientListMatrix5[[2, i]] +
CoefficientListMatrix5[[3, i]]];
terms6 = CoefficientListMatrix6 // Dimensions // Last; rcf6 =
Function[{i,
n}, (n + 3 - i) (n + 2 - i) CoefficientListMatrix6[[1,
i]] + (n + 3 - i) CoefficientListMatrix6[[2, i]] +
CoefficientListMatrix6[[3, i]]];
terms7 = CoefficientListMatrix7 // Dimensions // Last; rcf7 =
Function[{i,
n}, (n + 3 - i) (n + 2 - i) CoefficientListMatrix7[[1,
i]] + (n + 3 - i) CoefficientListMatrix7[[2, i]] +
CoefficientListMatrix7[[3, i]]];
NITMAX = 100;
\[Alpha] = rcf[#1, #2] &;
fs = Table[
If[i > j + 2, 0, \[Alpha][i, j]], {j, 0, NITMAX}, {i, 1, 10}];
Do[fs = FoldList[
If[PossibleZeroQ[#2 // Last], #2 //
Most, (Most[#2] - (#2 // Last)/(#1 // Last)*
Prepend[#1 // Most, 0]) // Together] &, fs // First // Most,
fs // Rest], {7}]
(*Gauss elimination*)
Attachments:
