Day 15

Example from Ackleh p547: check alleged exact solution

x = sp.symbols('x')
y = -100*x + 100*sp.log(1+x)/sp.log(2)
r = 1+x
s = 0
f = 100
de10_23 = -sp.diff( r*sp.diff(y,x), x) + s*y - f
sp.simplify(de10_23)

$\displaystyle 0$

Galerkin with pwlinear basis on example from Ackleh p547

from matplotlib import rcdefaults
rcdefaults()  # restore default matplotlib rc parameters
%config InlineBackend.figure_format='retina'   
import seaborn as sns  # wrapper for matplotlib that provides prettier styles and more
import matplotlib.pyplot as plt  # use matplotlib functionality directly
%matplotlib inline   
sns.set()

def make_phi(template,h,xi):
    def phi(x):
        return template( (x-xi)/h )
    return phi

def pwlinear_spline(x):
    f = np.zeros_like(x)
    subinterval = (x>=-1) & (x< 0); subx = x[subinterval]+1; f[subinterval] = subx
    subinterval = (x>= 0) & (x< 1); subx = 1-x[subinterval]; f[subinterval] = subx
    return f

fig, axs = plt.subplots(ncols=2, nrows=2, figsize=(15,10))
gs = axs[0,1].get_gridspec()
for ax in axs[:, -1]:
    ax.remove()
axbig = fig.add_subplot(gs[:, -1])

ns = []
errors = []
for np2 in [3,5,10,20,40,80]:  # various sizes of basis set

    # build the basis functions (pwlinear hat functions)
    gridx = np.linspace(0,1,np2)
    h = gridx[1] - gridx[0]
    phis = [make_phi(pwlinear_spline,h,xi) for xi in gridx[1:-1]]

    n = np2 - 2
    
    # integrals for right-hand side
    b = np.ones(n)*100*h

    # integrals for matrix elements
    A = np.zeros((n,n))
    def I(x): return (1+x)**2/2
    for i in range(1,n+1):
        A[i-1,i-1] = I(gridx[i+1]) - I(gridx[i-1])
        if i>1:
            A[i-1,i-2] = -I(gridx[i])   + I(gridx[i-1])
        if i<n:
            A[i-1,i-0] = -I(gridx[i+1]) + I(gridx[i])
    A /= h**2
    
    # solve for the coefficients
    a = np.linalg.solve(A,b)
    
    # make pictures
    x = np.linspace(0,1,500)
    approx = np.sum([ai*phi(x) for ai,phi in zip(a,phis)],axis=0)
    exact = -100*x + 100*np.log(1+x)/np.log(2)
    axs[0,0].plot(x, exact,'k')
    axs[0,0].plot(x, approx,alpha=0.5)
    axs[1,0].plot(x,approx-exact)
    ns.append(n)
    errors.append(np.abs(approx-exact).max())

axs[0,0].set_ylabel('exact and approximations')
axs[1,0].set_xlabel('x')
axs[1,0].set_ylabel('error(x)')
axbig.set_aspect(1)
axbig.loglog(ns,errors,'o')
axbig.set_xlabel('n')
axbig.set_ylabel('max error');

Plot on right has slope -2, showing max error goes to zero as $1/N^2$ or as $h^2$.