use Literate in usage examples (#818)

Author: jishnub

docs/make.jl

@@ -20,7 +20,8 @@ for (example, included) in [
             ("Sampling.jl", String[]),
             ("Periodic.jl", ["Periodic1.jl"]),
             ("Eigenvalue.jl", String[]),
-            ("NonlinearBVP.jl", ["NonlinearBVP1.jl"])
+            ("NonlinearBVP.jl", ["NonlinearBVP1.jl"]),
+            ("system_of_eqn.jl", String[])
             ]
     filename = joinpath(example_dir, example)
     Literate.markdown(filename, output_dir, documenter=true,
@@ -44,6 +45,7 @@ makedocs(
                     ],
                     "Examples" => [
                         "generated/ODE.md",
+                        "generated/system_of_eqn.md",
                         "generated/PDE.md",
                         "generated/Sampling.md",
                         "generated/Periodic.md",

docs/src/assets/Harmonic_eigs.pdf

@@ -89,69 +89,6 @@ u = [Dirichlet();
 u(0.1)
 ```
 
-## Eigenvalue Problems
-
-In analogy to linear algebra, many differential equations may be posed as eigenvalue problems. That is, for some differential operator ``\mathop{L}``, there are a family of functions ``\mathop{u}_i(x)`` such that
-
-```math
-\mathop{L} \mathop{u}_i(x) = λ_i \mathop{u}_i(x),
-```
-
-where ``λ_i`` is the ``i^{th}`` eigenvalue of the ``L`` and has a corresponding *eigenfunction* ``\mathop{u}_i(x)``. A classic eigenvalue problem is known as the quantum harmonic oscillator where
-
-```math
-\mathop{L} = -\frac{1}{2}\frac{\mathop{d}^2}{\mathop{dx}^2} + \frac{1}{2} x^2,
-```
-
-and one demands that ``\mathop{u}(∞) = \mathop{u}(-∞) = 0``. Because we expect the solutions to be exponentially suppressed for large ``x``, we can approximate this with Dirichlet boundary conditions at a 'reasonably large' ``x`` without much difference.
-
-We can express this in ApproxFun as the following:
-
-```julia
-x = Fun(-8..8)
-L = -𝒟^2/2 + x^2/2
-S = space(x)
-B = Dirichlet(S)
-λ, v = ApproxFun.eigs(B, L, 500,tolerance=1E-10)
-```
-
-Note that boundary conditions must be specified in the call to `eigs`. Plotting the first ``20`` eigenfunctions offset vertically by their eigenvalue, we see
-
-![harmonic_eigs](../assets/Harmonic_eigs.pdf)
-
-If the solutions are not relatively constant near the boundary then one should push the boundaries further out.
-
-For problems with different contraints or boundary conditions, `B` can be any zero functional constraint, e.g., `DefiniteIntegral()`.
-
-## Systems of equations
-
-Systems of equations can be handled by creating a matrix of operators and functionals.  For example, we can solve the system
-
-```math
-\begin{gathered}
-    \mathop{u}'' - \mathop{u} + 2 \mathop{v} = \mathop{e}^x,  \\
-    \mathop{u} + \mathop{v}' + \mathop{v} = \cos{x}, \\
-    \mathop{u}(-1) = \mathop{u}'(-1) = \mathop{v}(-1) = 0
-\end{gathered}
-```
-
-using the following code:
-
-```@repl using-pkgs
-x = Fun(); B = Evaluation(Chebyshev(),-1);
-A = [B      0;
-     B*𝒟    0;
-     0      B;
-     𝒟^2-I  2I;
-     I      𝒟+I];
-u,v = A \ [0;0;0;exp(x);cos(x)];
-u(-1),u'(-1),v(-1)
-norm(u''-u+2v-exp(x))
-norm(u+v'+v-cos(x))
-```
-
-In this example, the automatic space detection failed and so we needed to specify explicitly that the domain space for `B` is `Chebyshev()`.
-
 ## QR Factorization
 
 Behind the scenes, `A\b` where `A` is an `Operator` is implemented via an adaptive QR factorization.  That is, it is equivalent to `qr(A)\b`.  (There is a subtlety here in space inferring: `A\b` can use both `A` and `b` to determine the domain space, while `qr(A)` only sees the operator `A`.)  Note that `qr` adaptively caches a partial QR Factorization
@@ -184,24 +121,3 @@ Many PDEs have weak singularities at the corners, in which case it is beneficial
 ```julia
 \(A, [zeros(∂(d));f]; tolerance=1E-6)
 ```
-
-## Nonlinear equations
-
-There is preliminary support for nonlinear equations, via Newton iteration in function space.  Here is a simple two-point boundary value problem:
-
-```math
-\begin{gathered}
-    ϵ \mathop{u}'' + 6(1-x^2) \mathop{u}' + \mathop{u}^2=1, \\
-    \mathop{u}(-1) = \mathop{u}(1) = 0.
-\end{gathered}
-```
-
-This can be solved as follows:
-
-```@repl using-pkgs
-ϵ = 1/70; x = Fun();
-N = u -> [u(-1); u(1); ϵ*u''+6*(1-x^2)*u'+u^2-1];
-u0 = x; u = newton(N,u0);
-u(-1), u(1)
-norm(ϵ*u''+6*(1-x^2)*u'+u^2-1)
-```

docs/src/usage/equations.md

@@ -1,4 +1,52 @@
-# # Self-adjoint Eigenvalue Problem
+# # Eigenvalue problem
+
+# ## Standard eigenvalue problem
+
+# In analogy to linear algebra, many differential equations may be posed as eigenvalue problems.
+# That is, for some differential operator ``\mathop{L}``, there are a family of functions
+# ``\mathop{u}_i(x)`` such that
+# ```math
+# \mathop{L} \mathop{u}_i(x) = λ_i \mathop{u}_i(x),
+# ```
+# where ``λ_i`` is the ``i^{th}`` eigenvalue of the ``L`` and has a corresponding
+# *eigenfunction* ``\mathop{u}_i(x)``.
+# A classic eigenvalue problem is known as the quantum harmonic oscillator where
+# ```math
+# \mathop{L} = -\frac{1}{2}\frac{\mathop{d}^2}{\mathop{dx}^2} + \frac{1}{2} x^2,
+# ```
+# and one demands that ``\mathop{u}(∞) = \mathop{u}(-∞) = 0``.
+# Because we expect the solutions to be exponentially suppressed for large ``x``,
+# we can approximate this with Dirichlet boundary conditions at a 'reasonably large' ``x``
+# without much difference.
+
+# We can express this in ApproxFun as the following:
+using ApproxFun
+using LinearAlgebra
+
+x = Fun(-8..8)
+V = x^2/2
+L = -𝒟^2/2 + V
+S = space(x)
+B = Dirichlet(S)
+λ, v = ApproxFun.eigs(B, L, 500,tolerance=1E-10);
+
+# Note that boundary conditions must be specified in the call to `eigs`.
+# Plotting the first ``20`` eigenfunctions offset vertically by their eigenvalue, we see
+
+import Plots
+p = Plots.plot(V; legend=false, ylim=(-Inf, λ[22]))
+for k=1:20
+    Plots.plot!(real(v[k]/norm(v[k]) + λ[k]))
+end
+p
+
+# If the solutions are not relatively constant near the boundary then one should push
+# the boundaries further out.
+
+# For problems with different contraints or boundary conditions,
+# `B` can be any zero functional constraint, e.g., `DefiniteIntegral()`.
+
+# ## Self-adjoint Eigenvalue Problem
 # Ref:
 # [J. L. Aurentz & R. M. Slevinsky (2019), arXiv:1903.08538](https://arxiv.org/abs/1903.08538)
 
@@ -8,13 +56,10 @@
 # ```
 # where ``u(\pm 10) = 0``, ``V(x) = ωx^2 + x^4``, and ``ω = 25``.
 
-using ApproxFun
-using LinearAlgebra
-
 # Define parameters
 ω = 25.0
 d = -10..10;
-S = Legendre(d)
+S = Legendre(d) # Equivalently, Ultraspherical(0.5, d)
 NS = NormalizedPolynomialSpace(S) # NormalizedLegendre
 D1 = Conversion(S, NS)
 D2 = Conversion(NS, S);

examples/Eigenvalue.jl

@@ -1,15 +1,14 @@
 # # Nonlinear Boundary Value Problem
 
-# ## Non-linear differential equation
+# There is preliminary support for nonlinear equations, via Newton iteration in function space.
+# Here is a simple two-point boundary value problem:
+
 # We solve
 # ```math
 # Du = 0.001u^{\prime\prime} + 6(1-x^2)u^{\prime} + u^2 - 1 = 0,
 # ```
 # subject to the boundary conditions ``u(-1)=1`` and ``u(1)=-0.5``.
-# We collectively express the system as
-# ```math
-    # Nu = [u(-1)-1,\,u(1)+0.5,\,Du] = [0,0,0].
-# ```
+
 include("NonlinearBVP1.jl")
 
 # We plot the solution

examples/NonlinearBVP.jl

@@ -0,0 +1,33 @@
+# # Systems of equations
+
+# Systems of equations can be handled by creating a matrix of operators and functionals.
+# For example, we can solve the system
+
+# ```math
+# \begin{gathered}
+#     \mathop{u}'' - \mathop{u} + 2 \mathop{v} = \mathop{e}^x,  \\
+#     \mathop{u} + \mathop{v}' + \mathop{v} = \cos{x}, \\
+#     \mathop{u}(-1) = \mathop{u}'(-1) = \mathop{v}(-1) = 0
+# \end{gathered}
+# ```
+
+# using the following code:
+
+using ApproxFun
+using LinearAlgebra
+
+x = Fun();
+B = Evaluation(Chebyshev(),-1);
+A = [B      0;
+     B*𝒟    0;
+     0      B;
+     𝒟^2-I  2I;
+     I      𝒟+I];
+u,v = A \ [0;0;0;exp(x);cos(x)];
+
+import Plots
+Plots.plot(u, label="u", xlabel="x", legend=:topleft)
+Plots.plot!(v, label="v")
+
+# In this example, the automatic space detection failed and so we needed
+# to specify explicitly that the domain space for `B` is `Chebyshev()`.