I’m trying to solve a system of 24 non-linear equations, using Implicit Runge Kutta method. Code is working very smoothly with MachinePrecision. But, I want higher precision as much as I can. First, I change all initial conditions upto 34 decimal places, for example, using EE = N(976/1000, 40) for energy = EE = 0.976, and same for other initial conditions.
I’m calculating initial data using FindRoot, where I’m using WorkingPrecision -> 39, then I have precision 34 of some variables, and some have 36, etc.
I don’t know how to convert all initial conditions in the same precision, therefore I used like this.
Is there any command that we can use that convert all initial data in the same precision?
Moreover, I want to solve this system with precision higher than the MachinePrecision. I have tried using precision AccuracyGoal -> 18, PrecisionGoal -> 18, etc or using WorkingPrecision -> 40, but code gives errors.
Can anyone please help me to solve the system with precision higher than the MachinePrecision or quadruple-precision .
I’m posting my code here
n = 4;
AA(r_) := (1 - (2 M)/r); M = 1;
gtt(r_, θ_) := -AA(r); grr(r_, θ_) := 1/AA(r);
gθθ(r_, θ_) := r^2;
gϕϕ(r_, θ_) := r^2 Sin(θ)^2;
gUtt(r_, θ_) := 1/gtt(r, θ);
gUrr(r_, θ_) := 1/grr(r, θ);
gUθθ(r_, θ_) := 1/gθθ(r, θ);
gUϕϕ(r_, θ_) := 1/gϕϕ(r, θ);
glo = FullSimplify({ {gtt(r, θ), 0, 0, 0}, {0,
grr(r, θ), 0, 0}, {0, 0, gθθ(r, θ),
0}, {0, 0, 0, gϕϕ(r, θ)}});
gup = Simplify(Inverse(glo));
dglo = Simplify(Det(glo));
crd = {t, r, θ, ϕ};
Xup = {t(τ), r(τ), θ(τ), ϕ(τ)};
Vup = {Vt, Vr, Vθ, Vϕ};(* v^μ vector *)
Pup = {Pt(τ), Pr(τ), Pθ(τ), Pϕ(τ)};
Sup = {{Stt(τ), Str(τ), Stθ(τ),
Stϕ(τ)},
{Srt(τ), Srr(τ), Srθ(τ), Srϕ(τ)},
{Sθt(τ), Sθr(τ), Sθθ(τ),
Sθϕ(τ)},
{Sϕt(τ), Sϕr(τ), Sϕθ(τ),
Sϕϕ(τ)}};
christoffel =
Simplify(Table((1/2)*
Sum((gup((i, s)))*(D(glo((s, k)), crd((j)) ) +
D(glo((s, j)), crd((k)) ) - D(glo((j, k)), crd((s)) )), {s,
1, n}), {i, 1, n}, {j, 1, n}, {k, 1, n})
);
riemann = Simplify(
Table(
D(christoffel((i, j, l)), crd((k)) ) -
D(christoffel((i, j, k)), crd((l)) ) +
Sum(christoffel((s, j, l)) christoffel((i, k, s)) -
christoffel((s, j, k)) christoffel((i, l, s)),
{s, 1, n}), {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}) );
loriemann =
Simplify(Table(
Sum(glo((i, m))*riemann((m, j, k, l)), {m, 1, n}), {i, 1, n}, {j,
1, n}, {k, 1, n}, {l, 1, n}) );
EOM1 = Table( D(Xup((a)), τ) == Vup((a)) , {a, 1, n});
EOM2 = Table(
D(Pup((a)), τ) + !(
*UnderoverscriptBox((∑), (b = 1), (n))(
*UnderoverscriptBox((∑), (c =
1), (n))christoffel((()(a, b, c)()))*
Pup((()(b)()))*Vup((()(c)())))) == -(1/2) !(
*UnderoverscriptBox((∑), (b = 1), (n))(
*UnderoverscriptBox((∑), (c = 1), (n))(
*UnderoverscriptBox((∑), (d = 1), (n))riemann((()(a,
b, c, d)()))*Vup((()(b)()))*
Sup((()(c, d)()))))),
{a, 1, n});
EOM3 = Table(
D(Sup((a, b)), τ) + !(
*UnderoverscriptBox((∑), (c = 1), (n))(
*UnderoverscriptBox((∑), (d =
1), (n))christoffel((()(a, c, d)()))*
Sup((()(c, b)()))*Vup((()(d)())))) + !(
*UnderoverscriptBox((∑), (c = 1), (n))(
*UnderoverscriptBox((∑), (d =
1), (n))christoffel((()(b, c, d)()))*
Sup((()(a, c)()))*Vup((()(d)())))) ==
Pup((a))*Vup((b)) - Pup((b))*Vup((a)),
{a, 1, n}, {b, 1, n});
Wfactor = Simplify(4*μ^2 + !(
*UnderoverscriptBox((∑), (i = 1), (4))(
*UnderoverscriptBox((∑), (j = 1), (4))(
*UnderoverscriptBox((∑), (k = 1), (4))(
*UnderoverscriptBox((∑), (l =
1), (4))((loriemann((()(i, j, k,
l)()))*((Sup((()(i, j)()))))* ((Sup((()(k,
l)())))))))))));
Wvec = Simplify(Table(2/(μ*Wfactor)*(!(
*UnderoverscriptBox((∑), (i = 1), (4))(
*UnderoverscriptBox((∑), (k = 1), (4))(
*UnderoverscriptBox((∑), (m = 1), (4))(
*UnderoverscriptBox((∑), (l = 1), (4))Sup((()(j,
i)()))*
Pup((()(k)()))*((loriemann((()(i, k, l,
m)()))))*((Sup((()(l, m)()))))))))), {j,
1, n}));
NN = 1/Sqrt(1 - !(
*UnderoverscriptBox((∑), (i = 1), (4))(
*UnderoverscriptBox((∑), (k =
1), (4))((glo((()(i, k)()))))*Wvec((()(i)()))*
Wvec((()(k)())))));
{Vt, Vr, Vθ, Vϕ} = NN (Wvec + Pup);
EOM = Flatten(
Join({EOM1, EOM2, EOM3} /.
r -> r(τ) /. θ -> θ(τ) /.
Derivative(1)(r(τ))(τ) -> Derivative(1)(r)(τ) /.
Derivative(1)(θ(τ))(τ) ->
Derivative(1)(θ)(τ)));
(**********************************************************************************************************************************)
Sθϕ0 = (
LL Cot(θ0))/r0^2; Srθ0 = -(pθ0/
r0); Srϕ0 = -(1/
r0) (-LL + pϕ0/Sin(θ0)^2); Str0 = -(1/(
r0*pt0)) (pθ0^2 + pϕ0^2/Sin(θ0)^2 - LL*pϕ0);
Stθ0 =
1/(r0*pt0)*(pr0*pθ0 + (LL * pϕ0)/r0*
Cot(θ0)); Stϕ0 = -(1/(
r0*pt0))*(LL*pr0 - (pr0*pϕ0)/
Sin(θ0)^2 + (LL*pθ0)/r0*Cot(θ0));
glo0 = Simplify(glo /. {r -> r0, θ -> θ0});
gup0 = Simplify(Inverse(glo0));
plo0 = {pt0, pr0, pθ0, pϕ0};
Sup0 = {{0, Str0, Stθ0, Stϕ0},
{-Str0, 0, Srθ0, Srϕ0},
{-Stθ0, -Srθ0, 0, Sθϕ0},
{-Stϕ0, -Srϕ0, -Sθϕ0, 0}};
F0 = Simplify(μ^2 + !(
*UnderoverscriptBox((∑), (a = 1), (4))(
*UnderoverscriptBox((∑), (b =
1), (4))((((gup0((()(a,
b)()))))*((plo0((()(a)()))))*((plo0((()(b)
())))))))));
F1 = EE + pt0 -
M/r0^2*(1/(
r0*pt0)*(pθ0^2 + pϕ0^2/Sin(θ0)^2 -
LL*pϕ0));
F2 = Simplify(SS^2 - 1/2*(!(
*UnderoverscriptBox((∑), (a = 1), (4))(
*UnderoverscriptBox((∑), (b = 1), (4))(
*UnderoverscriptBox((∑), (c = 1), (4))(
*UnderoverscriptBox((∑), (d =
1), (4))((glo0((()(a, b)()))*
glo0((()(c, d)()))*Sup0((()(a, c)()))*
Sup0((()(b, d)()))))))))));
pth = N(-(μ^2 + (gup0((1, 1)) )*(- EE)^2 +
gup0((2, 2)) * (0)^2 + (gup0((4, 4)) )* (LL)^2)/(gup0((3,
3)) ), 40);
pθGR = N(If(pth < 0, 1, Sqrt(pth)), 40);
GiveMePSpoints(SS_, r0_, θ0_, EE_, LL_, konec_, color_) := (
Clear(Str0, Stθ0, Stϕ0, Srθ0, Srϕ0,
Sθϕ0, pt0, pr0, pθ0, pϕ0);
μ = N(1, 40); pr0 = N(0, 40);
{pt0, pr0, pθ0,
pϕ0} = {pt0, pr0, pθ0, pϕ0} /.
FindRoot({F0 == 0, F1 == 0,
F2 == 0}, {{pt0, -EE}, {pθ0, pθGR}, {pϕ0,
LL}}, WorkingPrecision -> 39);
Sθϕ0 = (LL Cot(θ0))/r0^2;
Srθ0 = -(pθ0/r0);
Srϕ0 = -(1/r0) (-LL + pϕ0/Sin(θ0)^2);
Str0 = -(1/(
r0*pt0)) (pθ0^2 + pϕ0^2/Sin(θ0)^2 -
LL*pϕ0);
Stθ0 =
1/(r0*pt0)*(pr0*pθ0 + (LL * pϕ0)/r0*Cot(θ0));
Stϕ0 = -(1/(
r0*pt0))*(LL*pr0 - (pr0*pϕ0)/
Sin(θ0)^2 + (LL*pθ0)/r0*Cot(θ0));
INT1 = {t(0) == 0,
r(0) == r0, θ(0) == θ0, ϕ(0) == 0};
INT2 = {Pt(0) == gUtt(r0, θ0) pt0,
Pr(0) == gUrr(r0, θ0) pr0,
Pθ(0) == gUθθ(r0, θ0) pθ0,
Pϕ(0) == gUϕϕ(r0, θ0) pϕ0};
INT3 = {{Stt(0) == 0, Str(0) == Str0, Stθ(0) == Stθ0,
Stϕ(0) == Stϕ0},
{Srt(0) == -Str0, Srr(0) == 0, Srθ(0) == Srθ0,
Srϕ(0) == Srϕ0},
{Sθt(0) == -Stθ0, Sθr(0) == -Srθ0,
Sθθ(0) == 0,
Sθϕ(0) == Sθϕ0},
{Sϕt(0) == -Stϕ0, Sϕr(0) == -Srϕ0,
Sϕθ(0) == -Sθϕ0, Sϕϕ(0) == 0}};
INT = Flatten(Join({INT1, INT2, INT3}));
SetSystemOptions(
"NDSolveOptions" -> "DefaultSolveTimeConstraint" -> 100.`);
data =
Reap(NDSolve(
Flatten(Join({EOM, INT})), {t, r, θ, ϕ, Pt, Pr,
Pθ, Pϕ, Stt, Str, Stθ, Stϕ, Srt, Srr,
Srθ, Srϕ,
Sθt, Sθr, Sθθ, Sθϕ,
Sϕt, Sϕr, Sϕθ, Sϕϕ}, {τ,
0, konec},
StartingStepSize ->
0.1, {Method -> {"EventLocator",
"Event" -> θ(τ) - Pi/2,
"EventAction" :> Sow({r(τ), Pr(τ)}),
"Method" -> {"FixedStep",
"Method" -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10,
"ImplicitSolver" -> {"Newton",
AccuracyGoal -> MachinePrecision,
PrecisionGoal -> MachinePrecision,
"IterationSafetyFactor" -> 1}}}}}));
psdata = Take(data((-1, 1)));
ListPlot(psdata, PlotStyle -> {{PointSize(0.003), color}},
Frame -> True, Axes -> False, BaseStyle -> 15, ImageSize -> 350,
AspectRatio -> 1, PlotRange -> range)
);
θ0 =
N(Pi/2, 40); r0 =
4.2521489655172413793103448275862068965517241379310344827585`40.;
EE = N(976/1000, 40); LL = N(38/10, 40); SS = N(10^-4, 40);
konec = 7 10^5; range = {{4.25213, 4.25218}, {-15 10^-6, 15 10^-6}};
GiveMePSpoints(SS, r0, θ0, EE, LL, konec, Black)
