From 5889bd51fc7a6228d76cdca92fc289f46d1c7d39 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Mon, 7 Oct 2024 20:52:42 +0200 Subject: [PATCH] Add basic Taylor series generation --- src/Symbolics.jl | 3 +++ src/taylor.jl | 31 +++++++++++++++++++++++++++++++ test/taylor.jl | 25 +++++++++++++++++++++++++ 3 files changed, 59 insertions(+) create mode 100644 src/taylor.jl create mode 100644 test/taylor.jl diff --git a/src/Symbolics.jl b/src/Symbolics.jl index 2b6a45912..561311878 100644 --- a/src/Symbolics.jl +++ b/src/Symbolics.jl @@ -140,6 +140,9 @@ export symbolic_linear_solve, solve_for include("groebner_basis.jl") export groebner_basis, is_groebner_basis +include("taylor.jl") +export taylor + import Libdl include("build_function.jl") export build_function diff --git a/src/taylor.jl b/src/taylor.jl new file mode 100644 index 000000000..c0fdf16f0 --- /dev/null +++ b/src/taylor.jl @@ -0,0 +1,31 @@ +# TODO: around arbitrary point x = x0 ≠ 0 +# TODO: error if x is not a "pure variable" +# TODO: optimize for multiple orders with loop/recursion +# TODO: get rational coefficients, not floats +""" + taylor(f, x, n) + +Calculate the `n`-th order term(s) in the Taylor series of the expression `f(x)` around `x = 0`. + +Examples +======== +```julia +julia> @variables x +1-element Vector{Num}: + x + +julia> taylor(exp(x), x, 0:3) +1.0 + x + 0.5(x^2) + 0.16666666666666666(x^3) +``` +""" +function taylor(f, x, n::Int) + D = Differential(x) + n! = factorial(n) + c = (D^n)(f) / n! + c = expand_derivatives(c) + c = substitute(c, x => 0) + return c * x^n +end +function taylor(f, x, n::AbstractArray{Int}) + return sum(taylor.(f, x, n)) +end diff --git a/test/taylor.jl b/test/taylor.jl new file mode 100644 index 000000000..4ce3edbe0 --- /dev/null +++ b/test/taylor.jl @@ -0,0 +1,25 @@ +@testset "Taylor series" begin + # https://en.wikipedia.org/wiki/Taylor_series#List_of_Maclaurin_series_of_some_common_functions + @variables x + @test taylor(exp(x), x, 0:9) - sum(1 / factorial(n) * x^n for n in 0:9) == 0 + @test taylor(log(1-x), x, 0:9) - sum(-1/n * x^n for n in 1:9) == 0 + @test taylor(log(1+x), x, 0:9) - sum((-1)^(n+1) * 1/n * x^n for n in 1:9) == 0 + + @test taylor(1/(1-x), x, 0:9) - sum(x^n for n in 0:9) == 0 + @test taylor(1/(1-x)^2, x, 0:8) - sum(n * x^(n-1) for n in 1:9) == 0 + @test taylor(1/(1-x)^3, x, 0:7) - sum((n-1)*n/2 * x^(n-2) for n in 2:9) == 0 + for α in (-1//2, 0, 1//2, 1, 2, 3) + @test taylor((1+x)^α, x, 0:7) - sum(binomial(α, n) * x^n for n in 0:7) == 0 + end + + @test taylor(sin(x), x, 0:7) - sum((-1)^n/factorial(2*n+1) * x^(2*n+1) for n in 0:3) == 0 + @test taylor(cos(x), x, 0:7) - sum((-1)^n/factorial(2*n) * x^(2*n) for n in 0:3) == 0 + @test taylor(tan(x), x, 0:7) - taylor(taylor(sin(x), x, 0:7) / taylor(cos(x), x, 0:7), x, 0:7) == 0 + @test taylor(asin(x), x, 0:7) - sum(factorial(2*n)/(4^n*factorial(n)^2*(2*n+1)) * x^(2*n+1) for n in 0:3) == 0 + @test taylor(acos(x), x, 0:7) - taylor(π/2 - asin(x), x, 0:7) == 0 + @test taylor(atan(x), x, 0:7) - taylor(asin(x/√(1+x^2)), x, 0:7) == 0 + + @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 +end