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_, (Theta)_) := -AA(r); grr(r_, (Theta)_) := 1/AA(r);
g(Theta)(Theta)(r_, (Theta)_) := r^2;
g(Phi)(Phi)(r_, (Theta)_) := r^2 Sin((Theta))^2;
gUtt(r_, (Theta)_) := 1/gtt(r, (Theta));
gUrr(r_, (Theta)_) := 1/grr(r, (Theta));
gU(Theta)(Theta)(r_, (Theta)_) := 1/g(Theta)(Theta)(r, (Theta));
gU(Phi)(Phi)(r_, (Theta)_) := 1/g(Phi)(Phi)(r, (Theta));
glo = FullSimplify({ {gtt(r, (Theta)), 0, 0, 0}, {0,
grr(r, (Theta)), 0, 0}, {0, 0, g(Theta)(Theta)(r, (Theta)),
0}, {0, 0, 0, g(Phi)(Phi)(r, (Theta))}});
gup = Simplify(Inverse(glo));
dglo = Simplify(Det(glo));
crd = {t, r, (Theta), (Phi)};
Xup = {t((Tau)), r((Tau)), (Theta)((Tau)), (Phi)((Tau))};
Vup = {Vt, Vr, V(Theta), V(Phi)};(* v^(Mu) vector *)
Pup = {Pt((Tau)), Pr((Tau)), P(Theta)((Tau)), P(Phi)((Tau))};
Sup = {{Stt((Tau)), Str((Tau)), St(Theta)((Tau)),
St(Phi)((Tau))},
{Srt((Tau)), Srr((Tau)), Sr(Theta)((Tau)), Sr(Phi)((Tau))},
{S(Theta)t((Tau)), S(Theta)r((Tau)), S(Theta)(Theta)((Tau)),
S(Theta)(Phi)((Tau))},
{S(Phi)t((Tau)), S(Phi)r((Tau)), S(Phi)(Theta)((Tau)),
S(Phi)(Phi)((Tau))}};
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)), (Tau)) == Vup((a)) , {a, 1, n});
EOM2 = Table(
D(Pup((a)), (Tau)) + !(
*UnderoverscriptBox(((Sum)), (b = 1), (n))(
*UnderoverscriptBox(((Sum)), (c =
1), (n))christoffel((()(a, b, c)()))*
Pup((()(b)()))*Vup((()(c)())))) == -(1/2) !(
*UnderoverscriptBox(((Sum)), (b = 1), (n))(
*UnderoverscriptBox(((Sum)), (c = 1), (n))(
*UnderoverscriptBox(((Sum)), (d = 1), (n))riemann((()(a,
b, c, d)()))*Vup((()(b)()))*
Sup((()(c, d)()))))),
{a, 1, n});
EOM3 = Table(
D(Sup((a, b)), (Tau)) + !(
*UnderoverscriptBox(((Sum)), (c = 1), (n))(
*UnderoverscriptBox(((Sum)), (d =
1), (n))christoffel((()(a, c, d)()))*
Sup((()(c, b)()))*Vup((()(d)())))) + !(
*UnderoverscriptBox(((Sum)), (c = 1), (n))(
*UnderoverscriptBox(((Sum)), (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*(Mu)^2 + !(
*UnderoverscriptBox(((Sum)), (i = 1), (4))(
*UnderoverscriptBox(((Sum)), (j = 1), (4))(
*UnderoverscriptBox(((Sum)), (k = 1), (4))(
*UnderoverscriptBox(((Sum)), (l =
1), (4))((loriemann((()(i, j, k,
l)()))*((Sup((()(i, j)()))))* ((Sup((()(k,
l)())))))))))));
Wvec = Simplify(Table(2/((Mu)*Wfactor)*(!(
*UnderoverscriptBox(((Sum)), (i = 1), (4))(
*UnderoverscriptBox(((Sum)), (k = 1), (4))(
*UnderoverscriptBox(((Sum)), (m = 1), (4))(
*UnderoverscriptBox(((Sum)), (l = 1), (4))Sup((()(j,
i)()))*
Pup((()(k)()))*((loriemann((()(i, k, l,
m)()))))*((Sup((()(l, m)()))))))))), {j,
1, n}));
NN = 1/Sqrt(1 - !(
*UnderoverscriptBox(((Sum)), (i = 1), (4))(
*UnderoverscriptBox(((Sum)), (k =
1), (4))((glo((()(i, k)()))))*Wvec((()(i)()))*
Wvec((()(k)())))));
{Vt, Vr, V(Theta), V(Phi)} = NN (Wvec + Pup);
EOM = Flatten(
Join({EOM1, EOM2, EOM3} /.
r -> r((Tau)) /. (Theta) -> (Theta)((Tau)) /.
Derivative(1)(r((Tau)))((Tau)) -> Derivative(1)(r)((Tau)) /.
Derivative(1)((Theta)((Tau)))((Tau)) ->
Derivative(1)((Theta))((Tau))));
(**********************************************************************************************************************************)
S(Theta)(Phi)0 = (
LL Cot((Theta)0))/r0^2; Sr(Theta)0 = -(p(Theta)0/
r0); Sr(Phi)0 = -(1/
r0) (-LL + p(Phi)0/Sin((Theta)0)^2); Str0 = -(1/(
r0*pt0)) (p(Theta)0^2 + p(Phi)0^2/Sin((Theta)0)^2 - LL*p(Phi)0);
St(Theta)0 =
1/(r0*pt0)*(pr0*p(Theta)0 + (LL * p(Phi)0)/r0*
Cot((Theta)0)); St(Phi)0 = -(1/(
r0*pt0))*(LL*pr0 - (pr0*p(Phi)0)/
Sin((Theta)0)^2 + (LL*p(Theta)0)/r0*Cot((Theta)0));
glo0 = Simplify(glo /. {r -> r0, (Theta) -> (Theta)0});
gup0 = Simplify(Inverse(glo0));
plo0 = {pt0, pr0, p(Theta)0, p(Phi)0};
Sup0 = {{0, Str0, St(Theta)0, St(Phi)0},
{-Str0, 0, Sr(Theta)0, Sr(Phi)0},
{-St(Theta)0, -Sr(Theta)0, 0, S(Theta)(Phi)0},
{-St(Phi)0, -Sr(Phi)0, -S(Theta)(Phi)0, 0}};
F0 = Simplify((Mu)^2 + !(
*UnderoverscriptBox(((Sum)), (a = 1), (4))(
*UnderoverscriptBox(((Sum)), (b =
1), (4))((((gup0((()(a,
b)()))))*((plo0((()(a)()))))*((plo0((()(b)
())))))))));
F1 = EE + pt0 -
M/r0^2*(1/(
r0*pt0)*(p(Theta)0^2 + p(Phi)0^2/Sin((Theta)0)^2 -
LL*p(Phi)0));
F2 = Simplify(SS^2 - 1/2*(!(
*UnderoverscriptBox(((Sum)), (a = 1), (4))(
*UnderoverscriptBox(((Sum)), (b = 1), (4))(
*UnderoverscriptBox(((Sum)), (c = 1), (4))(
*UnderoverscriptBox(((Sum)), (d =
1), (4))((glo0((()(a, b)()))*
glo0((()(c, d)()))*Sup0((()(a, c)()))*
Sup0((()(b, d)()))))))))));
pth = N(-((Mu)^2 + (gup0((1, 1)) )*(- EE)^2 +
gup0((2, 2)) * (0)^2 + (gup0((4, 4)) )* (LL)^2)/(gup0((3,
3)) ), 40);
p(Theta)GR = N(If(pth < 0, 1, Sqrt(pth)), 40);
GiveMePSpoints(SS_, r0_, (Theta)0_, EE_, LL_, konec_, color_) := (
Clear(Str0, St(Theta)0, St(Phi)0, Sr(Theta)0, Sr(Phi)0,
S(Theta)(Phi)0, pt0, pr0, p(Theta)0, p(Phi)0);
(Mu) = N(1, 40); pr0 = N(0, 40);
{pt0, pr0, p(Theta)0,
p(Phi)0} = {pt0, pr0, p(Theta)0, p(Phi)0} /.
FindRoot({F0 == 0, F1 == 0,
F2 == 0}, {{pt0, -EE}, {p(Theta)0, p(Theta)GR}, {p(Phi)0,
LL}}, WorkingPrecision -> 39);
S(Theta)(Phi)0 = (LL Cot((Theta)0))/r0^2;
Sr(Theta)0 = -(p(Theta)0/r0);
Sr(Phi)0 = -(1/r0) (-LL + p(Phi)0/Sin((Theta)0)^2);
Str0 = -(1/(
r0*pt0)) (p(Theta)0^2 + p(Phi)0^2/Sin((Theta)0)^2 -
LL*p(Phi)0);
St(Theta)0 =
1/(r0*pt0)*(pr0*p(Theta)0 + (LL * p(Phi)0)/r0*Cot((Theta)0));
St(Phi)0 = -(1/(
r0*pt0))*(LL*pr0 - (pr0*p(Phi)0)/
Sin((Theta)0)^2 + (LL*p(Theta)0)/r0*Cot((Theta)0));
INT1 = {t(0) == 0,
r(0) == r0, (Theta)(0) == (Theta)0, (Phi)(0) == 0};
INT2 = {Pt(0) == gUtt(r0, (Theta)0) pt0,
Pr(0) == gUrr(r0, (Theta)0) pr0,
P(Theta)(0) == gU(Theta)(Theta)(r0, (Theta)0) p(Theta)0,
P(Phi)(0) == gU(Phi)(Phi)(r0, (Theta)0) p(Phi)0};
INT3 = {{Stt(0) == 0, Str(0) == Str0, St(Theta)(0) == St(Theta)0,
St(Phi)(0) == St(Phi)0},
{Srt(0) == -Str0, Srr(0) == 0, Sr(Theta)(0) == Sr(Theta)0,
Sr(Phi)(0) == Sr(Phi)0},
{S(Theta)t(0) == -St(Theta)0, S(Theta)r(0) == -Sr(Theta)0,
S(Theta)(Theta)(0) == 0,
S(Theta)(Phi)(0) == S(Theta)(Phi)0},
{S(Phi)t(0) == -St(Phi)0, S(Phi)r(0) == -Sr(Phi)0,
S(Phi)(Theta)(0) == -S(Theta)(Phi)0, S(Phi)(Phi)(0) == 0}};
INT = Flatten(Join({INT1, INT2, INT3}));
SetSystemOptions(
"NDSolveOptions" -> "DefaultSolveTimeConstraint" -> 100.`);
data =
Reap(NDSolve(
Flatten(Join({EOM, INT})), {t, r, (Theta), (Phi), Pt, Pr,
P(Theta), P(Phi), Stt, Str, St(Theta), St(Phi), Srt, Srr,
Sr(Theta), Sr(Phi),
S(Theta)t, S(Theta)r, S(Theta)(Theta), S(Theta)(Phi),
S(Phi)t, S(Phi)r, S(Phi)(Theta), S(Phi)(Phi)}, {(Tau),
0, konec},
StartingStepSize ->
0.1, {Method -> {"EventLocator",
"Event" -> (Theta)((Tau)) - Pi/2,
"EventAction" :> Sow({r((Tau)), Pr((Tau))}),
"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)
);
(Theta)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, (Theta)0, EE, LL, konec, Black)
```