Expand Taylor series about nonzero points

src/taylor.jl

@@ -1,10 +1,9 @@
-# TODO: around arbitrary point x = x0 ≠ 0
 # TODO: error if x is not a "pure variable"
 # TODO: optimize for multiple orders with loop/recursion
 """
-    taylor(f, x, n; rationalize=true)
+    taylor(f, x, [x0=0,] n; rationalize=true)
 
-Calculate the `n`-th order term(s) in the Taylor series of the expression `f(x)` around `x = 0`.
+Calculate the `n`-th order term(s) in the Taylor series of the expression `f(x)` around `x = x0`.
 If `rationalize`, float coefficients are approximated as rational numbers (this can produce unexpected results for irrational numbers, for example).
 
 Examples
@@ -19,6 +18,9 @@ julia> taylor(exp(x), x, 0:3)
 
 julia> taylor(exp(x), x, 0:3; rationalize=false)
 1.0 + x + 0.5(x^2) + 0.16666666666666666(x^3)
+
+julia> taylor(√(x), x, 1, 0:3)
+1 + (1//2)*(-1 + x) - (1//8)*((-1 + x)^2) + (1//16)*((-1 + x)^3)
 ```
 """
 function taylor(f, x, n::Int; rationalize=true)
@@ -37,3 +39,15 @@ end
 function taylor(f, x, n::AbstractArray{Int}; kwargs...)
     return sum(taylor.(f, x, n; kwargs...))
 end
+function taylor(f, x, x0, n; kwargs...)
+    # 1) substitute dummy x′ = x - x0
+    name = Symbol(String(nameof(x)) * "′") # e.g. Symbol("x′")
+    x′ = only(@variables $name)
+    f = substitute(f, x => x′ + x0)
+
+    # 2) expand f around x′ = 0
+    s = taylor(f, x′, n; kwargs...)
+
+    # 3) substitute back x = x′ + x0
+    return substitute(s, x′ => x - x0)
+end

test/taylor.jl

@@ -23,3 +23,6 @@ end
 @test taylor(sinh(x), x, 0:7) - sum(1/factorial(2*n+1) * x^(2*n+1) for n in 0:3) == 0
 @test taylor(cosh(x), x, 0:7) - sum(1/factorial(2*n) * x^(2*n) for n in 0:3) == 0
 @test taylor(tanh(x), x, 0:7) - (x - x^3/3 + 2/15*x^5 - 17/315*x^7) == 0
+
+# around x ≠ 0
+@test substitute(taylor(√(x), x, 1, 0:6), x => x + 1) - taylor(√(1+x), x, 0:6) == 0