Poisson example using Literate (#822)

Author: jishnub

docs/src/usage/equations.md

@@ -93,31 +93,3 @@ u(0.1)
 
 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
 as it is applied to different right-hand sides, so the same operator can be inverted much more efficiently in subsequent problems.
-
-## Partial differential equations
-
-Partial differential operators are also supported.  Here's an example
-of solving the Poisson equation with zero boundary conditions:
-
-```julia
-d = Domain(-1..1)^2
-x,y = Fun(d)
-f = exp.(-10(x+0.3)^2-20(y-0.2)^2)  # use broadcasting as exp(f) not implemented in 2D
-A = [Dirichlet(d);Δ]  # Δ is an alias for Laplacian()
-@time u = A \ [zeros(∂(d));f]  # 4s for ~3k coefficients
-```
-
-Using a QR Factorization reduces the cost of subsequent calls substantially:
-
-```julia
-QR = qr(A)
-@time QR \ [zeros(∂(d));f]  # 4s
-g = exp.(-10(x+0.2)^2-20(y-0.1)^2)
-@time QR \ [zeros(∂(d));g]  # 0.09s
-```
-
-Many PDEs have weak singularities at the corners, in which case it is beneficial to specify a tolerance to reduce the time:
-
-```julia
-\(A, [zeros(∂(d));f]; tolerance=1E-6)
-```

examples/PDE.jl

@@ -1,5 +1,38 @@
 # # Solving PDEs
 
+# ## Poisson equation
+
+# We solve the Poisson equation
+# ```math
+# Δu(x,y) = 0
+# ```
+# on the rectangle `-1..1 × -1..1`,
+# with zero boundary conditions:
+
+using ApproxFun
+using LinearAlgebra
+
+d = (-1..1)^2
+x,y = Fun(d)
+f = exp.(-10(x+0.3)^2-20(y-0.2)^2)  # use broadcasting as exp(f) not implemented in 2D
+A = [Dirichlet(d); Laplacian()]
+u = A \ [zeros(∂(d)); f];
+
+# Using a QR Factorization reduces the cost of subsequent calls substantially
+
+QR = qr(A)
+u = QR \ [zeros(∂(d)); f];
+
+import Plots
+Plots.plot(u; st=:surface, legend=false)
+
+# Many PDEs have weak singularities at the corners,
+# in which case it is beneficial to specify a `tolerance` to reduce the time, as:
+
+\(A, [zeros(∂(d)); f]; tolerance=1E-6);
+
+# ## Helmholtz Equation
+
 # We solve
 # ```math
 # (Δ + 100)u(x,y) = 0

examples/Sampling.jl

@@ -11,3 +11,4 @@ import Plots
 Plots.histogram(x; bins=100, normalize = :pdf,
     title = "Histogram", xlabel="value", ylabel="frequency",
     xlim = (-5,5), ylim=(0, Inf), legend = false)
+Plots.plot!(f/sum(f), linewidth=3)