In the following working example
s = NDSolve[{F'[r] == Sin[200*r]*x[r], x'[r] == F[r]*r^2, F[0] == 1,
x[0] == 11/10}, {F, x}, {r, 0, 1},
Method -> "ExplicitRungeKutta",
WorkingPrecision -> 100, AccuracyGoal -> 31, PrecisionGoal -> 31,
InterpolationOrder -> All, MaxSteps -> 10^6];
FF = First[F /. s];
xx = First[x /. s];
NDSolve
returns two functions F[r]
and x[r]
with an accuracy of about 30 decimal digits
.
x[r]
is a monotonically increasing function of r
so in principle one can combine F[r]
and x[r]
to get F[x]
.
I know that x[r]
is something like {{r1,x1},{r2,x2},{r3,x3},...}
and f[r]
is something like {{r1,f1},{r2,f2},{r3,f3},...}
. So I have to get the second column from each function and construct something like {{x1,f1},{x2,f2},{x3,f3},...}
.
Then if fx[x]
is the new function it should be fx[x[ro]] = FF[ro]
for any ro
in {0,1}
.
However I have not yet found out how.
My problem is that I have already some such interpolating F[r]
's and x[r]
's stored in .txt
form with and accuracy of 30 decimal digits
and I need F[x]
with the same accuracy.
If there is no simple solution to this then I will have to start from the beggining.
Alternative: An alternative would be to find the inverse of x[r]
i.e. r[x]
and then define fx[x_]:=F[r[x]]
when preserving accuracy at the same time. I have not yet managed this.
fx = First[F@*x /. s]
give what you need? $\endgroup$N[ FF[1/2] - fx[xx[1/2]], 40]
i.e. withr=1/2
but only got anInterpolatingFunction::dmval
error. Probably I am doing something wrong. $\endgroup$x[0] == 11/10
, no value ofx
falls within the domain ofF
, which is{0,1}
. That is the source of the error message. $\endgroup$fx
? The problem arises with respect to the domain offx
notF
! $\endgroup$xx[r]
. Sincexx'[0] == 0
, this won't be possible without a loss of accuracy locally. It shouldn't be a big drawback (unlessr == 0+
is especially important to you), but it is a technical difficulty. The bigger problem is that the highly accurate interpolants are stored as a piecewise sequence of degree-9 polynomials in Chebyshev series form. There's no easy way to invert this forxx
. Doing it is not hard but it would take some coding and time to do. I've usedchebInterpolation
(search this site) for similar work. $\endgroup$