Implement special functions: erf, erfc, etc.

[ChristopherMayes]

How can we implement special functions, such as erf?

using TaylorSeries, SpecialFunctions 
t = Taylor1(Float64, 5)
erf(t)

returns:

MethodError: no method matching erf(::Taylor1{Float64})

[lbenet]

Thanks for reaching out!

The basic approach to implement the Taylor series of any function begins on the differential equation that the function satisfies. The equation is solved by assuming that the solution has a polynomial (Taylor) expansion, and the coefficients of that polynomial are obtained recursively. The initial condition of the differential equation is related to the value of the function you are interested at the point around which you want to Taylor expand. See the docs here for more details. My recommendation is that, as an exercise, you try to obtain the coefficients of exp(f(t)), where f(t) is a known polynomial, from the differential equation it has to satisfy. (You only need to use the product of polynomials at some point, and differentiate polynomials.)

If you do not want to implement the most general form, implement the Taylor series around 0 for erf and then use update! in small-enough steps, until you get to the point around which you want the Taylor expansion. But be aware that floating point error may accumulate.

I hope this helps.

[ChristopherMayes]

The derivative of erf is defined in terms of elementary functions in closed form:

https://en.wikipedia.org/wiki/Error_function

It looks like Julia ultimately calls a C function:
https://specialfunctions.juliamath.org/stable/functions_list/

There are many such functions. It would be elegant if just this knowledge could be coded - is that possible?

[lbenet]

I do think it is possible to code it; the next two weeks I'm too busy. If you want to give it a try, I can help.

[ChristopherMayes]

@lbenet thank you for the quick replies! I don't have the skills to do this at the moment. But, to be clear, what I mean is is to write some sort of code generation function that takes an existing scalar function (e.g. SpecialFunctions.erf) and its derivative function (e.g. derivative_erf(x) = ... or better if there were an existing SpecialFunctions.derivative_erf) that's defined in terms of functions that TaylorSeries.jl already understands. Then many special functions could be defined with one-liners.