Let user pass coefficients to series(); show both this and the version with automatic coefficients in tutorial

Author: hersle

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′)

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
 
 """

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