Problem 1 Poisson equation (Sparse matrices)

Consider a 2D Poisson equation $$ \Delta u \equiv \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f, \quad (x,y)\in [0,1]^2 $$ with boundary conditions as follows $$ u\big|_{x=0} = 0, \quad u\big|_{x=1} = 0, \quad u\big|_{y=0} = 0, \quad u\big|_{y=1} = 0, $$ with known function $f$ and unknown $u$. Choose $f(x,y)$ such that $\sin\pi x \sin \pi y $ is a solution.

  • (2 pts) Introduce a uniform grid with $N$ steps in each direction. Discretize equation via the 5-point stencil. Write down matrix $A_h$ of the linear system using Kronecker product properties ($h$ is a mesh step)
  • (1 pts) Solve the system via the direct sparse solver (use standard functions in Python). Plot the solution
  • (2 pts) Plot error as a function of $h$. Find numerically order of convergence
  • (3 pts) Solve the system via an appropriate Krylov method with and without ILU preconditioner (use standard functions in Python). Compare their timings
  • (2 pts) Compare timings for the direct solver and the preconditioned iterative Krylov method. Vary $N$ up to $512$
  • (bonus) Find analytically eigenvalues of the matrix $A_h$ and prove that $\text{cond}( A_h )= \mathcal{O}\left(\frac{1}{h^2}\right)$

In [1]:

Problem 2 First kind Fredholm equation (Structured matrices)

Consider a 1D first kind Fredholm equation $$ \left( Au \right)(x)\equiv \int_0^1 K(x,y) u(y) dy = f(x),\quad x\in [0,1] \quad (1) $$ with known functions $K(x,y)$, $f(x)$ and unknown $u(x)$.

  • (2 pts) Set $K(x,y) = \sqrt{|x-y|}$ and choose $f(x)$ such that $u(x) = 1$ is a solution. Discretize this equation on a uniform grid by Nystrom method. What type of matrix does this linear system have? Note: Nystrom method: $\int_0^1 K(x_i,y) u(y) dy \approx \sum_j h K(x_i,y_j) u(y_j)$, where $x_i, y_j$ are points of 1D grids and $h$ is a grid step
  • (3 pts) Implement matvec function which multiplies Toeplitz matrix by vector with $\mathcal{O}(N \log N)$ complexity.
  • (3 pts) By giving this function instead of the whole matrix to an appropriate Krylov solver in Python find the solution of (1) and compare timings with Krylov solver without fast matvec (vary $N$ up to 2048). Note: set tolerance level $10^{-9}$ in the Krylov solver
  • (1 pt) Add random noise with amplitude $\delta = 10^{-4}h$ to the right hand side $f$ and plot the solution. Explain the result
  • (2 pts) Regularize (1) by truncating singular values. Plot solutions for different truncation levels via widgets animation (find example below). Comment on the results

In [3]:
#http://nbviewer.ipython.org/github/ipython/ipython/blob/master/examples/Interactive%20Widgets/Using%20Interact.ipynb
%matplotlib inline

from IPython.html.widgets import interact
from IPython.display import display

import matplotlib.pyplot as plt
import numpy as np

fig, ax = plt.subplots()
plt.close(fig)

def set_cursor(x):
    ax.clear()
    ax.set_xlim([0, np.pi])
    ax.set_ylim([-1, 1])
    x1 = np.linspace(0,np.pi,50)
    ax.plot(x1, np.sin(x * x1))
    display(fig)

interact(set_cursor, x=(1., 2., 0.1))


Out[3]:
<function __main__.set_cursor>

Problem 3 Matrix exponential

Consider the following time dependent equation in 2D $$ \frac{\partial u}{\partial t} = \Delta u, \quad (x,y)\in \mathbb{R}^2, \quad t\in [0,1] $$ $$ u(x,y)\big|_{t=0} = \exp{(- x^2 - y^2)}, $$ with the boundary conditions (here we replace the whole $\mathbb{R}^2$ by $[-10,10]^2$ square ) $$ u\big|_{x=-10} = 0, \quad u\big|_{x=10} = 0, \quad u \big|_{y=-10} = 0, \quad u\big|_{y=10} = 0, $$ Spatial discretization looks as follows $$ \frac{du_h}{dt} = A_h u_h, $$ where $A_h\in \mathbb{R}^{256^2 \times 256^2}$ is a matrix from the Problem 1. To solve this equation one can use several approaches:

Explicit scheme (Euler scheme) $$ \frac{u_h^{i+1} - u_h^{i}}{\tau} = A_h u_h^{i}, \quad i=1,\dots,N_\tau, \quad N_\tau \tau = 1 $$

Implicit scheme $$ \frac{u_h^{i+1} - u_h^{i}}{\tau} = A_h u_h^{i+1}, \quad i=1,\dots,N_\tau, \quad N_\tau \tau = 1 $$

Matrix exponential $$ u_h(t) = e^{At}u_h(0) $$

  • (1 pt) Name advantages and disadvantages of the decribed methods
  • (2 pts) Implement the explicit scheme. Plot $u_h(1)$ for several $\tau$ (stable and unstable cases). Find maximum possible $\tau$ when the scheme is still stable
  • (2 pts) Implement the implicit scheme. Plot $u_h(1)$ for several $\tau$. Is this scheme absolute stable?
  • (1 pt) Find $u_h(1)$ via matrix exponential ($\verb|expm_multiply|$ in Python). Plot $u_h(1)$
  • (1 pt) Compare timings for these 3 schemes

In [ ]: