From 0b1bff915ef40e57afad539b8718525309c7f65b Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Mon, 7 Oct 2024 20:52:42 +0200 Subject: [PATCH 01/17] Add basic Taylor series generation --- src/Symbolics.jl | 3 +++ src/taylor.jl | 31 +++++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/taylor.jl | 25 +++++++++++++++++++++++++ 4 files changed, 60 insertions(+) create mode 100644 src/taylor.jl create mode 100644 test/taylor.jl diff --git a/src/Symbolics.jl b/src/Symbolics.jl index cd84fd16a..4e20e91a0 100644 --- a/src/Symbolics.jl +++ b/src/Symbolics.jl @@ -138,6 +138,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/runtests.jl b/test/runtests.jl index 8db2d5c9d..b42a8cae4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,6 +61,7 @@ if GROUP == "All" || GROUP == "Core" @safetestset "Utility Function Test" begin include("utils.jl") end @safetestset "RootFinding solver" begin include("solver.jl") end @safetestset "Function inverses test" begin include("inverse.jl") end + @safetestset "Taylor Series Test" begin include("taylor.jl") end end end diff --git a/test/taylor.jl b/test/taylor.jl new file mode 100644 index 000000000..3c0fd99f9 --- /dev/null +++ b/test/taylor.jl @@ -0,0 +1,25 @@ +using Symbolics + +# 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 From 45fdba019fdd685e9cd6fe6bacc0997a8964cc1d Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 8 Oct 2024 18:56:10 +0200 Subject: [PATCH 02/17] Rationalize float coefficients --- src/taylor.jl | 18 +++++++++++++----- test/taylor.jl | 12 ++++++------ 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/taylor.jl b/src/taylor.jl index c0fdf16f0..4b31c99f5 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -1,11 +1,11 @@ # 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) + taylor(f, x, n; rationalize=true) Calculate the `n`-th order term(s) in the Taylor series of the expression `f(x)` around `x = 0`. +If `rationalize`, float coefficients are approximated as rational numbers (this can produce unexpected results for irrational numbers, for example). Examples ======== @@ -15,17 +15,25 @@ julia> @variables x x julia> taylor(exp(x), x, 0:3) +1 + x + (1//2)*(x^2) + (1//6)*(x^3) + +julia> taylor(exp(x), x, 0:3; rationalize=false) 1.0 + x + 0.5(x^2) + 0.16666666666666666(x^3) ``` """ -function taylor(f, x, n::Int) +function taylor(f, x, n::Int; rationalize=true) D = Differential(x) n! = factorial(n) c = (D^n)(f) / n! c = expand_derivatives(c) c = substitute(c, x => 0) + if rationalize && unwrap(c) isa Number + # TODO: make rational coefficients "organically" and not using rationalize (see https://github.com/JuliaSymbolics/Symbolics.jl/issues/1299) + c = unwrap(c) + c = c % 1.0 == 0.0 ? Int(c) : Base.rationalize(c) # avoid nameclash with keyword argument + end return c * x^n end -function taylor(f, x, n::AbstractArray{Int}) - return sum(taylor.(f, x, n)) +function taylor(f, x, n::AbstractArray{Int}; kwargs...) + return sum(taylor.(f, x, n; kwargs...)) end diff --git a/test/taylor.jl b/test/taylor.jl index 3c0fd99f9..6fcb5cb99 100644 --- a/test/taylor.jl +++ b/test/taylor.jl @@ -2,22 +2,22 @@ using Symbolics # 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(exp(x), x, 0:9) - sum(x^n//factorial(n) for n in 0:9) == 0 +@test taylor(log(1-x), x, 0:9) - sum(-x^n/n for n in 1:9) == 0 +@test taylor(log(1+x), x, 0:9) - sum((-1)^(n+1)*x^n/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 +@test taylor(1/(1-x)^3, x, 0:7) - sum((n-1)*n*x^(n-2)/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 + @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(acos(x), x, 0:7) - (π/2 - taylor(asin(x), x, 0:7)) == 0 # TODO: make π/2 a proper fraction (like Num(π)/2) @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 From 44024511f48393b0b4300572d53009fa60bf80c1 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 8 Oct 2024 19:39:47 +0200 Subject: [PATCH 03/17] Expand Taylor series about nonzero points --- src/taylor.jl | 20 +++++++++++++++++--- test/taylor.jl | 3 +++ 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/src/taylor.jl b/src/taylor.jl index 4b31c99f5..0fc5e7a8b 100644 --- a/src/taylor.jl +++ b/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(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 diff --git a/test/taylor.jl b/test/taylor.jl index 6fcb5cb99..de927e158 100644 --- a/test/taylor.jl +++ b/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 From a8646f86cb9d878ca64dfc78ba260c019116c688 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 8 Oct 2024 20:23:53 +0200 Subject: [PATCH 04/17] Compute Taylor series coefficients and full Taylor series separately --- src/taylor.jl | 39 ++++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/src/taylor.jl b/src/taylor.jl index 0fc5e7a8b..444702c00 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -1,5 +1,23 @@ -# TODO: error if x is not a "pure variable" -# TODO: optimize for multiple orders with loop/recursion +""" + taylor_coeff(f, x, n; rationalize=true) + +Calculate the `n`-th order coefficient(s) in the Taylor series of the expression `f(x)` around `x = 0`. +""" +function taylor_coeff(f, x, n; rationalize=true) + # TODO: error if x is not a "pure variable" + D = Differential(x) + n! = factorial(n) + c = (D^n)(f) / n! # TODO: optimize the implementation for multiple n with a loop that avoids re-differentiating the same expressions + c = expand_derivatives(c) + c = substitute(c, x => 0) + if rationalize && unwrap(c) isa Number + # TODO: make rational coefficients "organically" and not using rationalize (see https://github.com/JuliaSymbolics/Symbolics.jl/issues/1299) + c = unwrap(c) + c = isinteger(c) ? Int(c) : Base.rationalize(c) # convert non-integer floats to rational numbers; avoid name clash between rationalize and Base.rationalize() + end + return c +end + """ taylor(f, x, [x0=0,] n; rationalize=true) @@ -23,21 +41,8 @@ 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) - D = Differential(x) - n! = factorial(n) - c = (D^n)(f) / n! - c = expand_derivatives(c) - c = substitute(c, x => 0) - if rationalize && unwrap(c) isa Number - # TODO: make rational coefficients "organically" and not using rationalize (see https://github.com/JuliaSymbolics/Symbolics.jl/issues/1299) - c = unwrap(c) - c = c % 1.0 == 0.0 ? Int(c) : Base.rationalize(c) # avoid nameclash with keyword argument - end - return c * x^n -end -function taylor(f, x, n::AbstractArray{Int}; kwargs...) - return sum(taylor.(f, x, n; kwargs...)) +function taylor(f, x, ns; kwargs...) + return sum(taylor_coeff(f, x, n; kwargs...) * x^n for n in ns) end function taylor(f, x, x0, n; kwargs...) # 1) substitute dummy x′ = x - x0 From 496db61afcb5c9d9e4718754f73cee9d0c794939 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 8 Oct 2024 23:20:37 +0200 Subject: [PATCH 05/17] Taylor expand equations --- src/taylor.jl | 20 ++++++++++++++++++-- test/taylor.jl | 7 +++++++ 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/taylor.jl b/src/taylor.jl index 444702c00..e45615f33 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -1,9 +1,21 @@ """ taylor_coeff(f, x, n; rationalize=true) -Calculate the `n`-th order coefficient(s) in the Taylor series of the expression `f(x)` around `x = 0`. +""" + taylor_coeff(f, x[, n]; rationalize=true) + +Calculate the `n`-th order coefficient(s) in the Taylor series of `f` around `x = 0`. """ function taylor_coeff(f, x, n; rationalize=true) + if n isa AbstractArray + # return array of expressions/equations for each order + return taylor_coeff.(Ref(f), Ref(x), n; rationalize) + elseif f isa Equation + # return new equation with coefficients of each side + return taylor_coeff(f.lhs, x, n; rationalize) ~ taylor_coeff(f.rhs, x, n; rationalize) + end + end + # TODO: error if x is not a "pure variable" D = Differential(x) n! = factorial(n) @@ -21,7 +33,7 @@ end """ 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 = x0`. +Calculate the `n`-th order term(s) in the Taylor series of `f` around `x = x0`. If `rationalize`, float coefficients are approximated as rational numbers (this can produce unexpected results for irrational numbers, for example). Examples @@ -42,6 +54,10 @@ julia> taylor(√(x), x, 1, 0:3) ``` """ function taylor(f, x, ns; kwargs...) + if f isa Equation + return taylor(f.lhs, x, ns; kwargs...) ~ taylor(f.rhs, x, ns; kwargs...) + end + return sum(taylor_coeff(f, x, n; kwargs...) * x^n for n in ns) end function taylor(f, x, x0, n; kwargs...) diff --git a/test/taylor.jl b/test/taylor.jl index de927e158..de6faf811 100644 --- a/test/taylor.jl +++ b/test/taylor.jl @@ -26,3 +26,10 @@ end # around x ≠ 0 @test substitute(taylor(√(x), x, 1, 0:6), x => x + 1) - taylor(√(1+x), x, 0:6) == 0 + +# equations +eq = sin(2*x) ~ 2*sin(x)*cos(x) +eq = taylor(eq, x, 0:7) +eqs = taylor_coeff(eq, x) # should automatically expand to 7th order +@test length(eqs) == 7+1 && all(isequal(eq.lhs, eq.rhs) for eq in eqs) + From 2cadd89d55030d6a61a39c64f1f6cbcbfc8fcfe7 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Wed, 9 Oct 2024 22:17:44 +0200 Subject: [PATCH 06/17] If not set, determine Taylor expansion order automatically from variable degree --- src/Symbolics.jl | 2 +- src/taylor.jl | 15 ++++++++++++--- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/Symbolics.jl b/src/Symbolics.jl index 4e20e91a0..ed95b8338 100644 --- a/src/Symbolics.jl +++ b/src/Symbolics.jl @@ -139,7 +139,7 @@ include("groebner_basis.jl") export groebner_basis, is_groebner_basis include("taylor.jl") -export taylor +export taylor, taylor_coeff import Libdl include("build_function.jl") diff --git a/src/taylor.jl b/src/taylor.jl index e45615f33..2a6e2f08e 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -6,14 +6,23 @@ Calculate the `n`-th order coefficient(s) in the Taylor series of `f` around `x = 0`. """ -function taylor_coeff(f, x, n; rationalize=true) +function taylor_coeff(f, x, n = missing; rationalize=true) if n isa AbstractArray # return array of expressions/equations for each order return taylor_coeff.(Ref(f), Ref(x), n; rationalize) elseif f isa Equation - # return new equation with coefficients of each side - return taylor_coeff(f.lhs, x, n; rationalize) ~ taylor_coeff(f.rhs, x, n; rationalize) + if ismissing(n) + # assume user wants maximum order in the equation + n = 0:max(degree(f.lhs, x), degree(f.rhs, x)) + return taylor_coeff(f, x, n; rationalize) + else + # return new equation with coefficients of each side + return taylor_coeff(f.lhs, x, n; rationalize) ~ taylor_coeff(f.rhs, x, n; rationalize) end + elseif ismissing(n) + # assume user wants maximum order in the expression + n = 0:degree(f, x) + return taylor_coeff(f, x, n; rationalize) end # TODO: error if x is not a "pure variable" From 15d3310c9ebc2c6a0577e5a4f6c229b01224b5f0 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Wed, 9 Oct 2024 22:30:53 +0200 Subject: [PATCH 07/17] Add convenience method for expanding a variable in a series --- src/Symbolics.jl | 2 +- src/taylor.jl | 12 +++++++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/src/Symbolics.jl b/src/Symbolics.jl index ed95b8338..146841893 100644 --- a/src/Symbolics.jl +++ b/src/Symbolics.jl @@ -139,7 +139,7 @@ include("groebner_basis.jl") export groebner_basis, is_groebner_basis include("taylor.jl") -export taylor, taylor_coeff +export series, taylor, taylor_coeff import Libdl include("build_function.jl") diff --git a/src/taylor.jl b/src/taylor.jl index 2a6e2f08e..34150cbf3 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -1,5 +1,15 @@ """ - taylor_coeff(f, x, n; rationalize=true) + series(y, x, [x0=0,] ns; name = nameof(y)) + +Expand the variable `y` in a power series in the variable `x` around `x0` to orders `ns`. +""" +function series(y, x, x0, ns; name = nameof(y)) + c, = @variables $name[ns] + return sum(c[n] * (x - x0)^n for n in ns) +end +function series(y, x, ns; kwargs...) + return series(y, x, 0, ns; kwargs...) +end """ taylor_coeff(f, x[, n]; rationalize=true) From 99db4e87185717c5c616d98625f4884a71c9338e Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Wed, 9 Oct 2024 22:31:27 +0200 Subject: [PATCH 08/17] Add quintic equation test --- test/taylor.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/taylor.jl b/test/taylor.jl index de6faf811..3d6363e66 100644 --- a/test/taylor.jl +++ b/test/taylor.jl @@ -33,3 +33,13 @@ eq = taylor(eq, x, 0:7) eqs = taylor_coeff(eq, x) # should automatically expand to 7th order @test length(eqs) == 7+1 && all(isequal(eq.lhs, eq.rhs) for eq in eqs) + +# expand quintic equation around x=1 +@variables ϵ +x_series = series(x, ϵ, 0:3) +x_coeffs = taylor_coeff(x_series, ϵ) +eq = x^5 + ϵ*x ~ 1 +eqs = taylor_coeff(substitute(eq, x => x_series), ϵ, 0:3) +sol = x_coeffs .=> [1, -1//5, -1//25, -1//125] # e.g. https://ekamperi.github.io/mathematics/2020/06/21/perturbation-theory.html#MathJax-Element-39-Frame +eqs = substitute(eqs, Dict(sol)) +@test all(isequal(eq.lhs, eq.rhs) for eq in eqs) From 27b33366e09e3c716c923f06882fce1e933f7212 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Wed, 9 Oct 2024 22:35:47 +0200 Subject: [PATCH 09/17] Add Taylor series documentation section --- docs/make.jl | 1 + docs/src/manual/taylor.md | 7 +++++++ 2 files changed, 8 insertions(+) create mode 100644 docs/src/manual/taylor.md diff --git a/docs/make.jl b/docs/make.jl index c99454081..22a8c6dd7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -41,6 +41,7 @@ makedocs( "manual/variables.md", "manual/expression_manipulation.md", "manual/derivatives.md", + "manual/taylor.md", "manual/groebner.md", "manual/solver.md", "manual/arrays.md", diff --git a/docs/src/manual/taylor.md b/docs/src/manual/taylor.md new file mode 100644 index 000000000..775e9420e --- /dev/null +++ b/docs/src/manual/taylor.md @@ -0,0 +1,7 @@ +# Taylor Series + +```@docs +series +taylor +taylor_coeff +``` From 305afc4310901bbf58c8e58329b7d4e242889989 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Wed, 9 Oct 2024 22:43:14 +0200 Subject: [PATCH 10/17] Add examples to series() and taylor() docstrings --- src/taylor.jl | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/src/taylor.jl b/src/taylor.jl index 34150cbf3..9b93a0ab8 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -2,6 +2,19 @@ series(y, x, [x0=0,] ns; name = nameof(y)) Expand the variable `y` in a power series in the variable `x` around `x0` to orders `ns`. + +Examples +======== + +```julia +julia> @variables z ϵ +2-element Vector{Num}: + z + ϵ + +julia> series(z, ϵ, 2, 0:3) +z[0] + z[1]*(-2 + ϵ) + z[2]*((-2 + ϵ)^2) + z[3]*((-2 + ϵ)^3) +``` """ function series(y, x, x0, ns; name = nameof(y)) c, = @variables $name[ns] @@ -15,6 +28,21 @@ end taylor_coeff(f, x[, n]; rationalize=true) Calculate the `n`-th order coefficient(s) in the Taylor series of `f` around `x = 0`. + +Examples +======== +```julia +julia> @variables x y +2-element Vector{Num}: + x + y + +julia> taylor_coeff(series(y, x, 0:5), x, 0:2:4) +3-element Vector{Num}: + y[0] + y[2] + y[4] +``` """ function taylor_coeff(f, x, n = missing; rationalize=true) if n isa AbstractArray @@ -70,6 +98,9 @@ julia> taylor(exp(x), x, 0:3; rationalize=false) julia> taylor(√(x), x, 1, 0:3) 1 + (1//2)*(-1 + x) - (1//8)*((-1 + x)^2) + (1//16)*((-1 + x)^3) + +julia> isequal(taylor(exp(im*x), x, 0:5), taylor(exp(im*x), x, 0:5)) +true ``` """ function taylor(f, x, ns; kwargs...) From 5e3f60d299a6490d0c20b86c64262a33757032fe Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Wed, 9 Oct 2024 22:49:33 +0200 Subject: [PATCH 11/17] Taylor expand systems of equations --- src/taylor.jl | 4 +++- test/taylor.jl | 3 +++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/taylor.jl b/src/taylor.jl index 9b93a0ab8..aa9ebc94c 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -104,7 +104,9 @@ true ``` """ function taylor(f, x, ns; kwargs...) - if f isa Equation + if f isa AbstractArray + return taylor.(f, Ref(x), Ref(ns); kwargs...) + elseif f isa Equation return taylor(f.lhs, x, ns; kwargs...) ~ taylor(f.rhs, x, ns; kwargs...) end diff --git a/test/taylor.jl b/test/taylor.jl index 3d6363e66..6f995dda2 100644 --- a/test/taylor.jl +++ b/test/taylor.jl @@ -43,3 +43,6 @@ eqs = taylor_coeff(substitute(eq, x => x_series), ϵ, 0:3) sol = x_coeffs .=> [1, -1//5, -1//25, -1//125] # e.g. https://ekamperi.github.io/mathematics/2020/06/21/perturbation-theory.html#MathJax-Element-39-Frame eqs = substitute(eqs, Dict(sol)) @test all(isequal(eq.lhs, eq.rhs) for eq in eqs) + +# system of equations +@test taylor(exp(im*x) ~ 0, x, 0:5) == taylor([cos(x) ~ 0, sin(x) ~ 0], x, 0:5) From 066a0358c25efd4bb74519472f4b037533da0e52 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Wed, 9 Oct 2024 22:50:10 +0200 Subject: [PATCH 12/17] Rewrite symbolic-numeric perturbation example --- docs/src/examples/perturbation.md | 281 +++++++++--------------------- 1 file changed, 86 insertions(+), 195 deletions(-) diff --git a/docs/src/examples/perturbation.md b/docs/src/examples/perturbation.md index c57ff8bcb..c6e3f74fe 100644 --- a/docs/src/examples/perturbation.md +++ b/docs/src/examples/perturbation.md @@ -1,250 +1,141 @@ # [Mixed Symbolic-Numeric Perturbation Theory](@id perturb_alg) -## Background +[**Symbolics.jl**](https://github.com/JuliaSymbolics/Symbolics.jl) is a fast and modern Computer Algebra System (CAS) written in Julia. It is an integral part of the [SciML](https://sciml.ai/) ecosystem of differential equation solvers and scientific machine learning packages. While **Symbolics.jl** is primarily designed for modern scientific computing (e.g. automatic differentiation and machine learning), it is also a powerful CAS that can be used for *classic* scientific computing. One such application is using *perturbation theory* to solve algebraic and differential equations. -[**Symbolics.jl**](https://github.com/JuliaSymbolics/Symbolics.jl) is a fast and modern Computer Algebra System (CAS) written in the Julia Programming Language. It is an integral part of the [SciML](https://sciml.ai/) ecosystem of differential equation solvers and scientific machine learning packages. While **Symbolics.jl** is primarily designed for modern scientific computing (e.g., auto-differentiation, machine learning), it is a powerful CAS that can also be useful for *classic* scientific computing. One such application is using the *perturbation* theory to solve algebraic and differential equations. +Perturbation methods are a collection of techniques to solve hard problems that generally don't have a closed solution, but depend on a tunable parameter and have closed or easy solutions for some values of this parameter. The main idea is to assume a solution that is a power series in the tunable parameter (say $ϵ$), such that $ϵ = 0$ corresponds to an easy solution, and then solve iteratively for higher-order corrections. -Perturbation methods are a collection of techniques to solve intractable problems that generally don't have a closed solution but depend on a tunable parameter and have closed or easy solutions for some values of the parameter. The main idea is to assume a solution as a power series in the tunable parameter (say $ϵ$), such that $ϵ = 0$ corresponds to an easy solution. +The hallmark of perturbation theory is the generation of long and convoluted intermediate equations by this process. These are subjected to algorithmic and mechanical manipulations, which makes perturbation methods an excellent fit for a CAS. In fact, CAS software have been used for perturbation calculations since the early 1970s. -We will discuss the general steps of the perturbation methods to solve algebraic (this tutorial) and [differential equations](https://docs.sciml.ai/ModelingToolkit/stable/examples/perturbation/) +This tutorial shows how to mix symbolic manipulations and numerical methods to solve algebraic equations with perturbation theory. [Another tutorial applies it to differential equations](https://docs.sciml.ai/ModelingToolkit/stable/examples/perturbation/). -The hallmark of the perturbation method is the generation of long and convoluted intermediate equations, which are subjected to algorithmic and mechanical manipulations. Therefore, these problems are well suited for CAS. In fact, CAS software packages have been used to help with the perturbation calculations since the early 1970s. - -In this tutorial, our goal is to show how to use a mix of symbolic manipulations (**Symbolics.jl**) and numerical methods to solve simple perturbation problems. - -## Solving the Quintic - -We start with the “hello world!” analog of the perturbation problems, solving the quintic (fifth-order) equations. We want to find a real valued $x$ such that $x^5 + x = 1$. According to the Abel's theorem, a general quintic equation does not have a closed form solution. Of course, we can easily solve this equation numerically; for example, by using the Newton's method. We use the following implementation of the Newton's method: +## Solving the quintic equation +The “hello world!” analog of perturbation problems is to find a real solution $x$ to the quintic (fifth-order) equation ```@example perturb -using Symbolics, SymbolicUtils - -function solve_newton(f, x, x₀; abstol=1e-8, maxiter=50) - xₙ = Float64(x₀) - fₙ₊₁ = x - f / Symbolics.derivative(f, x) +using Symbolics +@variables x +quintic = x^5 + x ~ 1 +``` +According to Abel's theorem, a general quintic equation does not have a closed form solution. But we can easily solve it numerically using Newton's method (here implemented for simplicity, and not performance): +```@example perturb +function solve_newton(eq, x, x₀; abstol=1e-8, maxiters=50) + # symbolic expressions for f(x) and f′(x) + f = eq.lhs - eq.rhs # want to find root of f(x) + f′ = Symbolics.derivative(f, x) - for i = 1:maxiter - xₙ₊₁ = substitute(fₙ₊₁, Dict(x => xₙ)) + xₙ = x₀ # numerical value of the initial guess + for i = 1:maxiters + # calculate new guess by numerically evaluating symbolic expression at previous guess + xₙ₊₁ = substitute(x - f / f′, x => xₙ) if abs(xₙ₊₁ - xₙ) < abstol - return xₙ₊₁ + return xₙ₊₁ # converged else xₙ = xₙ₊₁ end end - return xₙ₊₁ + error("Newton's method failed to converge") end -``` - -In this code, `Symbolics.derivative(eq, x)` does exactly what it names implies: it calculates the symbolic derivative of `eq` (a **Symbolics.jl** expression) with respect to `x` (a **Symbolics.jl** variable). We use `Symbolics.substitute(eq, D)` to evaluate the update formula by substituting variables or sub-expressions (defined in a dictionary `D`) in `eq`. It should be noted that `substitute` is the workhorse of our code and will be used multiple times in the rest of these tutorials. `solve_newton` is written with simplicity and clarity in mind, and not performance. - -Let's go back to our quintic. We can define a Symbolics variable as `@variables x` and then solve the equation `solve_newton(x^5 + x - 1, x, 1.0)` (here, `x₀ = 1.0` is our first guess). The answer is 0.7549. Now, let's see how we can solve the same problem using the perturbation methods. - -We introduce a tuning parameter $\epsilon$ into our equation: $x^5 + \epsilon x = 1$. If $\epsilon = 1$, we get our original problem. For $\epsilon = 0$, the problem transforms to an easy one: $x^5 = 1$ which has an exact real solution $x = 1$ (and four complex solutions which we ignore here). We expand $x$ as a power series on $\epsilon$: - -```math - x(\epsilon) = a_0 + a_1 \epsilon + a_2 \epsilon^2 + O(\epsilon^3) -``` - -$a_0$ is the solution of the easy equation, therefore $a_0 = 1$. Substituting into the original problem, - -```math - (a_0 + a_1 \epsilon + a_2 \epsilon^2)^5 + \epsilon (a_0 + a_1 \epsilon + a_2 \epsilon^2) - 1 = 0 -``` - -Expanding the equations, we get - -```math - \epsilon (1 + 5 a_1) + \epsilon^2 (a_1 + 5 a_2 + 10 a_1^2) + 𝑂(\epsilon^3) = 0 -``` - -This equation should hold for each power of $\epsilon$. Therefore, - -```math - 1 + 5 a_1 = 0 -``` - -and -```math - a_1 + 5 a_2 + 10 a_1^2 = 0 +x_newton = solve_newton(quintic, x, 1.0) +println("Newton's method solution: x = ", x_newton) ``` -This system of equations does not initially seem to be linear because of the presence of terms like $10 a_1^2$, but upon closer inspection is found to be linear (this is a feature of the perturbation methods). In addition, the system is in a triangular form, meaning the first equation depends only on $a_1$, the second one on $a_1$ and $a_2$, such that we can replace the result of $a_1$ from the first one into the second equation and remove the non-linear term. We solve the first equation to get $a_1 = -\frac{1}{5}$. Substituting in the second one and solve for $a_2$: - -```math - a_2 = \frac{(-\frac{1}{5} + 10(-(\frac{1}{5})²)}{5} = -\frac{1}{25} -``` - -Finally, - -```math - x(\epsilon) = 1 - \frac{\epsilon}{5} - \frac{\epsilon^2}{25} + O(\epsilon^3) -``` - -Solving the original problem, $x(1) = 0.76$, compared to 0.7548 calculated numerically. We can improve the accuracy by including more terms in the expansion of $x$. However, the calculations, while straightforward, become messy and intractable to do manually very quickly. This is why a CAS is very helpful to solve perturbation problems. - -Now, let's see how we can do these calculations in Julia. Let $n$ be the order of the expansion. We start by defining the symbolic variables: - +To solve the same problem with perturbation theory, we must introduce an expansion variable $ϵ$ in the equation: ```@example perturb -n = 2 -@variables ϵ a[1:n] +@variables ϵ # expansion variable +quintic = x^5 + ϵ*x ~ 1 ``` - -Then, we define - +If $ϵ = 1$, we get our original problem. With $ϵ = 0$, the problem transforms to the easy quintic equation $x^5 = 1$ with the trivial real solution $x = 1$ (and four complex solutions which we ignore). Next, expand $x$ as a power series in $ϵ$: ```@example perturb -x = 1 + a[1]*ϵ + a[2]*ϵ^2 +x_taylor = series(x, ϵ, 0:3) # expand x to third order ``` - -The next step is to substitute `x` in the problem equation - +Then insert this into the quintic equation and expand it, too, to the same order: ```@example perturb - eq = x^5 + ϵ*x - 1 +quintic_taylor = substitute(quintic, x => x_taylor) +quintic_taylor = taylor(quintic_taylor, ϵ, 0:3) ``` - -The expanded form of `eq` is - +This equation must hold for each power of $ϵ$, so we can separate it into one equation per order: ```@example perturb -expand(eq) +quintic_eqs = taylor_coeff(quintic_taylor, ϵ) ``` - -We need a way to get the coefficients of different powers of `ϵ`. Function `collect_powers(eq, x, ns)` returns the powers of variable `x` in expression `eq`. Argument `ns` is the range of the powers. - +These equations show three important features of perturbation theory: +1. The $0$-th order equation is *trivial* in $x_0$: here $x_0^5 = 1$ has the trivial real solution $x_0 = 1$. +2. The $n$-th order equation is *linear* in $x_n$ (except the trivial $0$-th order equation). +3. The equations are *triangular* in $x_n$: the $n$-th order equation can be solved for $x_n$ given only $x_m$ for $m x_taylor) # expand equation in taylor series + eqs = taylor_coeff(eq_taylor, ϵ, 0:order) # separate into order-by-order equations + sol = Dict(x_coeffs[1] => x₀) # store solutions in a symbolic-numeric map -```@example perturb -eqs = collect_powers(eq, ϵ, 1:2) -``` - -Having the coefficients of the powers of `ϵ`, we can set each equation in `eqs` to 0 (remember, we rearrange the problem such that `eq` is 0) and solve the system of linear equations to find the numerical values of the coefficients. **Symbolics.jl** has a function `symbolic_linear_solve` that can solve systems of linear equations. However, the presence of higher-order terms in `eqs` prevents `symbolic_linear_solve(eqs, a)` from workings properly. Instead, we can exploit the fact that our system is in a triangular form and start by solving `eqs[1]` for `a₁` and then substitute this in `eqs[2]` and solve for `a₂`, and so on. This *cascading* process is done by function `solve_coef(eqs, ps)`: - -```@example perturb -function solve_coef(eqs, ps) - vals = Dict() + # verify that x₀ is a solution of the 0-th order equation + eq0 = substitute(eqs[1], sol) + if !isequal(eq0.lhs, eq0.rhs) + error("$sol is not a 0-th order solution of $(eqs[1])") + end - for i = 1:length(ps) - eq = substitute(eqs[i], vals) - vals[ps[i]] = symbolic_linear_solve(eq, ps[i]) + # solve higher-order equations order-by-order + for i in 2:length(eqs) + x_coeff = Symbolics.symbolic_linear_solve(eqs[i], x_coeffs[i]) # solve linear n-th order equation for x_n + x_coeff = substitute(x_coeff, sol) # substitute lower-order solutions to get numerical value + sol = merge(sol, Dict(x_coeffs[i] => x_coeff)) # store solution end - vals -end -``` -Here, `eqs` is an array of expressions (assumed to be equal to 0) and `ps` is an array of variables. The result is a dictionary of *variable* => *value* pairs. We apply `solve_coef` to `eqs` to get the numerical values of the parameters: + return substitute(x_taylor, sol) # evalaute series with solved coefficients +end -```@example perturb -vals = solve_coef(eqs, a) +x_pert = solve_perturbed(quintic, x, 1, ϵ, 7) ``` - -Finally, we substitute back the values of `a` in the definition of `x` as a function of `𝜀`. Note that `𝜀` is a number (usually Float64), whereas `ϵ` is a symbolic variable. - +The $n$-th order solution of our original quintic equation is the sum up to the $ϵ^n$-th order term, evaluated at $ϵ=1$: ```@example perturb -X = 𝜀 -> 1 + vals[a[1]]*𝜀 + vals[a[2]]*𝜀^2 -``` - -Therefore, the solution to our original problem becomes `X(1)`, which is equal to 0.76. We can use larger values of `n` to improve the accuracy of estimations. - -| n | x | -|---|----------------| -|1 |0.8 | -|2 |0.76| -|3 |0.752| -|4 |0.752| -|5 |0.7533| -|6 |0.7543| -|7 |0.7548| -|8 |0.7550| - -Remember, the numerical value is 0.7549. The two functions `collect_powers` and `solve_coef(eqs, a)` are used in all the examples in this and the next tutorial. - -## Solving the Kepler's Equation - -Historically, the perturbation methods were first invented to solve orbital calculations of the Moon and the planets. In homage to this history, our second example has a celestial theme. Our goal is solving the Kepler's equation: - -```math - E - e\sin(E) = M +for n in 0:7 + println("$n-th order solution: x = ", substitute(taylor(x_pert, ϵ, 0:n), ϵ => 1.0)) +end ``` +This is close to the solution from Newton's method! -where $e$ is the *eccentricity* of the elliptical orbit, $M$ is the *mean anomaly*, and $E$ (unknown) is the *eccentric anomaly* (the angle between the position of a planet in an elliptical orbit and the point of periapsis). This equation is central to solving two-body Keplerian orbits. - -Similar to the first example, it is easy to solve this problem using the Newton's method. For example, let $e = 0.01671$ (the eccentricity of the Earth) and $M = \pi/2$. We have `solve_newton(x - e*sin(x) - M, x, M)` equals to 1.5875 (compared to π/2 = 1.5708). Now, we try to solve the same problem using the perturbation techniques (see function `test_kepler`). - -For $e = 0$, we get $E = M$. Therefore, we can use $e$ as our perturbation parameter. For consistency with other problems, we also rename $e$ to $\epsilon$ and $E$ to $x$. - -From here on, we use the helper function `def_taylor` to define Taylor's series by calling it as `x = def_taylor(ϵ, a, 1)`, where the arguments are, respectively, the perturbation variable, which is an array of coefficients (starting from the coefficient of $\epsilon^1$), and an optional constant term. +## Solving Kepler's Equation +Historically, perturbation methods were first invented to calculate orbits of the Moon and the planets. In homage to this history, our second example is to solve [Kepler's equation](https://en.wikipedia.org/wiki/Kepler's_equation), which is central to solving two-body Keplerian orbits: ```@example perturb -def_taylor(x, ps) = sum([a*x^i for (i,a) in enumerate(ps)]) -def_taylor(x, ps, p₀) = p₀ + def_taylor(x, ps) +@variables e E M +kepler = E - e * sin(E) ~ M ``` - -We start by defining the variables (assuming `n = 3`): - +We want to solve for the *eccentric anomaly* $E$ given the *eccentricity* $e$ and *mean anomaly* $M$. +This is also easy with Newton's method. With Earth's eccentricity $e = 0.01671$ and $M = \pi/2$: ```@example perturb -n = 3 -@variables ϵ M a[1:n] -x = def_taylor(ϵ, a, M) +vals_earth = Dict(e => 0.01671, M => π/2) +E_newton = solve_newton(substitute(kepler, vals_earth), E, π/2) +println("Newton's method solution: E = ", E_newton) ``` -We further simplify by substituting `sin` with its power series using the `expand_sin` helper function: - +Next, let us solve the same problem with our perturbative solver. It is most common to expand Kepler's equation in $M$ (the trivial solution when $M=0$ is $E=0$): ```@example perturb -expand_sin(x, n) = sum([(isodd(k) ? -1 : 1)*(-x)^(2k-1)/factorial(2k-1) for k=1:n]) +E_pert = solve_perturbed(kepler, E, 0, M, 5) ``` - -To test, - +Numerically, this gives almost the same answer as Newton's method: ```@example perturb -expand_sin(0.1, 10) ≈ sin(0.1) +for n in 0:5 + println("$n-th order solution: E = ", substitute(taylor(E_pert, M, 0:n), vals_earth)) +end ``` - -The problem equation is - +But unlike Newtons method, perturbation theory also gives us the power to work with the full *symbolic* series solution for $E$ (*before* numbers for $e$ and $M$ are inserted). Our series matches [this result from Wikipedia](https://en.wikipedia.org/wiki/Kepler%27s_equation#Inverse_Kepler_equation): ```@example perturb -eq = x - ϵ * expand_sin(x, n) - M +E_wiki = 1/(1-e)*M - e/(1-e)^4*M^3/factorial(3) + (9e^2+e)/(1-e)^7*M^5/factorial(5) ``` -We follow the same process as the first example. We collect the coefficients of the powers of `ϵ` - +Alternatively, we can expand Kepler's equation in $e$ instead of $M$ (the trivial solution when $e = 0$ is $E=M$): ```@example perturb -eqs = collect_powers(eq, ϵ, 1:n) +E_pert′ = solve_perturbed(kepler, E, M, e, 5) ``` - -and then solve for `a`: - +We can expand the trigonometric functions in $M$: ```@example perturb -vals = solve_coef(eqs, a) +E_pert′ = taylor(E_pert′, M, 0:5) ``` - -Finally, we substitute `vals` back in `x`: - +Up to order $e^5 M^5$, we see that this two-variable $(e,M)$-series also matches the result from Wikipedia: ```@example perturb -x′ = substitute(x, vals) -X = (𝜀, 𝑀) -> substitute(x′, Dict(ϵ => 𝜀, M => 𝑀)) -X(0.01671, π/2) -``` - -The result is 1.5876, compared to the numerical value of 1.5875. It is customary to order `X` based on the powers of `𝑀` instead of `𝜀`. We can calculate this series as `collect_powers(x′, M, 0:5)`. The result (after manual cleanup) is - +E_wiki′ = taylor(taylor(E_wiki, e, 0:5), M, 0:5) ``` -(1 + 𝜀 + 𝜀^2 + 𝜀^3)*𝑀 -- (𝜀 + 4*𝜀^2 + 10*𝜀^3)*𝑀^3/6 -+ (𝜀 + 16*𝜀^2 + 91*𝜀^3)*𝑀^5/120 -``` - -Comparing the formula to the one for 𝐸 in the [Wikipedia article on the Kepler's equation](https://en.wikipedia.org/wiki/Kepler%27s_equation): - -```math - E = \frac{1}{1-\epsilon}M - -\frac{\epsilon}{(1-\epsilon)^4} \frac{M^3}{3!} + \frac{(9\epsilon^2 - + \epsilon)}{(1-\epsilon)^7}\frac{M^5}{5!}\cdots -``` - -The first deviation is in the coefficient of $\epsilon^3 M^5$. From a31eebbdca9e8db0f9b286aae22e9b32d9ba3062 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Thu, 10 Oct 2024 13:38:10 +0200 Subject: [PATCH 13/17] Move perturbation example to tutorials (instead of having one lonely example) --- docs/make.jl | 6 ++---- docs/src/{examples => tutorials}/perturbation.md | 0 2 files changed, 2 insertions(+), 4 deletions(-) rename docs/src/{examples => tutorials}/perturbation.md (100%) diff --git a/docs/make.jl b/docs/make.jl index 22a8c6dd7..855ff3421 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -32,10 +32,8 @@ makedocs( "getting_started.md", "Tutorials" => Any[ "tutorials/auto_parallel.md", - "tutorials/converting_to_C.md" - ], - "Examples" => Any[ - "examples/perturbation.md" + "tutorials/converting_to_C.md", + "tutorials/perturbation.md" ], "Manual" => Any[ "manual/variables.md", diff --git a/docs/src/examples/perturbation.md b/docs/src/tutorials/perturbation.md similarity index 100% rename from docs/src/examples/perturbation.md rename to docs/src/tutorials/perturbation.md From 779810e3416293259674b08277bedda0ee25b107 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Thu, 10 Oct 2024 13:49:33 +0200 Subject: [PATCH 14/17] Link tutorial from Taylor series manual page --- docs/src/manual/taylor.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/manual/taylor.md b/docs/src/manual/taylor.md index 775e9420e..0139e5d43 100644 --- a/docs/src/manual/taylor.md +++ b/docs/src/manual/taylor.md @@ -1,5 +1,7 @@ # Taylor Series +For a real example of how to use the Taylor series functionality, see [this tutorial](@ref perturb_alg). + ```@docs series taylor From 5ade47f9123773ab76ab43882c1b42b1034a51ac Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Thu, 10 Oct 2024 14:24:59 +0200 Subject: [PATCH 15/17] Generate fractions in Taylor series coefficients a bit more robustly --- docs/src/tutorials/perturbation.md | 2 +- src/taylor.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/perturbation.md b/docs/src/tutorials/perturbation.md index c6e3f74fe..02d2d0122 100644 --- a/docs/src/tutorials/perturbation.md +++ b/docs/src/tutorials/perturbation.md @@ -79,8 +79,8 @@ function solve_perturbed(eq, x, x₀, ϵ, order) # solve higher-order equations order-by-order for i in 2:length(eqs) + eqs[i] = substitute(eqs[i], sol) # substitute lower-order solutions x_coeff = Symbolics.symbolic_linear_solve(eqs[i], x_coeffs[i]) # solve linear n-th order equation for x_n - x_coeff = substitute(x_coeff, sol) # substitute lower-order solutions to get numerical value sol = merge(sol, Dict(x_coeffs[i] => x_coeff)) # store solution end diff --git a/src/taylor.jl b/src/taylor.jl index aa9ebc94c..7218543be 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -72,7 +72,7 @@ function taylor_coeff(f, x, n = missing; rationalize=true) if rationalize && unwrap(c) isa Number # TODO: make rational coefficients "organically" and not using rationalize (see https://github.com/JuliaSymbolics/Symbolics.jl/issues/1299) c = unwrap(c) - c = isinteger(c) ? Int(c) : Base.rationalize(c) # convert non-integer floats to rational numbers; avoid name clash between rationalize and Base.rationalize() + c = Base.rationalize(c) # convert integers/floats to rational numbers; avoid name clash between rationalize and Base.rationalize() end return c end From 5f01f6a26eadd328b5067caa5198e0f8c7c05752 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Mon, 4 Nov 2024 21:07:12 +0100 Subject: [PATCH 16/17] Make perturbation tutorial more understandable --- docs/src/tutorials/perturbation.md | 83 +++++++++++++++++------------- 1 file changed, 48 insertions(+), 35 deletions(-) diff --git a/docs/src/tutorials/perturbation.md b/docs/src/tutorials/perturbation.md index 02d2d0122..b686b25da 100644 --- a/docs/src/tutorials/perturbation.md +++ b/docs/src/tutorials/perturbation.md @@ -47,52 +47,55 @@ quintic = x^5 + ϵ*x ~ 1 ``` If $ϵ = 1$, we get our original problem. With $ϵ = 0$, the problem transforms to the easy quintic equation $x^5 = 1$ with the trivial real solution $x = 1$ (and four complex solutions which we ignore). Next, expand $x$ as a power series in $ϵ$: ```@example perturb -x_taylor = series(x, ϵ, 0:3) # expand x to third order +x_taylor = series(x, ϵ, 0:7) # expand x in a power series in ϵ +x_coeffs = taylor_coeff(x_taylor, ϵ) # TODO: get coefficients at series creation +x_taylor ``` Then insert this into the quintic equation and expand it, too, to the same order: ```@example perturb quintic_taylor = substitute(quintic, x => x_taylor) -quintic_taylor = taylor(quintic_taylor, ϵ, 0:3) +quintic_taylor = taylor(quintic_taylor, ϵ, 0:7) ``` -This equation must hold for each power of $ϵ$, so we can separate it into one equation per order: +This messy equation must hold for each power of $ϵ$, so we can separate it into one nicer equation per order: ```@example perturb quintic_eqs = taylor_coeff(quintic_taylor, ϵ) +quintic_eqs[1:5] # for readability, show only 5 shortest equations ``` These equations show three important features of perturbation theory: 1. The $0$-th order equation is *trivial* in $x_0$: here $x_0^5 = 1$ has the trivial real solution $x_0 = 1$. 2. The $n$-th order equation is *linear* in $x_n$ (except the trivial $0$-th order equation). 3. The equations are *triangular* in $x_n$: the $n$-th order equation can be solved for $x_n$ given only $x_m$ for $m x_taylor) # expand equation in taylor series - eqs = taylor_coeff(eq_taylor, ϵ, 0:order) # separate into order-by-order equations - sol = Dict(x_coeffs[1] => x₀) # store solutions in a symbolic-numeric map +function solve_cascade(eqs, xs, x₀, ϵ) + sol = Dict(xs[begin] => x₀) # store solutions in a map - # verify that x₀ is a solution of the 0-th order equation + # verify that x₀ is a solution of the first equation eq0 = substitute(eqs[1], sol) - if !isequal(eq0.lhs, eq0.rhs) - error("$sol is not a 0-th order solution of $(eqs[1])") - end + isequal(eq0.lhs, eq0.rhs) || error("$sol does not solve $(eqs[1])") - # solve higher-order equations order-by-order - for i in 2:length(eqs) - eqs[i] = substitute(eqs[i], sol) # substitute lower-order solutions - x_coeff = Symbolics.symbolic_linear_solve(eqs[i], x_coeffs[i]) # solve linear n-th order equation for x_n - sol = merge(sol, Dict(x_coeffs[i] => x_coeff)) # store solution + # solve remaining equations sequentially + for i in 2:lastindex(eqs) + eq = substitute(eqs[i], sol) # insert previous solutions + x = Symbolics.symbolic_linear_solve(eq, xs[i]) # solve current equation + sol = merge(sol, Dict(xs[i] => x)) # store solution end - return substitute(x_taylor, sol) # evalaute series with solved coefficients + return sol end - -x_pert = solve_perturbed(quintic, x, 1, ϵ, 7) +``` +Let us solve our order-separated quintics for the coefficients, and substitute them into the full series for $x$: +```@example perturb +x_coeffs_sol = solve_cascade(quintic_eqs, x_coeffs, 1, ϵ) +x_pert = substitute(x_taylor, x_coeffs_sol) ``` The $n$-th order solution of our original quintic equation is the sum up to the $ϵ^n$-th order term, evaluated at $ϵ=1$: ```@example perturb for n in 0:7 - println("$n-th order solution: x = ", substitute(taylor(x_pert, ϵ, 0:n), ϵ => 1.0)) + x_pert_sol = substitute(taylor(x_pert, ϵ, 0:n), ϵ => 1) + println("$n-th order solution: x = $x_pert_sol = $(x_pert_sol * 1.0)") end ``` This is close to the solution from Newton's method! @@ -112,30 +115,40 @@ E_newton = solve_newton(substitute(kepler, vals_earth), E, π/2) println("Newton's method solution: E = ", E_newton) ``` -Next, let us solve the same problem with our perturbative solver. It is most common to expand Kepler's equation in $M$ (the trivial solution when $M=0$ is $E=0$): +Next, let us solve the same problem with the perturbative method. It is most common to expand $E$ as a series in $M$. Repeating the procedure from the quintic example, we get these equations: ```@example perturb -E_pert = solve_perturbed(kepler, E, 0, M, 5) +E_taylor = series(E, M, 0:5) +E_coeffs = taylor_coeff(E_taylor, M) # TODO: get coefficients at series creation +kepler_eqs = taylor_coeff(substitute(kepler, E => E_taylor), M, 0:5) +kepler_eqs[1:4] # for readability ``` -Numerically, this gives almost the same answer as Newton's method: +The trivial $0$-th order solution (when $M=0$) is $E_0=0$. This gives this full perturbative solution: +```@example perturb +E_coeffs_sol = solve_cascade(kepler_eqs, E_coeffs, 0, M) +E_pert = substitute(E_taylor, E_coeffs_sol) +``` + +Numerically, the result again converges to that from Newton's method: ```@example perturb for n in 0:5 println("$n-th order solution: E = ", substitute(taylor(E_pert, M, 0:n), vals_earth)) end ``` -But unlike Newtons method, perturbation theory also gives us the power to work with the full *symbolic* series solution for $E$ (*before* numbers for $e$ and $M$ are inserted). Our series matches [this result from Wikipedia](https://en.wikipedia.org/wiki/Kepler%27s_equation#Inverse_Kepler_equation): +But unlike Newtons method, this example shows how perturbation theory also gives us the powerful *symbolic* series solution for $E$ (*before* numbers for $e$ and $M$ are inserted). Our series matches [this result from Wikipedia](https://en.wikipedia.org/wiki/Kepler%27s_equation#Inverse_Kepler_equation): ```@example perturb E_wiki = 1/(1-e)*M - e/(1-e)^4*M^3/factorial(3) + (9e^2+e)/(1-e)^7*M^5/factorial(5) ``` -Alternatively, we can expand Kepler's equation in $e$ instead of $M$ (the trivial solution when $e = 0$ is $E=M$): -```@example perturb -E_pert′ = solve_perturbed(kepler, E, M, e, 5) -``` -We can expand the trigonometric functions in $M$: +Alternatively, we can expand $E$ in $e$ instead of $M$, giving the solution (the trivial solution when $e = 0$ is $E_0=M$): ```@example perturb -E_pert′ = taylor(E_pert′, M, 0:5) +E_taylor′ = series(E, e, 0:5) +E_coeffs′ = taylor_coeff(E_taylor′, e) # TODO: get at creation +kepler_eqs′ = taylor_coeff(substitute(kepler, E => E_taylor′), e, 0:5) +E_coeffs_sol′ = solve_cascade(kepler_eqs′, E_coeffs′, M, e) +E_pert′ = substitute(E_taylor′, E_coeffs_sol′) ``` -Up to order $e^5 M^5$, we see that this two-variable $(e,M)$-series also matches the result from Wikipedia: +This looks very different from our first series `E_pert`. If they are the same, we should get $0$ if we subtract and expand both as multivariate Taylor series in $(e,M)$. Indeed: ```@example perturb -E_wiki′ = taylor(taylor(E_wiki, e, 0:5), M, 0:5) +@assert taylor(taylor(E_pert′ - E_pert, e, 0:4), M, 0:4) == 0 # use this as a test # hide +taylor(taylor(E_pert′ - E_pert, e, 0:4), M, 0:4) ``` From 7b0546ff371c625c3f1e16e91bc6e0e03e658a99 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Mon, 4 Nov 2024 23:32:43 +0100 Subject: [PATCH 17/17] Let user pass coefficients to series(); show both this and the version with automatic coefficients in tutorial --- docs/src/tutorials/perturbation.md | 20 +++++++-------- src/taylor.jl | 40 +++++++++++++++++++++--------- test/taylor.jl | 18 ++++++++++++++ 3 files changed, 56 insertions(+), 22 deletions(-) diff --git a/docs/src/tutorials/perturbation.md b/docs/src/tutorials/perturbation.md index b686b25da..30963997e 100644 --- a/docs/src/tutorials/perturbation.md +++ b/docs/src/tutorials/perturbation.md @@ -47,9 +47,8 @@ quintic = x^5 + ϵ*x ~ 1 ``` If $ϵ = 1$, we get our original problem. With $ϵ = 0$, the problem transforms to the easy quintic equation $x^5 = 1$ with the trivial real solution $x = 1$ (and four complex solutions which we ignore). Next, expand $x$ as a power series in $ϵ$: ```@example perturb -x_taylor = series(x, ϵ, 0:7) # expand x in a power series in ϵ -x_coeffs = taylor_coeff(x_taylor, ϵ) # TODO: get coefficients at series creation -x_taylor +x_coeffs, = @variables a[0:7] # create Taylor series coefficients +x_taylor = series(x_coeffs, ϵ) # expand x in a power series in ϵ ``` Then insert this into the quintic equation and expand it, too, to the same order: ```@example perturb @@ -58,7 +57,7 @@ quintic_taylor = taylor(quintic_taylor, ϵ, 0:7) ``` This messy equation must hold for each power of $ϵ$, so we can separate it into one nicer equation per order: ```@example perturb -quintic_eqs = taylor_coeff(quintic_taylor, ϵ) +quintic_eqs = taylor_coeff(quintic_taylor, ϵ, 0:7) quintic_eqs[1:5] # for readability, show only 5 shortest equations ``` These equations show three important features of perturbation theory: @@ -77,10 +76,11 @@ function solve_cascade(eqs, xs, x₀, ϵ) isequal(eq0.lhs, eq0.rhs) || error("$sol does not solve $(eqs[1])") # solve remaining equations sequentially - for i in 2:lastindex(eqs) + for i in 2:length(eqs) eq = substitute(eqs[i], sol) # insert previous solutions - x = Symbolics.symbolic_linear_solve(eq, xs[i]) # solve current equation - sol = merge(sol, Dict(xs[i] => x)) # store solution + x = xs[begin+i-1] # need not be 1-indexed + xsol = Symbolics.symbolic_linear_solve(eq, x) # solve current equation + sol = merge(sol, Dict(x => xsol)) # store solution end return sol @@ -117,8 +117,8 @@ println("Newton's method solution: E = ", E_newton) Next, let us solve the same problem with the perturbative method. It is most common to expand $E$ as a series in $M$. Repeating the procedure from the quintic example, we get these equations: ```@example perturb -E_taylor = series(E, M, 0:5) -E_coeffs = taylor_coeff(E_taylor, M) # TODO: get coefficients at series creation +E_taylor = series(E, M, 0:5) # this auto-creates coefficients E[0:5] +E_coeffs = taylor_coeff(E_taylor, M) # get a handle to the coefficients kepler_eqs = taylor_coeff(substitute(kepler, E => E_taylor), M, 0:5) kepler_eqs[1:4] # for readability ``` @@ -142,7 +142,7 @@ E_wiki = 1/(1-e)*M - e/(1-e)^4*M^3/factorial(3) + (9e^2+e)/(1-e)^7*M^5/factorial Alternatively, we can expand $E$ in $e$ instead of $M$, giving the solution (the trivial solution when $e = 0$ is $E_0=M$): ```@example perturb E_taylor′ = series(E, e, 0:5) -E_coeffs′ = taylor_coeff(E_taylor′, e) # TODO: get at creation +E_coeffs′ = taylor_coeff(E_taylor′, e) kepler_eqs′ = taylor_coeff(substitute(kepler, E => E_taylor′), e, 0:5) E_coeffs_sol′ = solve_cascade(kepler_eqs′, E_coeffs′, M, e) E_pert′ = substitute(E_taylor′, E_coeffs_sol′) diff --git a/src/taylor.jl b/src/taylor.jl index 7218543be..cbb95dc4e 100644 --- a/src/taylor.jl +++ b/src/taylor.jl @@ -1,27 +1,43 @@ """ - series(y, x, [x0=0,] ns; name = nameof(y)) + series(cs, x, [x0=0,], ns=0:length(cs)-1) -Expand the variable `y` in a power series in the variable `x` around `x0` to orders `ns`. +Return the power series in `x` around `x0` to the powers `ns` with coefficients `cs`. + + series(y, x, [x0=0,] ns) + +Return the power series in `x` around `x0` to the powers `ns` with coefficients automatically created from the variable `y`. Examples ======== ```julia -julia> @variables z ϵ -2-element Vector{Num}: +julia> @variables x y[0:3] z +3-element Vector{Any}: + x + y[0:3] z - ϵ -julia> series(z, ϵ, 2, 0:3) -z[0] + z[1]*(-2 + ϵ) + z[2]*((-2 + ϵ)^2) + z[3]*((-2 + ϵ)^3) +julia> series(y, x, 2) +y[0] + (-2 + x)*y[1] + ((-2 + x)^2)*y[2] + ((-2 + x)^3)*y[3] + +julia> series(z, x, 2, 0:3) +z[0] + (-2 + x)*z[1] + ((-2 + x)^2)*z[2] + ((-2 + x)^3)*z[3] ``` """ -function series(y, x, x0, ns; name = nameof(y)) - c, = @variables $name[ns] - return sum(c[n] * (x - x0)^n for n in ns) +function series(cs::AbstractArray, x::Number, x0::Number, ns::AbstractArray = 0:length(cs)-1) + length(cs) == length(ns) || error("There are different numbers of coefficients and orders") + s = sum(c * (x - x0)^n for (c, n) in zip(cs, ns)) + return s +end +function series(cs::AbstractArray, x::Number, ns::AbstractArray = 0:length(cs)-1) + return series(cs, x, 0, ns) +end +function series(y::Num, x::Number, x0::Number, ns::AbstractArray) + cs, = @variables $(nameof(y))[ns] + return series(cs, x, x0, ns) end -function series(y, x, ns; kwargs...) - return series(y, x, 0, ns; kwargs...) +function series(y::Num, x::Number, ns::AbstractArray) + return series(y, x, 0, ns) end """ diff --git a/test/taylor.jl b/test/taylor.jl index 6f995dda2..770cdddbd 100644 --- a/test/taylor.jl +++ b/test/taylor.jl @@ -1,5 +1,23 @@ using Symbolics +# test all variations of series() input +ns = 0:3 +Y, = @variables y[ns] +@variables x y + +# 1) first argument is a variable +@test isequal(series(y, x, 7, ns), Y[0] + Y[1]*(x-7)^1 + Y[2]*(x-7)^2 + Y[3]*(x-7)^3) +@test isequal(series(y, x, ns), substitute(series(y, x, 7, ns), x => x + 7)) +@test_throws ErrorException series(2*y, x, ns) # 2*y is a meaningless variable name + +# 2) first argument is coefficients +@test isequal(series(Y, x, 7), series(y, x, 7, ns)) +@test isequal(series(Y, x, ns), substitute(series(Y, x, 7), x => x + 7)) +@test isequal(series([1,2,3], 8, 4, [5,6,7]), 1*(8-4)^5 + 2*(8-4)^6 + 3*(8-4)^7) +@test isequal(series([1,2,3], 8, 4), 1 + 2*(8-4)^1 + 3*(8-4)^2) +@test isequal(series([1,2,3], 4, [5,6,7]), series([1,2,3], 8, 4, [5,6,7])) +@test isequal(series([1,2,3], 4), 1^0 + 2*4^1 + 3*4^2) + # 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(x^n//factorial(n) for n in 0:9) == 0