diff --git a/docs/make.jl b/docs/make.jl index c99454081..855ff3421 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -32,15 +32,14 @@ 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", "manual/expression_manipulation.md", "manual/derivatives.md", + "manual/taylor.md", "manual/groebner.md", "manual/solver.md", "manual/arrays.md", diff --git a/docs/src/examples/perturbation.md b/docs/src/examples/perturbation.md deleted file mode 100644 index c57ff8bcb..000000000 --- a/docs/src/examples/perturbation.md +++ /dev/null @@ -1,250 +0,0 @@ -# [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 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 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. - -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/) - -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: - -```@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) - - for i = 1:maxiter - xₙ₊₁ = substitute(fₙ₊₁, Dict(x => xₙ)) - if abs(xₙ₊₁ - xₙ) < abstol - return xₙ₊₁ - else - xₙ = xₙ₊₁ - end - end - return xₙ₊₁ -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 -``` - -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: - -```@example perturb -n = 2 -@variables ϵ a[1:n] -``` - -Then, we define - -```@example perturb -x = 1 + a[1]*ϵ + a[2]*ϵ^2 -``` - -The next step is to substitute `x` in the problem equation - -```@example perturb - eq = x^5 + ϵ*x - 1 -``` - -The expanded form of `eq` is - -```@example perturb -expand(eq) -``` - -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. - -```@example perturb -function collect_powers(eq, x, ns) - eq = expand(eq) - [Symbolics.coeff(eq, x^i) for i in ns] -end -``` - -To return the coefficients of $ϵ$ and $ϵ^2$ in `eq`, we can write - -```@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() - - for i = 1:length(ps) - eq = substitute(eqs[i], vals) - vals[ps[i]] = symbolic_linear_solve(eq, ps[i]) - 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: - -```@example perturb -vals = solve_coef(eqs, a) -``` - -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. - -```@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 -``` - -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. - -```@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) -``` - -We start by defining the variables (assuming `n = 3`): - -```@example perturb -n = 3 -@variables ϵ M a[1:n] -x = def_taylor(ϵ, a, M) -``` - -We further simplify by substituting `sin` with its power series using the `expand_sin` helper function: - -```@example perturb -expand_sin(x, n) = sum([(isodd(k) ? -1 : 1)*(-x)^(2k-1)/factorial(2k-1) for k=1:n]) -``` - -To test, - -```@example perturb -expand_sin(0.1, 10) ≈ sin(0.1) -``` - -The problem equation is - -```@example perturb -eq = x - ϵ * expand_sin(x, n) - M -``` - -We follow the same process as the first example. We collect the coefficients of the powers of `ϵ` - -```@example perturb -eqs = collect_powers(eq, ϵ, 1:n) -``` - -and then solve for `a`: - -```@example perturb -vals = solve_coef(eqs, a) -``` - -Finally, we substitute `vals` back in `x`: - -```@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 - -``` -(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$. diff --git a/docs/src/manual/taylor.md b/docs/src/manual/taylor.md new file mode 100644 index 000000000..0139e5d43 --- /dev/null +++ b/docs/src/manual/taylor.md @@ -0,0 +1,9 @@ +# Taylor Series + +For a real example of how to use the Taylor series functionality, see [this tutorial](@ref perturb_alg). + +```@docs +series +taylor +taylor_coeff +``` diff --git a/docs/src/tutorials/perturbation.md b/docs/src/tutorials/perturbation.md new file mode 100644 index 000000000..30963997e --- /dev/null +++ b/docs/src/tutorials/perturbation.md @@ -0,0 +1,154 @@ +# [Mixed Symbolic-Numeric Perturbation Theory](@id perturb_alg) + +[**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. + +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. + +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. + +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/). + +## 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 +@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) + + 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ₙ₊₁ # converged + else + xₙ = xₙ₊₁ + end + end + error("Newton's method failed to converge") +end + +x_newton = solve_newton(quintic, x, 1.0) +println("Newton's method solution: x = ", x_newton) +``` + +To solve the same problem with perturbation theory, we must introduce an expansion variable $ϵ$ in the equation: +```@example perturb +@variables ϵ # expansion variable +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_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 +quintic_taylor = substitute(quintic, x => x_taylor) +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, ϵ, 0:7) +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₀) # store solutions in a map + + # verify that x₀ is a solution of the first equation + eq0 = substitute(eqs[1], sol) + isequal(eq0.lhs, eq0.rhs) || error("$sol does not solve $(eqs[1])") + + # solve remaining equations sequentially + for i in 2:length(eqs) + eq = substitute(eqs[i], sol) # insert previous solutions + 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 +end +``` +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 + 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! + +## 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 +@variables e E M +kepler = E - e * sin(E) ~ M +``` +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 +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) +``` + +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) # 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 +``` +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, 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 $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) +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′) +``` +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 +@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) +``` diff --git a/src/Symbolics.jl b/src/Symbolics.jl index cd84fd16a..146841893 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 series, taylor, taylor_coeff + 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..cbb95dc4e --- /dev/null +++ b/src/taylor.jl @@ -0,0 +1,142 @@ +""" + series(cs, x, [x0=0,], ns=0:length(cs)-1) + +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 x y[0:3] z +3-element Vector{Any}: + x + y[0:3] + z + +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(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::Num, x::Number, ns::AbstractArray) + return series(y, x, 0, ns) +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 + # return array of expressions/equations for each order + return taylor_coeff.(Ref(f), Ref(x), n; rationalize) + elseif f isa Equation + 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" + 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 = Base.rationalize(c) # convert integers/floats to rational numbers; avoid name clash between rationalize and Base.rationalize() + end + return c +end + +""" + taylor(f, x, [x0=0,] n; rationalize=true) + +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 +======== +```julia +julia> @variables x +1-element Vector{Num}: + 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) + +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...) + 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 + + 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 + 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/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..770cdddbd --- /dev/null +++ b/test/taylor.jl @@ -0,0 +1,66 @@ +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 +@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*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 +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) - (π/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 +@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 + +# 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) + + +# 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) + +# system of equations +@test taylor(exp(im*x) ~ 0, x, 0:5) == taylor([cos(x) ~ 0, sin(x) ~ 0], x, 0:5)