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/64fa889a6c4c0748db3450ac1ad763b1e51a05639965ea944db371117e90f39c.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/0dbef423f0d61c3c696af2b84be6565c1aa0ae4ab60f4835d3bbad6f8215604f.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/01e0692d06bbc28cbe2c05d7e1c5f484dcf03b81ef5587b3f7a42e303dab8366.svg
Contents