Merge pull request #285 from mkg33/master

Documentation changes

Author: ChrisRackauckas

README.md

@@ -6,74 +6,74 @@
 [![Coverage Status](https://coveralls.io/repos/JuliaDiffEq/DiffEqOperators.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/JuliaDiffEq/DiffEqOperators.jl?branch=master)
 [![codecov.io](http://codecov.io/github/shivin9/DiffEqOperators.jl/coverage.svg?branch=master)](http://codecov.io/github/JuliaDiffEq/DiffEqOperators.jl?branch=master)
 
-Construct high order finite-difference operators to discretize a partial
+Construct high-order finite-difference operators to discretize a partial
 differential equation and its boundary conditions by the method of
-lines.  This reduces it to a system of ordinary differential
+lines. This reduces it to a system of ordinary differential
 equations that can be solved by [`DifferentialEquations.jl`](https://juliadiffeq.org/).
 
 Both centered and [upwind](https://en.wikipedia.org/wiki/Upwind_scheme)
 operators are provided, for domains of any dimension, arbitrarily
-spaced grids, and for any order of accuracy.  The cases of 1, 2 and
+spaced grids, and for any order of accuracy. The cases of 1, 2, and
 3 dimensions with an evenly spaced grid are optimized with a
-convolution routine from `NNlib.jl`.  Care is taken to avoid
+convolution routine from `NNlib.jl`. Care is taken to avoid
 unnecessary allocations.
 
-Any operator can be concretised as an `Array`, a `BandedMatrix` or
+Any operator can be concretized as an `Array`, a `BandedMatrix` or
 a sparse matrix.
 
 ## The simplest case
 
 As shown in the figure, the operators act on a set of samples
-`f_j = f(x_j)` for a function `f` at a grid of points `x_j`.  The
+`f_j = f(x_j)` for a function `f` at a grid of points `x_j`. The
 grid has `n` interior points at `x_j = jh` for `j = 1` to `n`, and 2
-boundary points at `x_0 = 0` and `x_{n+1} = (n+1) h`.  The input to
+boundary points at `x_0 = 0` and `x_{n+1} = (n+1) h`. The input to
 the numerical operators is a vector `u = [f_1, f_2, …, f_N]`, and
 they output a vector of sampled derivatives `du ≈ [f'(x_1), f'(x_2),
-…, f'(x_N)]`, or a higher order derivative as requested.
+…, f'(x_N)]`, or a higher-order derivative as requested.
 
 A numerical derivative operator `D` of order `m` can be constructed
-for this grid with `D = CenteredDifference(1, m, h, n)`  The argument
-`1` indicates that this is a first derivative.  Order `m` means
+for this grid with `D = CenteredDifference(1, m, h, n)`. The argument
+`1` indicates that this is the first derivative. Order `m` means
 that the operator is exact up to rounding when `f` is a polynomial
 of degree `m` or lower.
 
 The derivative operator `D` is used along with a boundary condition
-operator `Q`, to compute derivatives at the interior points of the
-grid.  A simple boundary condition `f(x_0) = f(x_n+1) = 0` is
+operator `Q` to compute derivatives at the interior points of the
+grid. A simple boundary condition `f(x_0) = f(x_n+1) = 0` is
 constructed with `Q = Dirichlet0BC(eltype(u))`.
 
 Given these definitions, the derivatives are calculated as if the
-operators `D` and `Q` were matrices.  `du = D*Q*u`.  This is an
-abuse of notation!  The particular `Q` in this example is a linear
-operator, but in general boundary conditions are affine operators.
+operators `D` and `Q` were matrices, `du = D*Q*u`. This is an
+abuse of notation! The particular `Q` in this example is a linear
+operator but, in general, boundary conditions are affine operators.
 They have the form `Q(x) = M*x + c`, where `M` is a matrix and `c`
-is a constant vector.  As a consequence, `Q` can not be concretized
+is a constant vector. As a consequence, `Q` cannot be concretized
 to a matrix.
 
 ![Actions of DiffEqOperators on interior points and ghost points](action.svg)
 
 The operator `D` works by interpolating a polynomial of degree `m`
-through `m+1` adjacent points on the grid.  Near the middle of the
+through `m+1` adjacent points on the grid. Near the middle of the
 grid, the derivative is approximated at `x_j` by interpolating a
-polynomial of order `m` with `x_j` at its centre.  To define a
-order-`m` polynomial, values are required at `m+1` points.  When
+polynomial of order `m` with `x_j` at its centre. To define an
+order-`m` polynomial, values are required at `m+1` points. When
 `x_j` is too close to the boundary for that to fit, the polynomial
 is interpolated through the leftmost or rightmost `m+1` points,
 including two “ghost” points that `Q` appends on the boundaries.
 The numerical derivatives are linear combinations of the values
-through which the polynomials are interpolated.  The vectors of
+through which the polynomials are interpolated. The vectors of the
 coefficients in these linear combinations are called “stencils”.
 Because `D` takes values at the ghost points and returns values at
-the interior points, it is a `n×(n+2)` matrix.
+the interior points, it is an `n×(n+2)` matrix.
 
 The boundary condition operator `Q`
-acts as an `(n+2)×n` matrix.  The output `Q*u` is a vector of values
+acts as an `(n+2)×n` matrix. The output `Q*u` is a vector of values
 on the `n` interior and the 2 boundary points, `[a, f(x_1), …, f(x_N), b]`.
-The interior points take the values of `u`.  The values `a` and `b` are
-samples at “ghost” points on the grid boundaries.  As shown, these
+The interior points take the values of `u`. The values `a` and `b` are
+samples at “ghost” points on the grid boundaries. As shown, these
 values are assigned so that an interpolated polynomial `P(x)` satisfies
-the left hand boundary condition, and `Q(x)` satisfies the right hand
-boundary condition.  The boundary conditions provided by the
+the left hand boundary condition, and `Q(x)` satisfies the right-hand
+boundary condition. The boundary conditions provided by the
 library are precisely those for which the values `a` and `b` are affine
 functions of the interior values `f_j`, so that `Q` is an affine operator.
 
@@ -84,28 +84,28 @@ and the derivative and boundary condition operators are similar
 to matrices.
 
 In two dimensions, the values `f(x_j)` are naturally stored as a
-matrix.  Taking derivatives along the downwards axis is easy, because
-matrices act columnwise.  Horizontal derivatives can be taken by
-transposing the matrices.  The derivative along the rightward axis
-is `(D*F')' = F*D' `.  This is easy to code, but less easy to read
+matrix. Taking derivatives along the downwards axis is easy, because
+matrices act columnwise. Horizontal derivatives can be taken by
+transposing the matrices. The derivative along the rightward axis
+is `(D*F')' = F*D' `. This is easy to code, but less easy to read
 for those who haven't seen it before.
 
 When a function has three or more arguments, its values are naturally
-stored in a higher dimensional array.  Julia's multiplication
+stored in a higher-dimensional array. Julia's multiplication
 operator is only defined for `Vector` and `Matrix`, so applying an
 operator matrix to these arrays would require a complicated and
 error prone series of `reshape` and axis permutation functions.
 
 Therefore the types of derivative and boundary condition operators
-are parameterised by the axis along which the operator acts.  With
+are parameterised by the axis along which the operator acts. With
 derivative operators, the axis is supplied as a type parameter.
 The simple case `CenteredDifference(…)` is equivalent to
 `CenteredDifference{1}(…)`, row-wise derivatives are taken by
 `CenteredDifference{2}(…)`, sheet-wise by `CenteredDifference{3}(…)`,
 and along the `N`th axis by `CenteredDifference{N}(…)`.
 
-Boundary conditions are more complicated.  See `@doc MultiDimBC`
-for how they are supposed to work in multiple dimensions.  They
+Boundary conditions are more complicated. See `@doc MultiDimBC`
+for how they are supposed to work in multiple dimensions. They
 don't currently work that way.
 
 ## Constructors
@@ -125,15 +125,15 @@ UpwindDifference{N}(derivative_order::Int,
 The arguments are:
 
 - `N`: The directional dimension of the discretization. If `N` is not given,
-  it is assumed to be 1, i.e. differencing occurs along columns.
+  it is assumed to be 1, i.e., differencing occurs along columns.
 - `derivative_order`: the order of the derivative to discretize.
 - `approximation_order`: the order of the discretization in terms of O(dx^order).
 - `dx`: the spacing of the discretization. If `dx` is a `Number`, the operator
   is a uniform discretization. If `dx` is an array, then the operator is a
   non-uniform discretization.
 - `len`: the length of the discretization in the direction of the operator.
 - `coeff_func`: An operational argument for a coefficient function `f(du,u,p,t)`
-  which sets the coefficients of the operator. If `coeff_func` is a `Number`
+  which sets the coefficients of the operator. If `coeff_func` is a `Number`,
   then the coefficients are set to be constant with that number. If `coeff_func`
   is an `AbstractArray` with length matching `len`, then the coefficients are
   constant but spatially dependent.
@@ -159,11 +159,11 @@ The following concretizations are provided:
 
 Additionally, the function `sparse` is overloaded to give the most efficient
 matrix type for a given operator. For one-dimensional derivatives this is a
-`BandedMatrix`, while for higher dimensional operators this is a `BlockBandedMatrix`.
+`BandedMatrix`, while for higher-dimensional operators this is a `BlockBandedMatrix`.
 The concretizations are made to act on `vec(u)`.
 
-A contraction operator concretizes to an ordinary matrix, no matter which dimension 
-the contraction acts along, by doing the kroncker product formulation. I.e., the
+A contraction operator concretizes to an ordinary matrix, no matter which dimension
+the contraction acts along, by doing the Kronecker product formulation. I.e., the
 action of the built matrix will match the action on `vec(u)`.
 
 ## Boundary Condition Operators
@@ -281,4 +281,4 @@ will inherit the appropriate action.
 ### Efficiency of Composed Operator Actions
 
 Composed operator actions utilize NNLib.jl in order to do cache-efficient
-convolution operations in higher dimensional combinations.
+convolution operations in higher-dimensional combinations.

examples/heat_equation.jl

@@ -1,6 +1,6 @@
 # # Heat Equation
 #
-# This example demonstrates how to combine `OrdinaryDiffEq` with `DiffEqOperators` to solve a time dependent PDE.
+# This example demonstrates how to combine `OrdinaryDiffEq` with `DiffEqOperators` to solve a time-dependent PDE.
 # We consider the heat equation on the unit interval, with Dirichlet boundary conditions:
 # ∂ₜu = Δu
 # u(x=0,t)  = a
@@ -35,10 +35,9 @@ sol = solve(prob, alg)
 using Test
 @test u_analytic.(knots, t1) ≈ sol[end] rtol=1e-3
 
-# Because the creation of boundary conditions is cheap, we can 
+# Because the creation of boundary conditions is cheap, we can
 # implement time-dependent boundary conditions as follows:
 # function step(u,p,t)
 #    bc = DirichletBC(cos(t*1.5), sin(t))
 #    Δ*bc*u
 # end
-

examples/poisson.jl

@@ -1,8 +1,8 @@
 # # Poisson Equation
 #
 # We want to solve the [poisson equation](https://en.wikipedia.org/wiki/Poisson_equation) on the unit interval. It is given by
-# `Δu = f` with boundary conditions `u(0) = a` and `u(1) = b`
-# First of all let us choose some values for the parameters and remark, that there is an exact solution:
+# `Δu = f` with boundary conditions `u(0) = a` and `u(1) = b`.
+# First of all, let us choose some values for the parameters and remark that there is an exact solution:
 f = 1.0
 a = -1.0
 b = 2.0
@@ -20,16 +20,16 @@ ord_approx = 2
 Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)
 bc = DirichletBC(a, b)
 
-# Before solving the equation, lets take a look at Δ and bc:
+# Before solving the equation, let's take a look at Δ and bc:
 # display(Array(Δ))
 # display(bc*zeros(nknots))
-# We see that `Δ` is a (lazy) matrix with the laplace stencil extended over the boundaries.
+# We see that `Δ` is a (lazy) matrix with the Laplace stencil extended over the boundaries.
 # And `bc` acts by padding the values just outside the boundaries.
 
 u = (Δ*bc) \ fill(f, nknots)
 knots = range(h, step=h, length=nknots)
 
-# Since we used a second order approximation and the analytic solution itself was a second order
-# polynomial, we expect that they are equal up to rounding errors:
+# Since we used a second order approximation and the analytic solution itself was a second-order
+# polynomial, we expect them to be equal up to rounding errors:
 using Test
 @test u ≈ u_analytic.(knots)

src/MOL_discretization.jl

@@ -17,16 +17,16 @@ function get_bcs(bcs,tdomain,domain)
                 lhs_deriv_depvars_bcs[var] = Array{Expr}(undef,3)
             end
             j = 0
-            if isequal(bcs[i].lhs.args[1],tdomain.lower) # u(t=0,x) 
-                j = 1    
+            if isequal(bcs[i].lhs.args[1],tdomain.lower) # u(t=0,x)
+                j = 1
             elseif isequal(bcs[i].lhs.args[2],domain.lower) # u(t,x=x_init)
                 j = 2
             elseif isequal(bcs[i].lhs.args[2],domain.upper) # u(t,x=x_final)
                 j = 3
             end
             if bcs[i].rhs isa ModelingToolkit.Constant
                 lhs_deriv_depvars_bcs[var][j] = :(var=$(bcs[i].rhs.value))
-            else    
+            else
                 lhs_deriv_depvars_bcs[var][j] = Expr(bcs[i].rhs)
             end
         end
@@ -38,14 +38,14 @@ end
 # Recursively traverses the input expression (rhs), replacing derivatives by
 # finite difference schemes. It returns a time dependent expression (expr)
 # that will be evaluated in the "f" ODE function (in DiffEqBase.discretize),
-# Note: 'non-derived' dependent variables are inserted into the diff. equations
-#       E.g. Dx(u(t,x))=v(t,x)*Dx(u(t,x)), v(t,x)=t*x
+# Note: 'non-derived' dependent variables are inserted into the differential equations
+#       E.g., Dx(u(t,x))=v(t,x)*Dx(u(t,x)), v(t,x)=t*x
 #            =>  Dx(u(t,x))=t*x*Dx(u(t,x))
-                            
+
 function discretize_2(input,deriv_order,approx_order,dx,X,len,
                       deriv_var,dep_var_idx,indep_var_idx)
     if input isa ModelingToolkit.Constant
-        return :($(input.value)) 
+        return :($(input.value))
     elseif input isa Operation
         if input.op isa Variable
             expr = :(0.0)
@@ -64,7 +64,7 @@ function discretize_2(input,deriv_order,approx_order,dx,X,len,
                 if deriv_order == 0
                     expr = :(u[:,$j])
                 elseif deriv_order == 1
-                    # TODO: approx_order and forward/backward should be 
+                    # TODO: approx_order and forward/backward should be
                     #       input parameters of each derivative
                     approx_order = 1
                     L = UpwindDifference(deriv_order,approx_order,dx[i],len[i]-2,-1)
@@ -92,7 +92,7 @@ function discretize_2(input,deriv_order,approx_order,dx,X,len,
                 aux_2 = discretize_2(input.args[2],deriv_order,approx_order,dx,X,
                                      len,deriv_var,dep_var_idx,indep_var_idx)
                 return :(broadcast($(input.op), $aux_1, $aux_2))
-            end    
+            end
         end
     end
 end
@@ -105,22 +105,22 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     #   1) PDE System
     #        1.a) Transient
     #                There is more than one indep. variable, including  't'
-    #                E.g. du/dt = d2u/dx2 + d2u/dy2 + f(t,x,y)
+    #                E.g., du/dt = d2u/dx2 + d2u/dy2 + f(t,x,y)
     #        1.b) Stationary
     #                There is more than one indep. variable, 't' is not included
-    #                E.g. 0 = d2u/dx2 + d2u/dy2 + f(x,y)
+    #                E.g., 0 = d2u/dx2 + d2u/dy2 + f(x,y)
     #   2) ODE System
     #        't' is the only independent variable
     #        The ODESystem is packed inside a PDESystem
-    #        E.g. du/dt = f(t)
+    #        E.g., du/dt = f(t)
     #
     #   Note: regarding input format, lhs must be "du/dt" or "0".
     #
 
-    # The following code deals with 1.a case for 1D,
-    # i.e. only considering 't' and 'x'
+    # The following code deals with 1.a case for 1-D,
+    # i.e., only considering 't' and 'x'.
+
 
-    
     ### Declare and define independent-variable data structures ###############
 
     tdomain = 0.0
@@ -155,9 +155,9 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     approx_order = discretization.order
 
     lhs_deriv_depvars = Dict()
-    dep_var_idx = Dict() 
+    dep_var_idx = Dict()
     dep_var_disc = Dict() # expressions evaluated in the ODE function (f)
-    
+
     # if there is only one equation
     if pdesys.eq isa Equation
         eqs = [pdesys.eq]
@@ -181,7 +181,7 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
         # TODO: is there a better way to convert an Expr into a Function?
         dep_var_disc[var] = @eval (Q,u,t) -> $aux
     end
-    
+
     ### Declare and define boundary conditions ################################
 
     # TODO: extend to Neumann BCs and Robin BCs
@@ -191,7 +191,7 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     u_x0 = Array{Any}(undef,no_dep_vars)
     u_x1 = Array{Any}(undef,no_dep_vars)
     Q = Array{RobinBC}(undef,no_dep_vars)
-    
+
     for var in keys(dep_var_idx)
         j = dep_var_idx[var]
         bcs = lhs_deriv_depvars_bcs[var]
@@ -210,7 +210,7 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     ### Define the discretized PDE as an ODE function #########################
 
     function f(du,u,p,t)
-        
+
         # Boundary conditions can vary with respect to time
         for j in 1:no_dep_vars
             a = Base.invokelatest(u_x0[j],X[2][1],t)
@@ -233,4 +233,3 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     # Return problem ##########################################################
     return PDEProblem(ODEProblem(f,u_t0,(tdomain.lower,tdomain.upper),nothing),Q,X)
 end
-

src/boundary_padded_arrays.jl

@@ -33,7 +33,7 @@ function Base.getindex(Q::BoundaryPaddedVector,i)
 end
 
 """
-Higher dimensional generalization of BoundaryPaddedVector, pads an array of dimension N along the dimension D with 2 Arrays of dimension N-1, stored in lower and upper
+Higher-dimensional generalization of BoundaryPaddedVector pads an array of dimension N along the dimension D with 2 Arrays of dimension N-1, stored in lower and upper
 
 """
 struct BoundaryPaddedArray{T, D, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractDirectionalBoundaryPaddedArray{T,N, D}
@@ -60,9 +60,9 @@ Example:
 A = compose(Ax, Ay, Az) # 3D domain
 A = compose(Ax, Ay) # 2D Domain
 
-Composes BoundaryPaddedArrays that extend the same u for each different dimension that u has in to a ComposedBoundaryPaddedArray
+Composes BoundaryPaddedArrays that extend the same u for each different dimension that u has in to a ComposedBoundaryPaddedArray.
 
-Ax Ay and Az can be passed in any order, as long as there is exactly one BoundaryPaddedArray that extends each dimension.
+Ax, Ay, and Az can be passed in any order, as long as there is exactly one BoundaryPaddedArray that extends each dimension.
 """
 function compose(padded_arrays::BoundaryPaddedArray...)
     N = ndims(padded_arrays[1])
@@ -105,7 +105,7 @@ Ax, Ay,... = decompose(A::ComposedBoundaryPaddedArray)
 
 -------------------------------------------------------------------------------------
 
-Decomposes a ComposedBoundaryPaddedArray in to components that extend along each dimension individually
+Decomposes a ComposedBoundaryPaddedArray in to components that extend along each dimension individually.
 """
 decompose(A::ComposedBoundaryPaddedArray) = Tuple([BoundaryPaddedArray{gettype(A), i, ndims(A), ndims(A)-1, typeof(lower[1])}(A.lower[i], A.upper[i], A.u) for i in 1:ndims(A)])