Multiscale example in one dimension

This script applies the FEM to a one dimensional example of a multiscale problem. This problem was introduced by Peterseim in "Variational Multiscale Stabilization and the Exponential Decay of correctors, p.2".

$$ \begin{cases} - (A_{\varepsilon}(x)u_{\varepsilon}'(x))' &= 1, \qquad \text{ for }x \in (0,1)\\ u_{\varepsilon}(0)= u_{\varepsilon}(1) &= 0, \end{cases} $$

where, for $\varepsilon >0$, $$ A_{\varepsilon}(x) := \frac{1}{4}\left( 2 - \cos \left(\frac{2 \pi x}{\varepsilon}\right) \right)^{-1}. $$ The exact solution is given by $$ u_{\varepsilon}(x) = 4 (x-x^2)- 4 \varepsilon \left( \frac{1}{4 \pi} \sin(2 \pi \frac{x}{\varepsilon}) - \frac{1}{2 \pi}x \sin(2 \pi \frac{x}{\varepsilon}) - \frac{\varepsilon}{4 \pi^2} \cos(2 \pi \frac{x}{\varepsilon}) + \frac{\varepsilon}{4 \pi^2} \right). $$

In order to demonstrate the issue that comes with multiscale problems in terms of the FEM, we use several choices of the mesh size $h$ and compare the energy error.


In [1]:
import os
import sys
import numpy as np

%matplotlib notebook
import matplotlib.pyplot as plt

from gridlod import util, world, fem
from gridlod.world import World
import femsolverCoarse

First, we define the given functions and furthermore, we visualize the coefficient for two choices of epsilon.


In [2]:
fine = 4096
NFine = np.array([fine])
NpFine = np.prod(NFine+1)
NList = [2,4,8,16, 32, 64, 128, 256]

epsilon = 2**(-5)
epsilon1 = 2**(-6)

pi = np.pi
xt = util.tCoordinates(NFine).flatten()
xp = util.pCoordinates(NFine).flatten()

aFine = (2 - np.cos(2*pi*xt/epsilon))**(-1)
aFine1 = (2 - np.cos(2*pi*xt/epsilon1))**(-1)

uSol  = 4*(xp - xp**2) - 4*epsilon*(1/(4*pi)*np.sin(2*pi*xp/epsilon) -
                                    1/(2*pi)*xp*np.sin(2*pi*xp/epsilon) -
                                    epsilon/(4*pi**2)*np.cos(2*pi*xp/epsilon) +
                                    epsilon/(4*pi**2))

uSol = uSol/4

# plot the coefficient
plt.figure('Coefficient')
plt.plot(xt,aFine, label='$A_{\epsilon}(x)$')
plt.yticks((0,np.max(aFine)+np.min(aFine)),fontsize="small")
plt.ylabel('$y$', fontsize="small")
plt.xlabel('$x$', fontsize="small")
plt.legend(frameon=False,fontsize="large")
plt.title('$A_{\epsilon}(x)$ for $\epsilon=2^{-5}$.')
plt.show()

# plot the coefficient for a smaller epsilon
plt.figure('Coefficient1')
plt.plot(xt,aFine1, label='$A_{\epsilon}(x)$')
plt.yticks((0,np.max(aFine)+np.min(aFine1)),fontsize="small")
plt.ylabel('$y$', fontsize="small")
plt.xlabel('$x$', fontsize="small")
plt.legend(frameon=False,fontsize="large")
plt.title('$A_{\epsilon}(x)$ for $\epsilon=2^{-6}$.')
plt.show()


We compute the FEM approximation for each mesh size, plot the result and store the energy error. In order to apply the FEM, our computations are based on the 'gridlod' framework but with a coarse right hand side.


In [3]:
newErrorFine = []
x = []
y = []
for N in NList:
    NWorldCoarse = np.array([N])
    boundaryConditions = np.array([[0, 0]])

    NCoarseElement = NFine/NWorldCoarse
    world = World(NWorldCoarse, NCoarseElement, boundaryConditions)
    AFine = fem.assemblePatchMatrix(NFine, world.ALocFine, aFine)
    
    #grid nodes
    xpCoarse = util.pCoordinates(NWorldCoarse).flatten()
    NpCoarse = np.prod(NWorldCoarse+1)
    f = np.ones(NpCoarse)
    uCoarseFull = femsolverCoarse.solveCoarse_fem(world, aFine, f, boundaryConditions)
    
    basis = fem.assembleProlongationMatrix(NWorldCoarse, NCoarseElement)
    uLodCoarse = basis*uCoarseFull
    newErrorFine.append(np.sqrt(np.dot(uSol - uLodCoarse, AFine*(uSol - uLodCoarse))))
    x.append(N)
    y.append(1./N)
    if np.size(x)==1:
        plt.figure('FEM-Solutions')
        plt.subplots_adjust(left=0.01,bottom=0.04,right=0.99,top=0.95,wspace=0,hspace=0.2)
        plt.subplot(241)
        plt.plot(xp,uSol,'k', label='$u_{\epsilon}(x)$')
        plt.plot(xpCoarse,uCoarseFull,'o--', label= 'u_h(x)')
        plt.title('1/h= ' + str(N),fontsize="small")
        plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        plt.legend(frameon=False,fontsize="small")
    elif np.size(x)==2:
        plt.subplot(242)
        plt.plot(xp,uSol,'k', label='$u_{\epsilon}(x)$')
        plt.plot(xpCoarse,uCoarseFull,'o--', label= 'u_h(x)')
        plt.title('1/h= ' + str(N),fontsize="small")
        plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        plt.legend(frameon=False,fontsize="small")
    elif np.size(x)==3:
        plt.subplot(243)
        plt.plot(xp,uSol,'k', label='$u_{\epsilon}(x)$')
        plt.plot(xpCoarse,uCoarseFull,'o--', label= 'u_h(x)')
        plt.title('1/h= ' + str(N),fontsize="small")
        plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        plt.legend(frameon=False,fontsize="small")
    elif np.size(x)==4:
        plt.subplot(244)
        plt.plot(xp,uSol,'k', label='$u_{\epsilon}(x)$')
        plt.plot(xpCoarse,uCoarseFull,'o--', label= 'u_h(x)')
        plt.title('1/h= ' + str(N),fontsize="small")
        plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        plt.legend(frameon=False,fontsize="small")
    elif np.size(x)==5:
        plt.subplot(245)
        plt.plot(xp,uSol,'k', label='$u_{\epsilon}(x)$')
        plt.plot(xpCoarse,uCoarseFull,'--', label= 'u_h(x)')
        plt.title('1/h= ' + str(N),fontsize="small")
        plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        plt.legend(frameon=False,fontsize="small")
    elif np.size(x)==6:
        plt.subplot(246)
        plt.plot(xp,uSol,'k', label='$u_{\epsilon}(x)$')
        plt.plot(xpCoarse,uCoarseFull,'--', label= 'u_h(x)')
        plt.title('1/h= ' + str(N),fontsize="small")
        plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        plt.legend(frameon=False,fontsize="small")
    elif np.size(x)==7:
        plt.subplot(247)
        plt.plot(xp,uSol,'k', label='$u_{\epsilon}(x)$')
        plt.plot(xpCoarse,uCoarseFull,'--', label= 'u_h(x)')
        plt.title('1/h= ' + str(N),fontsize="small")
        plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        plt.legend(frameon=False,fontsize="small")
    elif np.size(x)==8:
        plt.subplot(248)
        plt.plot(xp,uSol,'k', label='$u_{\epsilon}(x)$')
        plt.plot(xpCoarse,uCoarseFull,'--', label= 'u_h(x)')
        plt.title('1/h= ' + str(N),fontsize="small")
        plt.tick_params(axis='both', which='both', bottom='off', top='off', labelbottom='off', right='off', left='off', labelleft='off')
        plt.legend(frameon=False,fontsize="small")
plt.show()


Lastly, we plot the energy error.


In [4]:
plt.figure("Error")
plt.loglog(x,newErrorFine,'o-', basex=2, basey=2)
plt.loglog(x,y,'--k',basex=2, basey=2, linewidth=1, alpha=0.3)
plt.ylabel('Energy error')
plt.xlabel('$1/h$')
plt.subplots_adjust(left=0.1,bottom=0.1,right=0.98,top=0.95,wspace=0.2,hspace=0.2)
plt.title('Energy error for the standard FEM')
plt.grid(True,which="both",ls="--")
plt.show()