Airy function in entire complex plane #51
Conversation
Codecov ReportBase: 94.08% // Head: 94.17% // Increases project coverage by
Additional details and impacted files@@ Coverage Diff @@
## master #51 +/- ##
==========================================
+ Coverage 94.08% 94.17% +0.08%
==========================================
Files 18 18
Lines 1539 1820 +281
==========================================
+ Hits 1448 1714 +266
- Misses 91 106 +15
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. ☔ View full report at Codecov. |
|
Ok I've updated the docs and added the appropriate limits for infinite arguments therefore no error is return they will simply return the branches defined by the exponential or 1 / z. There's one thing I'm still thinking about is how do we want to appropriately promote types in some of the series. For example the besseli series can always be computed without loss of precision but only for real arguments. For complex argument we will need to promote ComplexF32 to ComplexF64 and same happens here for the airy functions ( |
|
Ok I think this is ready to go. I've added tests for Float32 arguments. This function oscillates around the zero axis for negative arguments so it is just easy to promote these to Float64. I'm open to other ideas the best way to do this but went with a structure like. function airyai_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T
S = eltype(x)
...
...
for i in 0:MaxIter
...
abs(t) < tol * abs(ai1) && break
...
end
return ...
end
airyai_power_series(x::Float32) = Float32(airyai_power_series(Float64(x), tol=eps(Float32)))
airyai_power_series(x::ComplexF32) = ComplexF32(airyai_power_series(ComplexF64(x), tol=eps(Float32)))Here is a plot of the relative errors in the entire complex plane compared to the special functions implementation which is only accurate to like 1e-12-1e-13 in some domains anyway. |
|
@oscardssmith Ok with merging? |
This PR gives dedicated routines for the airy functions instead of relying on the relations to Bessel functions. The first thing this allows is significantly more accuracy than we had before in whole domain (for airybi these are roughly 4 more digits accurate). I've updated the tests to reflect that.
The second thing is they are much faster now.
So having dedicated routines allows for much better accuracy and speed. Of course, I have also added support for the entire complex plane. The motivation for this is some of the uniform asymptotic expansions for bessel functions require evalulation of the airy functions (https://dlmf.nist.gov/10.20) so this may come in handy as we look at those for complex plane.
Of course there are still some domains I'm not too happy with in the complex plane. I will add notes in the file on this. I chatted with @cgeoga and a potential use here is to use the hypergeometric functions or the exponentially improved functions to try and bridge domains. We can piece that together better in the future.
I still have a couple of things to finish up on this. Mainly tidying up and adding docs to this but also getting rid of some of those constant gamma calls. There is a lot of code overlap that could be simplified especially in the asymptotic expansions. However, what we want to overlap with will depend if we want to compute both airyai and airyaiprime or airyai and airybi. Just depends. It might be something we revisit in the future. Looking at the asymptotic expansions we will want airyai and airyprime so it might best to combine those....
Need to also go through and think more carefully about type promotions. I don't think we need to actually promote to float64 in a lot of the algorithms this is just left over from things I was trying.