Update examples

Extras/Dual Numbers.jl → Extras/Dual_Numbers.jl

@@ -3,8 +3,10 @@
 # changes in parameters
 ##
 
+using DualNumbers
+using ApproxFun
+using LinearAlgebra
 
-using DualNumbers, ApproxFun
 # What is the derivative of the function (without differentiating the Chebyshev expansion)?
 f=Fun(x->exp(dual(x,1)),-1..1)
 dualpart(f)

Extras/GalerkinOrthogonalPolynomials.jl

@@ -1,4 +1,6 @@
-using ApproxFun, Plots
+using ApproxFun
+using Plots
+using LinearAlgebra
 
 # In this example, we numerically compute with Galerkin orthogonal polynomials.
 # If H is a separable Hilbert space and B contains ν linear functional constraints,

Extras/Inverse_eigenvalue_problem.jl

@@ -0,0 +1,60 @@
+using ApproxFun, LinearAlgebra, DualNumbers
+import ApproxFun: bandwidth
+
+#
+# Consider the following Sturm--Liouville problem with polynomial p, q, and w.
+#
+# [(-𝒟)(p𝒟) + q] u = λ w u,  (u±bu')(±1) = 0.
+#
+# For which value of b, the Neumann part of the Robin boundary conditions,
+# is λ₀ equal to 1?
+#
+# To answer this question, we use a dual number to observe the sensitivity in λ₀
+# with respect to perturbations in b. The `realpart` of the objective function
+#
+#    f(λ(b)) = λ(b) - 1,
+#
+# extracts the difference, whereas the `dualpart` evaluates to ∂f/∂b (λ(b)).
+#
+
+function main()
+    n = 100
+    λ = dual(0.0, 0.0)
+    b = dual(0.0, 1.0)
+    v = zeros(Dual{Float64}, n)
+    v[1] = 1
+    v .+= 0.003.*randn(Float64, n)
+    for _ in 1:8
+        d = Segment(-1..1)
+        S = Ultraspherical(0.5, d)
+        NS = NormalizedPolynomialSpace(S)
+        D = Derivative(S)
+        p = Fun(x->1+x^2, S)
+        q = Fun(x->x^3, S)
+        L = -D*(p*D)+q
+        w = Fun(x->1+x/4, S)
+        M = Multiplication(w, S)
+        C = Conversion(domainspace(L), rangespace(L))
+        B = [Evaluation(S, -1, 0) - b*Evaluation(S, -1, 1); Evaluation(S, 1, 0) + b*Evaluation(S, 1, 1)]
+        QS = QuotientSpace(B)
+        Q = Conversion(QS, S)
+        D1 = Conversion(S, NS)
+        D2 = Conversion(NS, S)
+        R = D1*Q
+        P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
+        A = R'D1*P*L*D2*R
+        B = R'D1*M*D2*R
+        SA = Symmetric(A[1:n,1:n], :L)
+        SB = Symmetric(B[1:n,1:n], :L)
+        for _ = 1:7
+            global λ = dot(v, SA*v)/dot(v, SB*v)
+            ldiv!(ldlt!(SA-λ*SB), v)
+            global v = v/norm(v)
+            println("\tThis is λ: ", λ)
+        end
+        global b -= (realpart(λ)-1)/dualpart(λ)
+        println("This is b: ", b)
+    end
+end
+
+main()