Add pdfs to paper (#202)

* Add Fig1.pdf and Fig2.pdf, and include them in paper.md

* Experiment with the width specification of figures

* Adjustement on the width of figs

* Delete code blocks, and replace them by pdf figures.

* Crop images to adjust them

* Adjust width of figures

* Clarifications to paper

* Fix a minor typo

paper/paper.md

@@ -22,25 +22,26 @@ bibliography: paper.bib
 
 # Summary
 
-The purpose of `TaylorSeries.jl` is to provide a framework to exploit
-Taylor polynomials in one and more variables
-in the [Julia programming language](https://julialang.org) [@julia]. In
-certain cases, it provides a primitive CAS (computer algebra system),
+The purpose of the `TaylorSeries.jl` package is to provide a framework to exploit Taylor polynomials in one and several variables
+in the [Julia programming language](https://julialang.org) [@julia].
+It can be thought of as providing a primitive CAS (computer algebra system),
 which works numerically and not symbolically.
-The package allows to manipulate polynomials of a specified maximum
-degree, including power and composition, as well as series expansions
-of some elementary functions on polynomials, e.g. `exp`,
+The package allows the user to define dense polynomials $p(x)$ of one variable and $p(\mathbf{x})$ of several variables with a specified maximum degree, and perform operations on them, including powers and composition, as well as series expansions
+for elementary functions of polynomials, for example $\exp[p(x)]$,
 where techniques of automatic differentiation are used
 [@Tucker:ValidatedNumerics; @HaroEtAl:ParameterizMeth]. Differentiation and
 integration are also implemented.
 
-Two basic immutable types are defined, `Taylor1` and `TaylorN`,
-which represent the series expansions in one or several variables,
-respectively. These are essentially vectors of coefficients,
-ordered by increasing degree. In the case of `TaylorN`, the
-coefficients are `HomogeneousPolynomials`, which in turn are vectors
+Two basic immutable types are defined, `Taylor1{T}` and `TaylorN{T}`,
+which represent polynomials in one and several variables, respectively; the maximum degree is a field of the types. These types are parametrized by the type `T` of the polynomial coefficients; they essentially consist of one-dimensional arrays of coefficients, ordered by increasing degree.
+
+In the case of `TaylorN`, the
+coefficients are `HomogeneousPolynomial`s, which in turn are vectors
 of coefficients representing all monomials with a given number of variables
-and order (total degree), in some lexicographical order.
+and order (total degree), ordered lexicographically. Higher degree
+polynomials require more
+memory allocation, especially for several variables; while we have not extensively tested the limits of the degree of the polynomials that can be used, `Taylor1` polynomials up to degree 80 and
+`TaylorN` polynomials up to degree 60 in 4 variables have been successfully used. Note that the current implementation of multi-variable series builds up extensive tables in memory which allow to speed up the index calculations.
 
 Julia's parametrized type system allows the construction of Taylor series whose coefficient type is any subtype of the `Number` abstract type. Use cases include complex numbers,
 arbitrary precision `BigFloat`s [@MPFR],
@@ -67,73 +68,23 @@ examples, as well as a detailed user guide, can be found in the
 As a first example we describe how to generate the [Hermite polynomials][@HermitePols_wikipedia]
 ("physicist's" version) up to a given maximum order. Firstly we directly exploit the recurrence relation satisfied by the polynomials.
 
-```julia
-julia> using TaylorSeries
-
-julia> displayBigO(false)
-
-julia> function hermite_polynomials(::Type{T}, nmax::Int) where {T <: Integer}
-
-           x = Taylor1(T, nmax)    # Taylor variable
-           H = fill(x, nmax + 1)   # vector of Taylor series to be overwritten
-
-           H[1] = 1   # order 0
-           H[2] = 2x  # order 1
-
-           for n in 2:nmax
-               # recursion relation for order n:
-               H[n+1] = 2x * H[n] - 2(n-1) * H[n-1]
-           end
-
-           return H
-     end
-
-julia> hermite_polynomials(n) = hermite_polynomials(Int, n);
-
-julia> H = hermite_polynomials(10);
-
-julia> function hermite_polynomial(n::Int)
-    @assert 0 ≤ n ≤ length(H) "Not enough Hermite polynomials generated"
-    return H[n+1]
-end
-
-julia> hermite_polynomial(6)  # degree 6
-- 120 + 720 t² - 480 t⁴ + 64 t⁶
-
-```
-
+![Code to generate Hermite polynomials directly from the recursion relation; the last line displays the 6th Hermite polynomial.](Fig1.pdf){width=110%}  
 The example above can be slightly modified to compute, for example, the 100th Hermite polynomial.
 In this case, the coefficients will be larger than $2^{63}-1$, so the modular
-behavior under overflow of the standard `Int64` type will not suffice. Rather, the polynomials should
+behavior, under overflow of the standard `Int64` type, will not suffice. Rather, the polynomials should
 be generated with `hermite_polynomials(BigInt, 100)` to ensure
-the use of arbitrary length integers.
+the use of arbitrary-length integers.
 
 ## Using a generating function
 As a second example, we describe a numerical way of obtaining the
-Hermite polynomials from their generating function: the $n$-th Hermite polynomial
-corresponds to the $n$-th derivative of the function $\exp(2t \, x - t^2)$.
-
-```julia
-julia> 𝒢(x,t) = exp(2t * x - t^2)  # generating function; 𝒢 is typed as \scrG<TAB>
-
-julia> xn = set_variables("x", numvars=1, order=10)
-
-julia> x = xn[1]
-
-julia> t = Taylor1([zero(x), one(x)], 10)  # Taylor1{TaylorN{Float64}}
-
-julia> gf = 𝒢(x, t)  # Taylor1 expansion of 𝒢
-
-julia> HH(n::Int) = derivative(n, gf)  # n-th derivative of `gf`
-
-julia> HH(6)
-- 120.0 + 720.0 x₁² - 480.0 x₁⁴ + 63.99999999999999 x₁⁶
-```
+Hermite polynomials from their generating function: the $n$th Hermite polynomial
+corresponds to the $n$th derivative of the function $\exp(2t \, x - t^2)$.
 
+![Code to generate Hermite polynomials from the generating function $\exp(2t \, x - t^2)$; the last line displays the result for the 6th Hermite polynomial.](Fig2.pdf){width=110%}  
 This example shows that the calculations are performed numerically and not
 symbolically, using `TaylorSeries.jl` as a polynomial manipulator; this
- is manifested by the fact that the last coefficient of `HH(6)` is not
- identical to an integer.
+is manifested by the fact that the last coefficient of `HH(6)` is not
+identical to an integer.
 
 ## Taylor method for integrating ordinary differential equations
 As a final example, we give a simple implementation of Picard
@@ -144,41 +95,11 @@ We consider the initial-value problem $\dot{x} = x$,
 with initial condition $x(0) = 1$. One step of the integration corresponds
 to constructing the Taylor series of the solution $x(t)$ in powers of $t$:
 
- ```julia
-
-julia> ∫⬩dt(u::Taylor1) = integrate(u) # the symbol ∫ is obtained as \int<TAB>
-
-julia> function taylor_step(f, u0)
-
-           u = copy(u0)
-           unew = u0 + ∫⬩dt(f(u))
-
-           while unew != u
-               u = unew
-               unew = u0 + ∫⬩dt(f(u))   # Picard iteration
-           end
-
-           return u
-       end
-
-julia> f(x) = x  # Differential equation
-
-julia> order = 20  # maximum order of the Taylor expansion for the solution
-
-julia> u0 = Taylor1([1.0], order)  # initial condition given as a Taylor expansion
-
-julia> solution = taylor_step(f, u0);   # solution
-
-julia> solution(1.0) - exp(1.0) # compare the solution evaluated at t=1 with the exact value
-0.0
-
-```
-
+![Code to implement Picard iteration to integrate the initial value problem $\dot{x} = x$, $x(0) = 1$, using a 20th order local Taylor expansion.](Fig3.pdf){width=110%}  
 Thus this Taylor expansion of order 20 around $t_0=0$
 suffices to obtain the exact solution at $t=1$, while the error at time $t=2$
 from the same expansion is $4.53 \times 10^{-14}$.
-This indicates that a proper treatment should try to estimate what size step
-should be taken as a function of the behavior of the solution.
+This indicates that a proper treatment should estimate the size of the required step that should be taken as a function of the solution.
 
 ## Acknowledgements