1. Mandatory assignment due 11/10-2024

• Please organize you project in a fork of the repository at: MATMEK-4270/mandatory1. You should start by changing the badge in the README to reflect your own work.

• Collaboration/discussion among students is encouraged, but each student need an individual fork of the repository. Please send me an email about your github username if it is not obvious who you are from this username. I can see all the forks, but I cannot necessarily understand who is who.

• You need a green badge in order to pass the mandatory assignment. You also need to answer a few questions in a report and to place a movie (see Section 1.2.6) where I can easily see it on github.com (without downloading).

1.1. Poisson’s equation

Consider Poisson’s equation

(1.1)\[ \begin{equation} \nabla^2 u = f, \end{equation} \]

for the two-dimensional domain \(\Omega = [0, L]^2\). Use a uniform mesh of equal length for both directions, such that the mesh can be described as \(\boldsymbol{x} \times \boldsymbol{y}\), where

\[ \begin{equation} \boldsymbol{x} = (i h)_{i=0}^{N}, \quad \boldsymbol{y} = (j h)_{j=0}^{N} \end{equation} \]

and \(h = L/N\). There is already the start of an implementation of a solver in Poisson2D. The solver is created such that it takes a Sympy function ue as input when creating the solver class (see Poisson2D.__init__). This manufactured solution ue is used to create the right hand side function \(f(x, y)\), such that when you solve the Poisson problem numerically the resulting mesh function \(u_{ij}\) should converge towards ue.

1.1.1. Implement the Poisson solver

Modify the solver class Poisson2D such that the two tests test_convergence_poisson2d and test_interpolation passes.

1.2. Stationary waves

Consider the wave equation

(1.2)\[ \begin{equation} \frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u, \end{equation} \]

in a two-dimensional spatial domain \(\Omega = [0, 1]^2\) and in time \(t \in [0, T]\). Use a uniform mesh for the finite difference method in both space and time. Use the same spatial discretization in both \(x\) and \(y\)-directions, such that the mesh can be described as in Section 1.1 for \(L=1\). Time is discretized as \(t_n = n \Delta t\) for \(n \in 0, 1, \ldots, N_t\), and \(\Delta t = T/N_t\). The CFL number \(C\) is defined as

\[ \begin{equation} C = \frac{c \Delta t}{h}. \end{equation} \]

A mesh function will be denoted as

\[ u^n_{ij} = u(t_n, x_i, y_j) \]

and for a given timestep \(n\) it is useful to use matrix notation for the mesh function:

\[ U^n = (u^n_{ij})_{i,j=0}^N. \]

All internal points in the computational domain can be discretized using central differences in both time and space

(1.3)\[ \begin{equation} \frac{u^{n+1}_{i,j} - 2u^n_{i,j} + u^{n-1}_{i, j}}{\Delta t^2} = c^2 \left(\frac{u^n_{i+1,j} - 2u^n_{i,j} + u^n_{i-1, j}}{h^2} + \frac{u^n_{i,j+1} - 2u^n_{i,j} + u^n_{i, j-1}}{h^2}\right) \end{equation} \]

The wave equation admits two real stationary waves as solutions:

The Dirichlet problem:

(1.4)\[ \begin{equation} u(t, x, y) = \sin(k_x x) \sin(k_y y) \cos(\omega t). \end{equation} \]

The Neumann problem:

(1.5)\[ \begin{equation} u(t, x, y) = \cos(k_x x) \cos(k_y y) \cos(\omega t). \end{equation} \]

Here \(k_x=m_x\pi\) and \(k_y=m_y\pi\), where \(m_x\) and \(m_y\) are arbitrary integers, and \(\omega\) is the dispersion coefficient.

1.2.1. The Dirichlet problem

There is the beginning of a solver class Wave2D for the homogeneous Dirichlet problem using the standing wave in Eq. (1.4) as initial condition.

Hint

A solver for the Dirichlet problem was implemented in lecture 7.

Complete the implementation of the Wave2D solver and verify the implementation using the analytical solution in (1.4). Compute the error using the \(\ell^{2}\) norm

\[ E = \|e^n\|_{\ell^2} = \sqrt{h^2 \sum_{i=0}^N \sum_{j=0}^N (e_{ij}^{n})^2} \]

where \(e^n = U^n - U^n_e\), \(U^n\) is the solution at time \(n \Delta t\) and \(U^n_e\) here is the analytical solution Eq. (1.4) at the mesh at the same time. The solver is verified when test_convergence_wave2d passes.

Hint

Use either (1.5) or (1.4) as an exact solution of (1.2) and find the dispersion coefficient \(\omega\) as a function of \(c, k_x\) and \(k_y\). Implement the solution in the property method w.

1.2.2. The Neumann problem

There is also the beginning of a solver for the wave equation with Neumann boundary conditions in Wave2D_Neumann. This solver should use the standing wave in Eq. (1.5) as initial condition and exact solution. Modify the three class methods that are included in order to make test_convergence_wave2d_neumann pass.

For the remaining questions, please answer in a report placed in the report folder. You can use a notebook or a PDF.

1.2.3. Exact solution

The two stationary solutions to the wave equations are real and imaginary components of the waves

(1.6)\[ \begin{equation} u(t,x,y) = e^{\imath (k_x x + k_y y - \omega t)} \end{equation} \]

where the imaginary unit \(\imath = \sqrt{-1}\). Show that Eq. (1.6) satisfies the wave equation.

1.2.4. Dispersion coefficient

Assume that \(m_x=m_y\) such that \(k_x=k_y=k\). A discrete version of Eq. (1.6) will then read

(1.7)\[ \begin{equation} u^n_{ij} = e^{\imath (kh(i+j) - \tilde{\omega} n\Delta t)} \end{equation} \]

where \(\tilde{\omega}\) is a numerical dispersion coefficient, i.e., the numerical approximation of the exact \(\omega\). Insert Eq. (1.7) into the discretized Eq. (1.3) and show that for CFL number \(C=1/\sqrt{2}\) we get \(\tilde{\omega} = \omega\).

1.2.5. Test for exact solution

Modify the test test_exact_wave2d such that it verifies that \(m_x=m_y\) and \(C=1/\sqrt{2}\) leads to l2-errors less than \(10^{-12}\) for both the Dirichlet problem and the Neumann problem.

1.2.6. Create movie

Create a movie of the Neumann problem using \(m_x=m_y=2\) and \(C = 1/\sqrt{2}\). Call the move neumannwave.gif and save it to the report folder of your forked repository. Make sure that the movie is not larger than 1 MB.

Hint

See Lecture 7 about how to create an animation and a movie in Python.