p30: Spectral integration, ODE style

p30: Spectral integration, ODE style#

Compare to p12

%config InlineBackend.figure_format='svg'
from numpy import zeros,abs,dot,exp,sqrt,pi
from scipy.special import erf
from numpy.linalg import inv
from chebPy import cheb
from matplotlib.pyplot import figure,subplot,semilogy,grid,xlim,ylim,xlabel,ylabel,text

Nmax = 50; E = zeros((4,Nmax+1))
for N in range(1,Nmax+1):
    D, x = cheb(N)
    x = x[1:]; Di = inv(D[1:,1:]); w = Di[-1,:]
    f = abs(x)**3;             E[0,N] = abs(dot(w,f)-0.5)
    f = exp(-(x+1e-20)**(-2)); E[1,N] = abs(dot(w,f) - 2*(exp(-1)+sqrt(pi)*(erf(1)-1)))
    f = 1/(1+x**2);            E[2,N] = abs(dot(w,f) - pi/2)
    f = x**10;                 E[3,N] = abs(dot(w,f) - 2.0/11.0)

labels = ('$|x|^3$', '$exp(-x^{-2})$', '$1/(1+x^2)$', '$x^{10}$')
figure(figsize=(8,8))
for iplot in range(4):
    subplot(3,2,iplot+1)
    semilogy(range(1,Nmax+1), E[iplot,1:]+1e-100, '.-')
    xlabel('N'); ylabel('error'); xlim(1,Nmax); ylim(1e-18,1e3); grid(True)
    text(25,0.02,labels[iplot])
_images/1aa5932ac573a1f5227d466c6a7f04014d9a61d19733a5894171c20bae3b72f1.svg

Repeat with Clenshaw-Curtis#

from clencurt import *
Nmax = 50; E = zeros((4,Nmax+1))
for N in range(1,Nmax+1):
    x,w = clencurt(N+1)
    f = abs(x)**3;             E[0,N] = abs(dot(w,f)-0.5)
    f = exp(-(x+1e-20)**(-2)); E[1,N] = abs(dot(w,f) - 2*(exp(-1)+sqrt(pi)*(erf(1)-1)))
    f = 1/(1+x**2);            E[2,N] = abs(dot(w,f) - pi/2)
    f = x**10;                 E[3,N] = abs(dot(w,f) - 2.0/11.0)

labels = ('$|x|^3$', '$exp(-x^{-2})$', '$1/(1+x^2)$', '$x^{10}$')
figure(figsize=(8,8))
for iplot in range(4):
    subplot(3,2,iplot+1)
    semilogy(range(1,Nmax+1), E[iplot,1:]+1e-100, '.-')
    xlabel('N'); ylabel('error'); xlim(1,Nmax); ylim(1e-18,1e3); grid(True)
    text(25,0.02,labels[iplot])
_images/d29163bd95cfe99180310b1521b2e8c62da4c256b9df8112c769369b2343a26b.svg

Repeat with Gauss-Legendre#

from gauss import *
Nmax = 50; E = zeros((4,Nmax+1))
for N in range(1,Nmax+1):
    x,w = gauss(N+1)
    f = abs(x)**3;             E[0,N] = abs(dot(w,f)-0.5)
    f = exp(-(x+1e-20)**(-2)); E[1,N] = abs(dot(w,f) - 2*(exp(-1)+sqrt(pi)*(erf(1)-1)))
    f = 1/(1+x**2);            E[2,N] = abs(dot(w,f) - pi/2)
    f = x**10;                 E[3,N] = abs(dot(w,f) - 2.0/11.0)

labels = ('$|x|^3$', '$exp(-x^{-2})$', '$1/(1+x^2)$', '$x^{10}$')
figure(figsize=(8,8))
for iplot in range(4):
    subplot(3,2,iplot+1)
    semilogy(range(1,Nmax+1), E[iplot,1:]+1e-100, '.-')
    xlabel('N'); ylabel('error'); xlim(1,Nmax); ylim(1e-18,1e3); grid(True)
    text(25,0.02,labels[iplot])
_images/fec68d34fb61cb6a895a8a34929cf3f2404675dc0f4926ccc077f5559f65bd02.svg
Contents