diff --git a/Week 6/feedback/2022-03-18 09:30:41.125777 UTC/10 Hyperbolic PDEs.html b/Week 6/feedback/2022-03-18 09:30:41.125777 UTC/10 Hyperbolic PDEs.html new file mode 100644 index 0000000..b6899ca --- /dev/null +++ b/Week 6/feedback/2022-03-18 09:30:41.125777 UTC/10 Hyperbolic PDEs.html @@ -0,0 +1,345997 @@ + + + + + +10 Hyperbolic PDEs + + + + + + + + + + + + + + + + + +
+
+
+

10 Hyperbolic PDEs (Score: 11.0 / 11.0)

+
+
    + + + + + + + + + + + +
  1. Coding free-response (Score: 0.0 / 0.0)
    2. + + + + + + + +
  2. Coding free-response (Score: 3.0 / 3.0)
    3. + + +
  3. Comment
    4. + + + + + + +
  4. Coding free-response (Score: 0.0 / 0.0)
    5. + + +
  5. Comment
    6. + + + +
  6. Coding free-response (Score: 3.0 / 3.0)
    7. + + + + + + + +
  7. Coding free-response (Score: 2.0 / 2.0)
    8. + + +
  8. Comment
    9. + + + + + + +
  9. Coding free-response (Score: 3.0 / 3.0)
    10. + + +
  10. Comment
    11. + + + +
  11. Coding free-response (Score: 0.0 / 0.0)
    12. + + + + +
+
+
+
+
+
+ +
+
+
+
+

CDS: Numerical Methods Assignments

    +

  • See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.

    +
    • +

  • Solutions must be submitted via the Jupyter Hub.

    +
    • +

  • Make sure you fill in any place that says YOUR CODE HERE or "YOUR ANSWER HERE".

    +
    • +
+

Submission

    +
  1. Name all team members in the the cell below
    2. +
  2. make sure everything runs as expected
    3. +
  3. restart the kernel (in the menubar, select Kernel$\rightarrow$Restart)
    4. +
  4. run all cells (in the menubar, select Cell$\rightarrow$Run All)
    5. +
  5. Check all outputs (Out[*]) for errors and resolve them if necessary
    6. +
  6. submit your solutions in time (before the deadline)
    7. +
+ +
+
+team_members = "Koen Vendrig, Kees van Kempen" +
+
+
+
+

Dynamic Behavior from Hyperbolic PDEs

In the following you will derive and implement finite-difference methods to study generalized hyperbolic partial differential equations.

+ +
+
+ +
+
+
In [1]:
+
Student's answer + Score: 0.0 / 0.0 (Top) +
+
+
+ 
import numpy as np
+from matplotlib import pyplot as plt
+from mpl_toolkits import mplot3d
+
+ +
+
+ +
+
+ +
+
+
+
+
+

Task 1: finite-difference hyperbolic PDE solver

Our aim is to implement a Python function to find the solution of the following PDE:

+\begin{align*} + \frac{\partial^2}{\partial t^2} u(x,t) + - \alpha^2 \frac{\partial^2}{\partial x^2} u(x,t) = 0, + \qquad + 0 \leq x \leq l, + \qquad + 0 \leq t +\end{align*}

with the boundary conditions

+\begin{align*} + u(0,t) = u(l,t) &= 0, + &&\text{for } t > 0 \\ + u(x,0) &= f(x) \\ + \frac{\partial}{\partial t} u(x,0) &= g(x) + &&\text{for } 0 \leq x \leq l. +\end{align*}

By approximating the partial derivatives with finite differences we can recast the problem into the following form

+\begin{align*} + \vec{w}_{j+1} = \mathbf{A} \vec{w}_{j} + - \vec{w}_{j-1}. +\end{align*}

Here $\vec{w}_j$ is a vector of length $m$ in the discretized spatial coordinate $x_i$ at time step $t_j$. The spatial coordinates are defined as $x_i = i h$, where $i=0,1,\dots,m-1$ and $h = l/(m-1)$. The time steps are defined as $t_j = j k$, where $j = 0, 1, \dots, n-1$.

+

The tri-diagonal matrix $\mathbf{A}$ has size $(m-2)\times(m-2)$ and is defined by

+\begin{align*} + \mathbf{A} = + \left( \begin{array}{cccc} + 2(1-\lambda^2) & \lambda^2 & & 0\\ + \lambda^2 & \ddots & \ddots & \\ + & \ddots & \ddots & \lambda^2 \\ + 0 & & \lambda^2 & 2(1-\lambda^2) + \end{array} \right), +\end{align*}

where $\lambda = \alpha k / h$. This $(m-2)\times(m-2)$ structure accounts for the first set of boundary conditions. Note that the product $\mathbf{A} \vec{w}_{j}$ is thus only performed over the $m-2$ subset, i.e. $i=1,2,\dots,m-2$. The other boundary conditions are accounted for by initializing the first two time steps with

