Day 15¶

Example from Ackleh p547: check alleged exact solution

In [3]:
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)

Out[3]:
$\displaystyle 0$

Galerkin with pwlinear basis on example from Ackleh p547¶

In [4]:
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()

In [209]:
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()

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$.

In [ ]: