Exercise 5



In [1]:

    
import numpy
from matplotlib import pyplot
% matplotlib inline



In [2]:

    
import os, sys
sys.path.append(os.path.split(os.path.split(os.getcwd())[0])[0])



In [3]:

    
import utils.grids.one_d as assembly
import utils.quadrature as quad



In [4]:

    
import warnings
warnings.filterwarnings('error')

Consider the following problem:

$$ \frac{d^2 u}{d x^2} - \lambda u = f $$

where $x$ is defined as $x \in [-1, 1]$. The weak form with test function $\nu$ before lifting is

$$ \int_{-1}^{1}\frac{d\nu}{dx}\frac{du}{dx} dx + \lambda\int_{-1}^{1}\nu u dx = \nu(1)\frac{du}{dx}(1) - \nu(-1)\frac{du}{dx}(-1) - \int_{-1}^{1} \nu f dx $$

For continuous Galerkin method, the matrix form of the problem is:

$$ \left(\mathbb{L} + \lambda\mathbb{M}\right)\mathbf{u} = - \mathbf{f} + \mathbf{BC_N} $$

$\mathbf{BC_N}$ represents Newmann type boundary conditions, and

$$ \mathbf{BC_N} = \left\{ \begin{array}{ll} \left.-\frac{du}{dx}\right|_{x=-1} & \text{, if i represents the node }x=-1 \\ \left.\frac{du}{dx}\right|_{x=1} & \text{, if i represents the node }x=1 \\ 0 & \text{, else} \end{array} \right. $$

If any Dirchlet BC is given at any boundary point, then $\frac{du}{dx}$ is normally unknown at that point. However, after lifting the problem, the unknown $\frac{du}{dx}$ are eliminated. So those unknown $\frac{du}{dx}$ will not show up in the matrix system.

Here we define some wrapper functions first for our convenience.



In [5]:

    
def solver(N, P, f, DBCs, ends):
    """solvie 1D Helmholts PDE"""
    
    g = assembly.DecomposeAssembly(N, ends, "CommonJacobi", P)
    
    rhs = - g.weak_rhs(f)
    A = g.wL + g.M

    rhs -= A[:, g.l2g[0][0]].A.flatten() * DBCs[0]
    rhs -= A[:, g.l2g[-1][-1]].A.flatten() * DBCs[1]
    
    A = numpy.delete(A, [g.l2g[0][0], g.l2g[-1][-1]], 0)
    A = numpy.delete(A, [g.l2g[0][0], g.l2g[-1][-1]], 1)
    rhs = numpy.delete(rhs, [g.l2g[0][0], g.l2g[-1][-1]])
    
    ui = numpy.linalg.solve(A, rhs)
    ui = numpy.insert(ui, g.l2g[0][0], DBCs[0])
    ui = numpy.insert(ui, g.l2g[-1][-1], DBCs[1])
    g.set_coeffs(ui)
    
    return g



In [6]:

    
def H1(g, exact, d_exact, l):
    """calculate energy (H1) norm"""
    
    qd = quad.GaussLobattoJacobi(15)
    ans = 0
    for e in g.elems:
        f = lambda x: (e.derivative(x) - d_exact(x))**2 + (e(x) - exact(x))**2
        ans += numpy.sqrt(qd(f, xmin=e.ends[0], xMax=e.ends[1]))
    return ans / numpy.sqrt(l)



In [7]:

    
def ex_wrapper(N, P, f, exact, d_exact, ends):
    """exercise wrapper"""
    
    Ndof = {"h": numpy.zeros_like(N, dtype=numpy.int),
            "p": numpy.zeros_like(P, dtype=numpy.int)}
    err = {"h": numpy.zeros_like(N, dtype=numpy.float64),
           "p": numpy.zeros_like(P, dtype=numpy.float64)}
    
    DBCs = [exact(ends[0]), exact(ends[1])]
    L = ends[1] - ends[0]
    
    for i, Ni in enumerate(N):
        g = solver(Ni, 1, f, DBCs, ends)
        Ndof["h"][i] = g.nModes - 2  # values on two boundary nodes are know
        err["h"][i] = H1(g, exact, d_exact, L)
     
    for i, Pi in enumerate(P):
        g = solver(2, Pi, f, DBCs, ends)
        Ndof["p"][i] = g.nModes - 2  # values on two boundary nodes are know
        err["p"][i] = H1(g, exact, d_exact, L)
        
    return Ndof, err

Case 1

Now we consider a all-Dirichlet BCs case (with $\lambda=1$):

$$ f(x) = -(\pi^2 + \lambda)\sin(\pi x) $$

where the exact solution is $u(x) = \sin(\pi x)$

Plot the convergence order for both h- and p-type expansion.



In [8]:

    
N = numpy.array([2 + 4 * i for i in range(25)])
P = numpy.array([1, 3, 5, 7, 9, 11, 13, 15])

f = lambda x: - (1 + numpy.pi * numpy.pi) * numpy.sin(numpy.pi * x)
exact = lambda x: numpy.sin(numpy.pi * x)
d_exact = lambda x: numpy.pi * numpy.cos(numpy.pi * x)

Ndof, err = ex_wrapper(N, P, f, exact, d_exact, [-1, 1])



In [9]:

    
pyplot.loglog(Ndof["h"], err["h"], 'ko-', lw=1.5, label="h-type expansion")
pyplot.loglog(Ndof["p"], err["p"], 'k^-', lw=1.5, label="p-type expansion")
pyplot.title("Convergence of h- and p-type expansion, \nlog-log plot, and " + 
             r"$u(x)=\sin(\pi x)$")
pyplot.xlabel(r"$N_{dof}$")
pyplot.ylabel(r"Energy norm, $\left|\left|error\right|\right|_E$")
pyplot.legend(loc=0);



In [10]:

    
pyplot.semilogy(Ndof["h"], err["h"], 'ko-', lw=1.5, label="h-type expansion")
pyplot.semilogy(Ndof["p"], err["p"], 'k^-', lw=1.5, label="p-type expansion")
pyplot.title("Convergence of h- and p-type expansion, \n semi-log plot, and " + 
             r"$u(x)=\sin(\pi x)$")
pyplot.xlabel(r"$N_{dof}$")
pyplot.ylabel(r"Energy norm, $\left|\left|error\right|\right|_E$")
pyplot.legend(loc=0);

Case 2

Now we consider cases with non-smooth exact solutions: $u(x) = x^{\alpha}$, where $\alpha=0.6,\ 0.9,\ 1.2$. And the right-hand-side will be $f(x) = \alpha(\alpha-1)x^{\alpha-2}-x^{\alpha}$. We still use all-Dirichlet BCs here.

Plot the convergence order for both h- and p-type expansion.

Due to the singularity of the RHS is at $x=0$, so the integration for weak-form RHS on the first element will be Gauss-Radau-Jacobi quadrature (with $x=1$ included). Hence we have to modify the code.



In [11]:

    
def solver(N, P, f, DBCs, ends):
    """solvie 1D Helmholts PDE"""

    g = assembly.DecomposeAssembly(N, ends, "CommonJacobi", P)
    
    rhs = numpy.zeros(g.nModes, dtype=numpy.float64)
    rhs[g.l2g[0]] -= g.elems[0].weak_rhs(f, 25, "GaussRadauJacobi", end=1)
    for i, e in enumerate(g.elems[1:]):
        rhs[g.l2g[i+1]] -= e.weak_rhs(f, 20)
        
    A = g.wL + g.M

    rhs -= A[:, g.l2g[0][0]].A.flatten() * DBCs[0]
    rhs -= A[:, g.l2g[-1][-1]].A.flatten() * DBCs[1]
    
    A = numpy.delete(A, [g.l2g[0][0], g.l2g[-1][-1]], 0)
    A = numpy.delete(A, [g.l2g[0][0], g.l2g[-1][-1]], 1)
    rhs = numpy.delete(rhs, [g.l2g[0][0], g.l2g[-1][-1]])
    
    ui = numpy.linalg.solve(A, rhs)
    ui = numpy.insert(ui, g.l2g[0][0], DBCs[0])
    ui = numpy.insert(ui, g.l2g[-1][-1], DBCs[1])
    g.set_coeffs(ui)
    
    return g



In [12]:

    
def H1(g, exact, d_exact, l):
    """calculate energy (H1) norm"""
    
    qd = quad.GaussRadauJacobi(25, end=1)
    f = lambda x: ((g.elems[0].derivative(x) - d_exact(x))**2 +
                   (g.elems[0](x) - exact(x))**2)
    ans = numpy.sqrt(qd(f, xmin=g.elems[0].ends[0], xMax=g.elems[0].ends[1]))

    qd = quad.GaussLobattoJacobi(20)
    for e in g.elems[1:]:
        f = lambda x: (e.derivative(x) - d_exact(x))**2 + (e(x) - exact(x))**2
        ans += numpy.sqrt(qd(f, xmin=e.ends[0], xMax=e.ends[1]))
    return ans / numpy.sqrt(l)



In [13]:

    
N = numpy.array([2 + 4 * i for i in range(25)])
P = numpy.array([1 + 2 * i for i in range(11)])

Alpha = [0.6, 0.9, 1.2]
Ndof = numpy.empty(3, dtype=dict)
err = numpy.empty(3, dtype=dict)

for i, alpha in enumerate(Alpha):
    exact = lambda x: x**alpha
    d_exact = lambda x: alpha * (x**(alpha-1))
    f = lambda x: alpha * (alpha - 1) * (x**(alpha-2)) - x**alpha
    Ndof[i], err[i] = ex_wrapper(N, P, f, exact, d_exact, [0, 1])



In [14]:

    
pyplot.loglog(Ndof[0]["h"], err[0]["h"], 'ko-', lw=1.5, label=r"$\alpha=0.6$")
pyplot.loglog(Ndof[1]["h"], err[1]["h"], 'k^-', lw=1.5, label=r"$\alpha=0.9$")
pyplot.loglog(Ndof[2]["h"], err[2]["h"], 'ks-', lw=1.5, label=r"$\alpha=1.2$")
pyplot.title("Convergence of h-expansion, \nlog-log plot, and " + 
             r"$u(x)=x^{\alpha}$")
pyplot.xlabel(r"$N_{dof}$")
pyplot.ylabel(r"Energy norm, $\left|\left|error\right|\right|_E$")
pyplot.legend(loc=0);



In [15]:

    
pyplot.loglog(Ndof[0]["p"], err[0]["p"], 'ko-', lw=1.5, label=r"$\alpha=0.6$")
pyplot.loglog(Ndof[1]["p"], err[1]["p"], 'k^-', lw=1.5, label=r"$\alpha=0.9$")
pyplot.loglog(Ndof[2]["p"], err[2]["p"], 'ks-', lw=1.5, label=r"$\alpha=1.2$")
pyplot.title("Convergence of p-expansion, \nlog-log plot, and " + 
             r"$u(x)=x^{\alpha}$")
pyplot.xlabel(r"$N_{dof}$")
pyplot.ylabel(r"Energy norm, $\left|\left|error\right|\right|_E$")
pyplot.legend(loc=0);