Astroforum  don't use UiO password
AST 5220 / AST9420 => Milestone II => Topic started by: winther on 15. March 2021, 17:56:20

We donâ€™t have an initial condition for the ODE for the optical depth \(\tau\) if we solve it from the early Universe till today. However we can just pick any initial value, solve the ODE and then normalize it by subtracting off the value at \(x=0\) to get a function that has \( \tau(0) = 0 \). This works fine, but there is one issue to be aware of. The optical depth we subtract off doing this normalization can be huge, lets say, \(10^{10} \). The resulting function \(\tau = \tau_{\rm ODE}  10^{10}\) will then have \(\tau(0) \equiv 0\) as we want, but it will only have an (absolute) numerical accuracy of just \(\sim 10^{6}\) in this case (double precision means a relative precision of \(10^{16}\)). Thus close to \(x=0\) where \(\tau\) is small we might get some numerical noise popping up. This might not be seen in the function itself, but it will be seen in the derivatives which is more sensitive to this.
Solutions:
Fill an array with \(d\tau(x)/dx\) just as you do in the ODE solver using the formula for this. Make a spline of this and use this spline to get \(d\tau(x)/dx\) and its derivative to get \(d^2\tau(x)/dx^2\). This usually works very well. Another option is to just solve the ODE from today and back in time (I just updated the ODEsolver on GitHub to allow this) and then its no issue.

So I downloaded the ODEsolver to be able to do the backwards integration, but now the code won't compile, I get the following error:
"more than one instance of overloaded function "std::abs" matches the argument list"
the error appears to be at line 78 in the code.

Thanks, that is probably just a missing header file (some compilers add standard math functions automatically, but some require you to specify the header). Can you try to add
#include <cmath>
to the list on the top in ODESolver.h and see if it works?
And if that doesn't work try to change that call on line 78 in ODESolver.cpp to
std::fabs(hstart)
i.e. from abs to fabs.

Yes that fixed it :)
#include <cmath>
was all that I needed

Great, thanks. Good to fix some bugs, I'll upload the fix. Almost all of you have found small bugs in the template this year!