Add harmonic oscillator examples (#879)

* Add harmonic oscillator examples

* Add tunnelling eigenvalue example

Author: jishnub

examples/Eigenvalue.jl

@@ -7,15 +7,30 @@ include("Eigenvalue_standard.jl")
 
 import Plots
 using LinearAlgebra: norm
-p = Plots.plot(V; legend=false, ylim=(-Inf, λ[22]))
+p = Plots.plot(V, legend=false, ylim=(-Inf, λ[22]))
 for k=1:20
-    Plots.plot!(real(v[k]/norm(v[k]) + λ[k]), )
+    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.
 
+include("Eigenvalue_tunnelling.jl")
+
+# We plot the first few eigenfunctions
+p = Plots.plot(V, legend=false)
+Plots.vline!([-Lx/2, Lx/2], linecolor=:black)
+p_twin = Plots.twinx(p)
+for k=1:6
+    Plots.plot!(p_twin, real(v[k]/norm(v[k]) + λ[k]), label="$k")
+end
+p
+
+# Note that the parity symmetry isn't preserved exactly at finite matrix sizes.
+# In general, it's better to preserve the symmetry of the operator matrices (see section below),
+# and projecting them on to the appropriate subspaces.
+
 # For problems with different contraints or boundary conditions,
 # `B` can be any zero functional constraint, e.g., `DefiniteIntegral()`.
 

examples/Eigenvalue_standard.jl

@@ -8,6 +8,9 @@
 # ```
 # where ``λ_i`` is the ``i^{th}`` eigenvalue of the ``L`` and has a corresponding
 # *eigenfunction* ``\mathop{u}_i(x)``.
+
+### Quantum harmonic oscillator
+
 # 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,

examples/Eigenvalue_tunnelling.jl

@@ -0,0 +1,19 @@
+# ### Tunnelling
+
+# We solve the Schrodinger equation in an infinite square well in the domain `-Lx/2..Lx/2`, with a barrier at the center from `-d/2..d/2`
+# ```math
+# \mathop{L} = -\frac{1}{2}\frac{\mathop{d}^2}{\mathop{dx}^2} + V(x),
+# ```
+# where we choose a smoothed barrier potential ``V(x)``.
+
+Lx = 4 # size of the domain
+d = 1 # size of the barrier
+Δ = 0.1 # smoothing scale of the barrier
+x = Fun(-Lx/2..Lx/2)
+V = 5 * (1 + tanh((x + d/2)/Δ)) * (1 - tanh((x - d/2)/Δ));
+H = -𝒟^2/2 + V;
+
+S = space(x)
+B = Dirichlet(S);
+
+λ, v = ApproxFun.eigs(B, H, 1000, tolerance=1E-8);

examples/ODE.jl

@@ -16,6 +16,34 @@ import Plots
 Plots.plot(u, xlabel="x", ylabel="u(x)", legend=false)
 
 
+# ## Initial value problem
+# ### Undamped harmonic oscillator
+# ```math
+# \frac{d^2 u}{dt^2} + 4 u = 0,
+# ```
+# in an interval `0..2pi`, subject to the initial conditions ``u(0)=0``
+# and ``u'(0)=2``.
+
+include("ODE_IVP.jl")
+
+# We plot the solution
+t = a:0.1:b
+Plots.plot(t, u.(t), xlabel="t", label="u(t)")
+Plots.plot!(t, sin.(2t), label="Analytical", seriestype=:scatter)
+
+# ### Damped harmonic oscillator
+# ```math
+# \frac{d^2 y}{dt^2} + 2\zeta\omega_0\frac{dy}{dt} + \omega_0^2 y = 0,
+# ```
+# from ``t=0`` to ``t=T``, subject to the initial conditions ``y(0)=4``
+# and ``y'(0)=0``.
+
+include("ODE_IVP2.jl")
+
+# Plot the solution
+import Plots
+Plots.plot(y, xlabel="t", ylabel="y(t)", legend=false)
+
 # ## Increasing Precision
 
 # Solving differential equations with high precision types is available.

examples/ODE_IVP.jl

@@ -0,0 +1,23 @@
+using ApproxFun
+using LinearAlgebra
+
+# Construct the domain
+a, b = 0, 2pi
+d = a..b;
+
+# We construct the differential operator ``L = d^2/dx^2 + 4``:
+D = Derivative(d)
+L = D^2 + 4;
+
+# We have not chosen the space explicitly, and the solve chooses it to be `Chebyshev(d)` internally
+
+# We incorporate the initial conditions using `ivp`, which is a shortcut for evaluating the function
+# and it's derivative at the left boundary of the domain
+A = [L; ivp()];
+
+# Initial conditions
+u0 = 0;
+dtu0 = 2;
+
+# We solve the differential equation
+u = \(A, [0,u0,dtu0], tolerance=1e-6);

examples/ODE_IVP2.jl

@@ -0,0 +1,18 @@
+# Construct the domain
+T = 12
+d = 0..T;
+
+Dt = Derivative(d);
+ζ = 0.2 # damping ratio
+ω0 = 2 # oscillation frequency
+L = Dt^2 + 2ζ * ω0 * Dt + ω0^2;
+
+# initial conditions
+y0 = 4; # initial displacement
+dty0 = 3; # initial velocity
+
+# The differential operator along with the initial condition evaluation operator
+A = [L; ivp()];
+
+# We solve the differential equation
+y = \(A, [0,y0,dty0], tolerance=1e-6);