7
$\begingroup$

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.

Improve this question
$\endgroup$
14
  • 1
    $\begingroup$ does fx = First[F@*x /. s] give what you need? $\endgroup$
    – kglr
    Dec 21, 2019 at 11:47
  • $\begingroup$ Well, it should be $fx(xx(r)) = FF(r)$ which I tried to verify with N[ FF[1/2] - fx[xx[1/2]], 40] i.e. with r=1/2 but only got an InterpolatingFunction::dmval error. Probably I am doing something wrong. $\endgroup$
    – user67126
    Dec 21, 2019 at 12:23
  • 1
    $\begingroup$ @jheidk51 With x[0] == 11/10, no value of x falls within the domain of F, which is {0,1}. That is the source of the error message. $\endgroup$
    – bbgodfrey
    Dec 21, 2019 at 12:28
  • $\begingroup$ Then which is the domain of fx? The problem arises with respect to the domain of fx not F! $\endgroup$
    – user67126
    Dec 21, 2019 at 12:44
  • 1
    $\begingroup$ As observed in the comments, you need to invert xx[r]. Since xx'[0] == 0, this won't be possible without a loss of accuracy locally. It shouldn't be a big drawback (unless r == 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 for xx. Doing it is not hard but it would take some coding and time to do. I've used chebInterpolation (search this site) for similar work. $\endgroup$
    – Michael E2
    Dec 21, 2019 at 17:53

1 Answer 1

Reset to default
7
$\begingroup$ 
fx = Interpolation[Transpose[{xx["ValuesOnGrid"], FF["ValuesOnGrid"]}][[24 ;;]],
    InterpolationOrder -> 15]

provides FF[r[xx]] without loss of precision. The first 23 elements are dropped, because the first 24 are identical in xx. Here is what the function looks like.

Plot[fx[x], {x, 1.1, 1.435}, ImageSize -> Large, LabelStyle -> {15, Black, Bold}]

enter image description here

Addendum

Motivated by comments below by MichaelE2, I plotted the actual accuracy of fx for increasing values of InterpolationOrder. At 15, the computed accuracy is substantially better than the requested 30 digits

ParametricPlot[{xx[r], RealExponent[fx[xx[r]] - FF[r]]}, {r, 0, 1}, 
    WorkingPrecision -> 100, ImageSize -> Large, LabelStyle -> {15, Black, Bold}, 
    AspectRatio -> 1/GoldenRatio, AxesLabel -> {"xx", "fx"}]

enter image description here

except for xx less than about 1.100021410. Higher InterpolationOrdercan reduce but not eliminate this small domain of poor accuracy, because r'[xx] is singular at x = 1.1.

Improve this answer
$\endgroup$
13
  • $\begingroup$ Counting the points that are identical in x has to be done manually? $\endgroup$
    – user67126
    Dec 21, 2019 at 18:25
  • $\begingroup$ @jheidk51 No. The first time I ran Interpolation, it told me that 1.1``100 was duplicated, so I ran Position to find the duplicates. They were the first 24. I also could have used Union which might have been better. $\endgroup$
    – bbgodfrey
    Dec 21, 2019 at 18:29
  • $\begingroup$ Thanks a lot! An extra question: How can I use this method to get r[x]? $\endgroup$
    – user67126
    Dec 21, 2019 at 18:55
  • $\begingroup$ @jheidk51 Yes, you can. By the way, my suggestion just above to try Union is not a good idea. $\endgroup$
    – bbgodfrey
    Dec 21, 2019 at 19:10
  • $\begingroup$ @jheidk51 Specifically, see how I inverted a complicated function in a[s] to generate int in question 211194. $\endgroup$
    – bbgodfrey
    Dec 21, 2019 at 21:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.