+\begin{align*} + w_{i,j=0} &= f(x_i) \\ + w_{i,j=1} &= (1-\lambda^2) f(x_i) + + \frac{\lambda^2}{2} f(x_{i+1}) + + \frac{\lambda^2}{2} f(x_{i-1}) + + kg(x_i). +\end{align*}

Implement a Python function of the form $\text{pdeHyperbolic(a, x, t, f, g)}$, where $\text{a}$ represents the PDE parameter $\alpha$, $\text{x}$ and $\text{t}$ are the discretized spatial and time grids, and $\text{f}$ and $\text{g}$ are the functions defining the boundary conditions. This function should return a two-dimensional array $\text{w[:,:]}$, which stores the spatial vector $\vec{w}_j$ at each time coordinate $t_j$.

+ +
+
+ +
+
+
In [2]:
+
Student's answer + Score: 3.0 / 3.0 (Top) +
+
+
+ 
def pdeHyperbolic(a, x, t, f, g):
+    """
+    Numerically solves the hyperbolic differential equation ∂^2u(x,t)/∂t^2 - a*∂^2u(x,t)/∂x^2 = 0
+    with constant a for boundary conditions given by
+        u(0, t) = 0 = u(x[-1], t) for all t
+        y(x, 0) = f(x)
+        ∂u(x, t)/∂t = g(x) for t = 0 and all x
+
+    Args:
+        a: numerical constant in the PDE
+        x: array of evenly spaced space values x in the PDE
+        t: array of evenly spaced times values t in the PDE
+        f: callable function of numerical x giving a boundary condition to solution u of the PDE
+        g: callable function of numerical x giving a boundary condition to derivative ∂u/∂t of the PDE
+
+    Returns:
+        An |t| by |x| matrix w giving approximate solutions to the PDE over the grid
+        imposed by arrays x and t, such that w[j, i] corresponds to u[x[i], y[j]].
+    """
+    
+    n = len(t)
+    m = len(x)
+    
+    # TODO: The stepsize is defined by fixed h and, but the input x and t
+    #       could allow variable step size. What should it be?
+    #h = x[1] - x[0]
+    l = x[-1]
+    h = l/(m - 1)
+    k = t[1] - t[0]
+    λ = a*k/h
+    
+    # Create the tri-diagonal matrix A of size (m - 2) by (m - 2).
+    A = np.eye(m - 2)*2*(1 - λ**2) + ( np.eye(m - 2, m - 2, 1) + np.eye(m - 2, m - 2, -1) )*λ**2
+    
+    # Create empty matrix w for the result.
+    # The boundary values at x[0] and x[-1] are taken care of this way, too.
+    w = np.zeros((n, m))
+    
+    # Set initial values for w[0, i] and w[1, i] for all i except for the boundaries.
+    w[0] = f(x)
+    w[1, 1:m - 1] = (1 - λ**2)*f(x[1:m - 1]) + λ**2/2*f(x[2:m]) + λ**2/2*f(x[0:m - 2]) + k*g(x[1:m - 1])
+    
+    # Loop over all times t[j] to find w[j, i].
+    for j in range(2, n):
+        w[j, 1:m - 1] = np.dot(A, w[j - 1, 1:m - 1]) - w[j - 2, 1:m - 1]
+    
+    return w
+
+ +
+
+
+
+ Comments:

Please observe the correct ranges of variables in the boundary conditions.

+ +
+
+
+
+ +
+
+
+
+
+

Task 2

Use your implementation to solve the following problems. Compare in the first problem your numerical solution to the analytic one and use it to debug your code.

+

Problem 1:

\begin{align*} + \frac{\partial^2}{\partial t^2} u(x,t) + - \frac{\partial^2}{\partial x^2} u(x,t) &= 0, + &&\text{for }0 \leq x \leq 1 \text{ and } 0 \leq t \leq 1\\ + u(0,t) = u(l,t) &= 0, + &&\text{for } t > 0 \\ + u(x,0) &= \operatorname{sin}\left(2 \pi x \right), \\ + \frac{\partial}{\partial t} u(x,0) &= 2 \pi \operatorname{sin}\left(2 \pi x \right) + &&\text{for } 0 \leq x \leq 1. +\end{align*}

Compare your numerical solution to the analytic one and use it to debug your code. The corresponding analytic solution is

+\begin{equation} + u(x,t) = \operatorname{sin}(2 \pi x) (\operatorname{cos}(2 \pi t) + \operatorname{sin}(2 \pi t)). +\end{equation}

Then write a unit test for your function based on this problem.

+ +
+
+ +
+
+
In [3]:
+
Student's answer + Score: 0.0 / 0.0 (Top) +
+
+
+ 
# There is an 𝑙 in the boundary conditions we assume should be a 1.
+a1 = 1
+l1 = 1
+x1 = np.linspace(0, l1, 1200)
+t1 = np.linspace(0, 1, 1200)
+f1 = lambda x: np.sin(2*np.pi*x)
+g1 = lambda x: 2*np.pi*np.sin(2*np.pi*x)
+
+w1 = pdeHyperbolic(a1, x1, t1, f1, g1)
+u1 = lambda x, t: np.sin(2*np.pi*x)*(np.cos(2*np.pi*t) + np.sin(2*np.pi*t))
+
+# Plot the whole.
+fig = plt.figure(figsize=(24, 8))
+
+grid1 = np.meshgrid(x1, t1)
+
+# Plot w
+ax0 = fig.add_subplot(1, 3, 1, projection='3d')
+ax0.plot_surface(*grid1, w1, cmap="turbo")
+ax0.set_title('numerically found solution $w[t_j, x_i]$')
+ax0.set_zlabel('$w[t_j, x_i]$')
+
+# Plot the difference u - w
+ax1 = fig.add_subplot(1, 3, 2, projection='3d')
+ax1.plot_surface(*grid1, u1(*grid1) - w1, cmap="turbo")
+ax1.set_title('difference $u(x_i,t_j) - w[t_j, x_i]$')
+ax1.set_zlabel('$u(x_i,t_j) - w[t_j, x_i]$')
+
+# Plot u
+ax2 = fig.add_subplot(1, 3, 3, projection='3d')
+ax2.plot_surface(*grid1, u1(*grid1), cmap="turbo")
+ax2.set_title('exact solution $u(x_i,t_j)$')
+ax2.set_zlabel('$u(x_i, t_j)$')
+
+# Set shared settings
+for ax in [ax0, ax1, ax2]:
+    ax.set_xlabel('$x_i$')
+    ax.set_ylabel('$t_j$')
+
+fig.show()
+
+print("""
+On the left, the numerically found solution w is plotted against a grid of equidistant x and t values x_i and t_j.
+On the right, the exact solution 𝑢(𝑥,𝑡) = sin(2𝜋𝑥)(cos(2𝜋𝑡)+sin(2𝜋𝑡)) is plotted against the same grid.
+In the middle, their difference u - w is shown, again against the same grid.
+Do note the order of magnitude of the z axis in the middle as compared to the other plots. The error is quite small.
+Increasing the grid density lowers the errors.
+""")
+
+ +
+
+
+
+ Comments:

Nice figures!

+ +
+
+
+
+ +
+
+ + +
+ +
+ + +
+
+On the left, the numerically found solution w is plotted against a grid of equidistant x and t values x_i and t_j.
+On the right, the exact solution 𝑢(𝑥,𝑡) = sin(2𝜋𝑥)(cos(2𝜋𝑡)+sin(2𝜋𝑡)) is plotted against the same grid.
+In the middle, their difference u - w is shown, again against the same grid.
+Do note the order of magnitude of the z axis in the middle as compared to the other plots. The error is quite small.
+Increasing the grid density lowers the errors.
+
+
+
+
+ +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
In [4]:
+
Student's answer + Score: 3.0 / 3.0 (Top) +
+
+
+ 
def test_pdeHyperbolic():
+    """
+    Check that function pdeHyperbolic returns solutions for the PDE
+    ∂^2/∂𝑡^2 𝑢(𝑥,𝑡) − ∂^2/∂𝑥^2 𝑢(𝑥,𝑡) = 0 for boundary conditions given by
+        𝑢(0, t) = 0 = 𝑢(x[-1], t) for all 0 <= t <= 1
+        y(x, 0) = f(x) for all 0 <= x <= 1
+        ∂𝑢(x, t)/∂t = g(x) for t = 0 and all 0 <= x <= 1
+    that are within tolerance (TOL = 1e-5) of the exact solution
+    𝑢(𝑥,𝑡) = sin(2𝜋𝑥)(cos(2𝜋𝑡)+sin(2𝜋𝑡)),
+    using a grid of x and t values of each 1200 equidistant points.
+    """
+    
+    a = 1
+    l = 1
+    x = np.linspace(0, l, 1200)
+    t = np.linspace(0, 1, 1200)
+    f = lambda x: np.sin(2*np.pi*x)
+    g = lambda x: 2*np.pi*np.sin(2*np.pi*x)
+
+    w = pdeHyperbolic(a, x, t, f, g)
+    u = lambda x, t: np.sin(2*np.pi*x)*(np.cos(2*np.pi*t) + np.sin(2*np.pi*t))
+    grid = np.meshgrid(x, t)
+
+    TOL = 1e-5
+    assert np.all(np.abs(u(*grid) - w) < TOL)
+    
+test_pdeHyperbolic()
+
+ +
+
+ +
+
+ +
+
+
+
+
+

Problem 2:

\begin{align*} + \frac{\partial^2}{\partial t^2} u(x,t) + - \frac{\partial^2}{\partial x^2} u(x,t) &= 0, + &&\text{for }0 \leq x \leq 1 \text{ and } 0 \leq t \leq 2\\ + u(0,t) = u(l,t) &= 0, + &&\text{for } t > 0 \\ + u(x,0) &= \left\{ \begin{array}{cc} +1 & \text{for } x<0.5 \\ -1 & \text{for } x \geq 0.5 \end{array} \right., \\ + \frac{\partial}{\partial t} u(x,0) &= 0 + &&\text{for } 0 \leq x \leq 1. +\end{align*}

Use $m=200$ and $n=400$ to discretize the spatial and time grids, respectively.

+ +
+
+ +
+
+
In [5]:
+
Student's answer + Score: 2.0 / 2.0 (Top) +
+
+
+ 
# There is an 𝑙 in the boundary conditions we assume should be a 1.
+a2 = 1
+l2 = 1
+m2 = 200
+n2 = 400
+x2 = np.linspace(0, l2, m2)
+t2 = np.linspace(0, 2, n2)
+f2 = lambda x: -1 + 2.*(x < .5)
+g2 = lambda x: 0
+
+w2 = pdeHyperbolic(a2, x2, t2, f2, g2)
+
+# Plot the whole.
+fig = plt.figure(figsize=(16, 16))
+
+grid2 = np.meshgrid(x2, t2)
+
+# Plot w
+ax = fig.add_subplot(1, 1, 1, projection='3d')
+ax.plot_surface(*grid2, w2, cmap="turbo")
+ax.set_title('numerically found solution $w[t_j, x_i]$')
+ax.set_zlabel('$w[t_j, x_i]$')
+ax.set_xlabel('$x_i$')
+ax.set_ylabel('$t_j$')
+
+fig.show()
+
+# TODO: Are comments really needed?
+
+ +
+
+
+
+ Comments:

The comments are sometimes needed to show the depth of your understanding of your observations over the obtained results! +From that figure is very hard to see the patterns of x at fixed t. You can consider 2D plots or surface plots for better clarity.

+ +
+
+
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
+
+

Task 3

Animate your solutions! To this end you can use the following code:

+
# use matplotlib's animation package
+import matplotlib.pylab as plt
+import matplotlib
+import matplotlib.animation as animation
+# set the animation style to "jshtml" (for the use in Jupyter)
+matplotlib.rcParams['animation.html'] = 'jshtml'
+
+# create a figure for the animation
+fig = plt.figure()
+plt.grid(True)
+plt.xlim( ... )     # fix x limits
+plt.ylim( ... )     # fix y limits
+
+# Create an empty plot object and prevent its showing (we will fill it each frame)
+myPlot, = plt.plot([0], [0])
+plt.close()
+
+# This function is called each frame to generate the animation (f is the frame number)
+def animate(f):                   
+    myPlot.set_data( ... )  # update plot
+
+# Show the animation
+frames = np.arange(1, np.size(t))  # t is the time grid here
+myAnimation = animation.FuncAnimation(fig, animate, frames, interval = 20)
+myAnimation
+
+ +
+
+ +
+
+
In [6]:
+
Student's answer + Score: 3.0 / 3.0 (Top) +
+
+
+ 
# use matplotlib's animation package
+import matplotlib.pylab as plt
+import matplotlib
+import matplotlib.animation as animation
+# set the animation style to "jshtml" (for the use in Jupyter)
+matplotlib.rcParams['animation.html'] = 'jshtml'
+
+# create a figure for the animation
+fig = plt.figure()
+plt.grid(True)
+plt.xlim(0, 1)     # fix x limits
+plt.ylim(-1.5, 1.5)     # fix w limits
+plt.xlabel('$x_i$')
+plt.ylabel('$w[t_j, x_i]$')
+
+# Create an empty plot object and prevent its showing (we will fill it each frame)
+myPlot, = plt.plot([0], [0])
+plt.close()
+
+# This function is called each frame to generate the animation (f is the frame number)
+def animate(f):                   
+    myPlot.set_data(x1, w1[f])  # update plot
+
+# Show the animation
+frames = np.arange(1, np.size(t1))  # t is the time grid here
+myAnimation = animation.FuncAnimation(fig, animate, frames, interval = 20)
+myAnimation
+
+ +
+
+
+
+ Comments:

You can also animate both the numerical and the analytical solutions at the same time.

+ +
+
+
+
+ +
+
+ + +
+ +
Out[6]:
+ + + +
+ + + + + + +
+ +
+ +
+ + + + + + + + + +
+
+ + + + + + +
+
+
+ + + + +
+ +
+ +
+
+ +
+
+
+
In [7]:
+
Student's answer + Score: 0.0 / 0.0 (Top) +
+
+
+ 
# use matplotlib's animation package
+import matplotlib.pylab as plt
+import matplotlib
+import matplotlib.animation as animation
+# set the animation style to "jshtml" (for the use in Jupyter)
+matplotlib.rcParams['animation.html'] = 'jshtml'
+
+# create a figure for the animation
+fig = plt.figure()
+plt.grid(True)
+plt.xlim(0, 1)     # fix x limits
+plt.ylim(-1.5, 1.5)     # fix w limits
+plt.xlabel('$x_i$')
+plt.ylabel('$w[t_j, x_i]$')
+
+# Create an empty plot object and prevent its showing (we will fill it each frame)
+myPlot, = plt.plot([0], [0])
+plt.close()
+
+# This function is called each frame to generate the animation (f is the frame number)
+def animate(f):                   
+    myPlot.set_data(x2, w2[f])  # update plot
+
+# Show the animation
+frames = np.arange(1, np.size(t2))  # t is the time grid here
+myAnimation = animation.FuncAnimation(fig, animate, frames, interval = 20)
+myAnimation
+
+ +
+
+ +
+
+ +
+
+ + +
+ +
Out[7]:
+ + + +
+ + + + + + +
+ +
+ +
+ + + + + + + + + +
+
+ + + + + + +
+
+
+ + + + +
+ +
+ +
+
+ +
+
+
+
+
+
+ + diff --git a/Week 6/feedback/2022-03-18 09:30:41.125777 UTC/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.html b/Week 6/feedback/2022-03-18 09:30:41.125777 UTC/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.html new file mode 100644 index 0000000..83c3b37 --- /dev/null +++ b/Week 6/feedback/2022-03-18 09:30:41.125777 UTC/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.html @@ -0,0 +1,13739 @@ + + + + + +11 Parabolic PDEs: 1D Schrödinger Equation with Potential + + + + + + + + + + + + + + + + + +
+
+
+

11 Parabolic PDEs: 1D Schrödinger Equation with Potential (Score: 6.0 / 15.0)

+
+
    + + + + + + + + + + + +
  1. Coding free-response (Score: 0.0 / 0.0)
    2. + + + + + + + +
  2. Written response (Score: 1.0 / 1.0)
    3. + + +
  3. Comment
    4. + + + +
  4. Coding free-response (Score: 1.0 / 3.0)
    5. + + +
  5. Comment
    6. + + + + + + +
  6. Coding free-response (Score: 2.0 / 3.0)
    7. + + +
  7. Comment
    8. + + + + + + +
  8. Written response (Score: 1.0 / 2.0)
    9. + + +
  9. Comment
    10. + + + + + + +
  10. Coding free-response (Score: 1.0 / 3.0)
    11. + + +
  11. Comment
    12. + + + +
  12. Coding free-response (Score: 0.0 / 3.0)
    13. + + + + + + + +
  13. Coding free-response (Score: 0.0 / 0.0)
    14. + + + + + + + +
  14. Coding free-response (Score: 0.0 / 0.0)
    15. + + + + +
+
+
+
+
+
+ +
+
+
+
+

CDS: Numerical Methods Assignments

    +

  • See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.

    +
    • +

  • Solutions must be submitted via the Jupyter Hub.

    +
    • +

  • Make sure you fill in any place that says YOUR CODE HERE or "YOUR ANSWER HERE".

    +
    • +
+

Submission

    +
  1. Name all team members in the the cell below
    2. +
  2. make sure everything runs as expected
    3. +
  3. restart the kernel (in the menubar, select Kernel$\rightarrow$Restart)
    4. +
  4. run all cells (in the menubar, select Cell$\rightarrow$Run All)
    5. +
  5. Check all outputs (Out[*]) for errors and resolve them if necessary
    6. +
  6. submit your solutions in time (before the deadline)
    7. +
+ +
+
+team_members = "Koen Vendrig, Kees van Kempen" +
+
+
+
+

Parabolic PDEs: 1D Schrödinger Equation with Potential

Here we aim to solve the 1D (parabolic) Schrödinger equation (SEQ) for a particle exposed to some potential $V(x)$

+\begin{align*} + i \hbar \frac{\partial}{\partial t} \Psi(x, t) = \frac{-\hbar^2}{2m} \frac{\partial^2}{\partial x^2} \Psi(x, t) + V(x) \Psi(x, t). +\end{align*}

For simplicity we will set $\hbar = 1$ and $m=0.5$ in the following. First you will analyze the stationary solutions of the SEQ for a given potential, using finite differences to approximate all involved derivatives. Then you will need to derive and apply the Crank-Nicolson approach for this PDE to study the dynamics of the SEQ.

+ +
+
+ +
+
+
In [1]:
+
Student's answer + Score: 0.0 / 0.0 (Top) +
+
+
+ 
import numpy as np
+from matplotlib import pyplot as plt
+
+ +
+
+ +
+
+ +
+
+
+
+
+

Task 1

The stationary SEQ

+\begin{align*} + H \Psi(x) = \frac{-\hbar^2}{2m} \frac{\partial^2}{\partial x^2} \Psi(x) + V(x) \Psi(x) = E \Psi(x) +\end{align*}

can be numerically studied by discretizing the spatial coordinate $x$ as $x_i = ih$ and approximating the second derivative with finite differences. This results in an eigenvalue problem of the form $H \vec{\Psi} = E \vec{\Psi}$, where $\Psi$ is now a vector in the discretized positions $x_i$. By diagonalizing the matrix $H$ we simultaneously get both the eigenenergies $E_i$ and the corresponding eigenfunctions $\vec{\Psi}_i$.

+

Derive the matrix representation of the Hamiltonian $H$ and write your final result in the cell below. (Double click on "YOUR ANSWER HERE" to open the cell, and ctrl+enter to compile.)

+

Implement a simple Python function $\text{SEQStat(x, V)}$ which takes the discretized grid $x_i$ and the potential $V(x)$ as arguments and which returns all eigenvalues $E$ and eigenfunctions $\vec{\Psi}$.

+ +
+
+ +
+
+
+
Student's answer + Score: 1.0 / 1.0 (Top) +
+
+
+

The potential $V(x)$ is a scalar quantity, and will therefore only occur on the diagonal of the matrix $H$. +The second derivative operator will however be matrix valued, and we will determine it by looking at the finite differences approximation.

+$$ +\frac{\partial^2\Psi}{\partial{}x^2}(x_i) = \frac{\Psi(x_{i+1}) - 2\Psi(x_i) + \Psi(x_{i-1})}{h^2} - \frac{h^2}{12} \frac{\partial^4\Psi}{\partial x^4}(x_i) +$$

For small $h$, we can ignore the latter term, which we will indeed do. +We can encode this in a tridiagonal matrix, as each step depends on the step before and after it and not on anything else.

+

Our beautiful hamiltonian now has elements like

+$$ +H_{ij} = \frac{-\hbar^2}{2m} \left[ \delta_{i,j-1} - 2\delta_{ij} + \delta_{i,j+1} \right] + V(x) \delta_{ij} = -\left[ \delta_{i,j-1} - 2\delta_{ij} + \delta_{i,j+1} \right] + V(x) \delta_{ij} +$$

with $\delta_{ij} = 1 \iff i = j$, and using the approximations $\hbar = 1$, $m = 0.5$.

+$$ +%finite differences for second derivative to get a tridiagonal matrix (one lower and one higher than diagonal) +%then plot eigenvalues +%no time component -> not complex? +%energy should be real +%eigenvectors squared should be normalized +$$ +
+
+
+
+ Comments:

If you use finite differences, you need to provide the step size h for the numerical implementations. (-0.5 pt).

+ +
+
+
+ +
+
+
In [2]:
+
Student's answer + Score: 1.0 / 3.0 (Top) +
+
+
+ 
def SEQStat(x, V):
+    """
+    Numerically solves the stationary Schrödinger equation over
+    the grid of x values, in a potential given by V, using
+    the approximations ℏ = 1, 𝑚 = 0.5.
+    
+    Args:
+        x: array of evenly spaced space values x
+        V: function that takes in a location x and returns a scalar potential value
+
+    Returns:
+        A tuple of the eigenvalues E and eigenfunctions \Psi of the hamiltonian H.
+    """
+    
+    n = len(x)
+    h = x[1] - x[0]
+    
+    # In the stationary case, 
+    H = 1/h**2*(np.eye(n, n, -1) - 2*np.eye(n) + np.eye(n, n, 1)) + np.eye(n)*V(x) 
+    
+    # In principle, the eigenvalues and eigenvectors could be complex, so
+    # we opt for eigh instead of eig.
+    return np.linalg.eigh(H)
+
+ +
+
+
+
+ Comments:

Please check the superdiagonal and subdiagonal matrix elements they should be multiplied by -1. Also, check the main diagonal matrix elements (-2 pt) since 2/h^2 should be positive.

+ +
+
+
+
+ +
+
+
+
+
+

Task 2

Use your function to calculate and plot the first few eigenfunctions for $V(x) = 0.25 x^2$. Add the potential to your plot. Use $-5 \leq x \leq +5$ discretized in $100$ steps.

+ +
+
+ +
+
+
In [3]:
+
Student's answer + Score: 2.0 / 3.0 (Top) +
+
+
+ 
V = lambda x: 0.25*x**2
+x = np.linspace(-5, 5, 100)
+the_first_few = 4
+
+E, Psi = SEQStat(x, V)
+
+fig = plt.figure(figsize=(12,12))
+
+for n in range(the_first_few):
+    # Normalize the found eigenfunctions using Riemann midpoint sums
+    prob_mass = np.sum((x[1:] - x[:-1])*(Psi[1:, n]**2 + Psi[:-1, n]**2)/2)
+    Psi[:, n] /= np.sqrt(prob_mass)
+    
+    # Confirm normalization
+    prob_mass = np.sum((x[1:] - x[:-1])*Psi[:-1, n]**2)
+    print("Eigenstate n =", n, "yields total probability mass =", prob_mass)
+    
+    # Plot the eigenfunctions
+    ax = fig.add_subplot(the_first_few, 2, 2*n + 1)
+    ax.set_title("Eigenstate n = " + str(n))
+    ax.plot(x, Psi[:, n], label="Eigenstate n = " + str(n))
+    ax.set_xlabel("$x$")
+    ax.set_ylabel("$\Psi$")
+    
+    # Plot the squared eigenfunctions to give an idea about probability
+    ax = fig.add_subplot(the_first_few, 2, 2*n + 2)
+    ax.plot(x, Psi[:, n]**2, label="Eigenstate n = " + str(n))
+    ax.set_xlabel("$x$")
+    ax.set_ylabel("$\Psi^2$")
+
+plt.tight_layout()
+plt.show()
+
+ +
+
+
+
+ Comments:

The figures are nice! However, please check the comments in task 1 regarding the diagonal elements of the Hamiltonian matrix. (-1 pt)

+ +
+
+
+
+ +
+
+ + +
+ +
+ + +
+
Eigenstate n = 0 yields total probability mass = 0.9999999999999999
+Eigenstate n = 1 yields total probability mass = 1.0
+Eigenstate n = 2 yields total probability mass = 1.0
+Eigenstate n = 3 yields total probability mass = 1.0
+
+
+
+ +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
+
+

Task 3

Now we turn to the full dynamical problem, which can be numerically analyzed using the Crank-Nicolson approach, which averages over forward and backward scattering finite differences for the time derivative. Your task is to derive the corresponding formulas for the 1D SEQ from above, implement it, and apply it to study the time-dynamics for the 1D SEQ. Applying the Crank-Nicolson approach to the 1D SEQ from above results in an equation system of the form

+\begin{align*} + \mathbf{A} \vec{w}_{j+1} = \mathbf{B} \vec{w}_{j} +\end{align*}

where $\vec{w}_j$ is the space- and time-discretized approximation to $\Psi(x_i, t_j)$ with $j$ being the time-discretization index. $\mathbf{A}$ and $\mathbf{B}$ are $(m-2)\times(m-2)$ matrices of the form

+\begin{align*} +\mathbf{A} = + \left( \begin{array}{cccc} + 1 + \lambda_1 + \lambda_2 V_1& -\frac{\lambda_1}{2} & & 0\\ + -\frac{\lambda_1}{2} & \ddots & \ddots & \\ + & \ddots & \ddots & -\frac{\lambda_1}{2} \\ + 0 & & -\frac{\lambda_1}{2} & 1 + \lambda_1 + \lambda_2 V_{m-1} + \end{array} \right) +\end{align*}

and

+\begin{align*} +\mathbf{B} = + \left( \begin{array}{cccc} + 1 - \lambda_1 - \lambda_2 V_1& +\frac{\lambda_1}{2} & & 0\\ + +\frac{\lambda_1}{2} & \ddots & \ddots & \\ + & \ddots & \ddots & +\frac{\lambda_1}{2} \\ + 0 & & +\frac{\lambda_1}{2} & 1 - \lambda_1 - \lambda_2 V_{m-1} + \end{array} \right). +\end{align*}

Derive this form and thereby the definitions of $\lambda_1$ and $\lambda_2$. Start with recapping the derivation of the Crank-Nicolson approach (see Wikipedia or Numerical Analysis book by Burden and Faires).

+

Write your definitions of $\lambda_1$ and $\lambda_2$ in the Markdown cell below. (Double click on "YOUR ANSWER HERE" to open the cell, and ctrl+enter to compile.)

+ +
+
+ +
+
+
+
Student's answer + Score: 1.0 / 2.0 (Top) +
+
+
+

To do the derivation of the Crank-Nicolson approach following $\textit{Numerical Analysis}$ by Burden and Faires, we start from averaging our 'Forward-Difference method' at the jth step in t:

+$$ +\frac{w_{i,j+1}-w_{i,j}}{k}-\alpha^2\frac{w_{i+1,j}-2w_{i,j}+w_{i-1,j}}{h^2} = 0. +$$

The local truncation error is:

+$$ +\tau_F = \frac{k}{2}\frac{\partial^2 u}{\partial t^2}(x_i,\mu_j)+O(h^2). +$$

We do the same for the 'Backward-Difference method' to get:

+$$ +\frac{w_{i,j+1}-w_{i,j}}{k}-\alpha^2\frac{w_{i+1,j+1}-2w_{i,j+1}+w_{i-1,j+1}}{h^2} = 0, +$$

again with a local truncation error:

+$$ +\tau_B = -\frac{k}{2}\frac{\partial^2 u}{\partial t^2}(x_i,\hat{u} _j)+O(h^2). +$$

If we make the following assumption:

+$$ +\frac{\partial^2 u}{\partial t^2}(x_i,\hat{\mu}_j) \approx \frac{\partial^2 u}{\partial t^2}(x_i,\mu_j), +$$

Then we can get the 'averaged-difference method' to write:

+$$ +\frac{w_{i,j+1}-w_{i,j}}{k}-\frac{\alpha^2}{2}\left[\frac{w_{i+1,j}-2w_{i,j}+w_{i-1,j}}{h^2}+\frac{w_{i+1,j+1}-2w_{i,j+1}+w_{i-1,j+1}}{h^2}\right] = 0, +$$

which now has an truncation error of the order $O(k^2+h^2)$. This is known as the Crank-Nicolson method and is represented in the form described above.

+

For $\lambda_1$ and $\lambda_2$ we have the following expressions:

+$$ +\lambda_1 = \alpha^2\frac{k}{h^2}\\ +\lambda_2 = 1. +$$ +
+
+
+
+ Comments:

The derivation has small errors in it. First, the definitions of \lambda_1 and \lambda_2 aren't correct! The imaginary number is missing too from the definition (-1 pt). You started with two terms and ended up with three please check that.

+ +
+
+
+ +
+
+
+
+

Task 4

The equation system from above can be used to calculate the time propagation by multiplying the whole expression from the left by $\mathbf{A}^{-1}$. The only needed ingredients are the boundary conditions $w_{i=0,j} = w_{i=m-1,j} = 0$ for all $j$ (not explicitly needed, just used for the derivation of $\mathbf{A}$ and $\mathbf{B}$) and the initial value $w_{i,j=0} = f(x_i)$.

+

Implement the Crank-Nicolson approach for our 1D SEQ in a Python function $\text{SEQDyn(x, t, V, f)}$, where $\text{x}$ and $\text{t}$ are the discretized spatial and time grids, $\text{V}$ is the potential, and $\text{f}$ is the function defining the boundary conditions. The function should return a two-dimensional array, which stores the eigenfunction $\Psi(x_i, t_j)$ at all discretized spatial and time coordinates.

+

Apply it to the SEQ from task 2.2, with $V(x) = 0.25 x^2$ and $f(x) = \operatorname{exp}(-x^2) / 2.5$. Use $-5 \leq x \leq +5$ with $100$ steps and $0 \leq t \leq 50$ with $1000$ steps. Animate your solution (see exercise 9).

+ +
+
+ +
+
+
In [4]:
+
Student's answer + Score: 1.0 / 3.0 (Top) +
+
+
+ 
def SEQDyn(x, t, V, f):
+    """
+    Numerically solves the time-dependent Schrödinger equation over
+    the gives input parameters, in a potential given by V, using
+    the approximations ℏ = 1, 𝑚 = 0.5, using the Crank-Nicolson approach
+    
+    Args:
+        x: array of evenly spaced space values x
+        t: array of evenly spaced time values t
+        V: function that takes a location x and returns a scalar potential value
+        f: function that takes a location x and returns initial values to w for
+           the Crank-Nicolson approach
+
+    Returns:
+        A tuple of eigenfunctions \Psi of the time-dependent SE.
+    """
+    
+    n = len(x)
+    h = x[1] - x[0]
+    m = len(t)
+    k = t[1] - t[0]
+    
+    # Setup the Crank-Nicholson parameters
+    α  =
+    λ1 = α**2*k/h**2
+    λ2 = 1
+    
+    A  = -λ1/2*( np.eye(m - 2, m - 2, -1) + np.eye(m - 2, m - 2, 1) ) - (1 + λ1 + λ2*V)*np.eye(m - 2)
+    B  = +λ1/2*( np.eye(m - 2, m - 2, -1) + np.eye(m - 2, m - 2, 1) ) - (1 + λ1 + λ2*V)*np.eye(m - 2)
+
+    # TODO: Finish SEQDyn: implement algorithm, return the right things.
+
+
+V = lambda x: 0.25*x**2
+f = lambda x: np.exp(-x**2)/2.5
+x = np.linspace(-5, 5, 100)
+t = np.linspace(0, 50, 1000)
+
+ +
+
+
+
+ Comments:

The function implementation is missing \alpha (-2 pt). Also, the function has a problem. Please check the previous comments.

+ +
+
+
+
+ +
+
+ + +
+ +
+ + +
+
+  File "/tmp/ipykernel_53509/1568846030.py", line 24
+    α  =
+        ^
+SyntaxError: invalid syntax
+
+
+
+ +
+
+ +
+
+
+
In [5]:
+
Student's answer + Score: 0.0 / 3.0 (Top) +
+
+
+ 
# Animate your solution here ...
+
+# YOUR CODE HERE
+
+ +
+
+ +
+
+ +
+
+
+
+
+

[Optional] Task 5

What happens when your initial condition is set to one of the eigenfunctions you obtained from the static 1D SEQ?

+ +
+
+ +
+
+
In [6]:
+
Student's answer + Score: 0.0 / 0.0 (Top) +
+
+
+ 
# YOUR CODE HERE
+
+ +
+
+ +
+
+ +
+
+
+
+
+

[Optional] Task 6

Change the potential to a staggered one of the form

+\begin{align*} + V(x) = \left\{ + \begin{array}{cc} + +15 & -4.0 < x < -3.5, \, -2.5 < x < -2.0, \, -1.0 < x < -0.5, \\ + +15 & +0.5 < x < +1.0, \, +2.0 < x < +2.5, \, +3.5 < x < +4.0, \\ + 0 & \text{else} + \end{array} \right. +\end{align*} +
+
+ +
+
+
In [7]:
+
Student's answer + Score: 0.0 / 0.0 (Top) +
+
+
+ 
def V(x):
+    if -4.0 < x < -3.5 or -2.5 < x < -2.0 or -1.0 < x < -0.5:
+        return 15.
+    if  0.5 < x <  1.0 or  2.0 < x <  2.5 or  3.5 < x <  4.0:
+        return 15.
+    return 0
+
+ +
+
+ +
+
+ +
+
+
+
+
+
+ +