I would like to know if anyone has a clear understanding on the following parameters used to control the precision/accuracy of NDSolve.
I have been working on the solution of the following differential equation
x'[s] = Cos[theta[s]]
y'[s] = Sin[theta[s]]
theta'[s] = 2 b - Sin[theta[s]] / x[s] (for s >0, and b for s = 0)
with the initial condition x = 0, y = 0, theta = 0.
This differential equation has an exact solution x[s] = (1/b) Sin[b s], y[s] = (1/b) (1 - Cos[b s]), theta[s] = b s. Therefore, I can check the precision of the result. So I launch Mathematica and solve it with b = 1, PrecisionGoal -> 20, WorkingPrecision -> 30. But I can't seem to get something better than 9 digits of precision.
I have attached a Mathematica notebook to prove that.
I kind of understand that there is no way an algorithm can guarantee a given precision. But the documentation around PrecisionGoal looks pretty unclear to me. So here are my questions :
From the Mathematica 10 documentation:
In this case, the internal rule is not really wise. As you could ask NDSolve[..., PrecisionGoal -> 30, ... ] without having any guarantee on the precision of the solution outside the computations points which are hidden in the "interpolation" object. I find this behaviour disturbing.
According to the documentation, the default interpolating order is "Automatic". This probably means that Mathematica has some internal rule to chose the order.
I finally found the solution. NDSolve only guarantees the precision at the points where the solution is computed (which you don't really now expect if you inspect the interpolation object). A linear interpolation is made in between the points is made if no option is given. Therefore, the precision at some points could be much larger than the expected precision.
In order to have something better than a linear approximation, one should use the option InterpolationOrder->All . The precision seems to be the expected one everywhere on the interval. This is very misleading and it is strange that the option is not the default one.
Try setting an AccuracyGoal as well.