Cheb2 demo

Little demo how to use cheb2mats and Chebyshev nodes to approximate functions with polynomials, for Chebyshev polynomials of the 2nd kind

Some of the outputs in cheb2mats come from the advanced toolbox chebfun by Trefethen et al, some are didactic implementations from their accompanying papers and manuals. http://www.chebfun.org Chebfun is completely object oriented permitting one to treat functions like vectors and ODEs like algebraic equations.

Contents

Investigate functions on interval

addpath('tools')
clear;format compact;
order=11;                          % order of polynomials to be used
interval=[0,1];                    % interval for BVP

Obtain basic quantities

[xi,D,intw,cdat]=cheb2mats(order,interval,'ynodes',1:order+1);
xd=linspace(interval(1),interval(2),200); % dense set of test points

General interpolation matrix:

for every (vectorized) function f:R->R, Ex*f(xi') approximates f(xd) on the dense grid, similarly Ex*D*f(xi') approximates the derivative on the dense grid, and intw*f(xi') approximates the integral of f over the interval.

Ex=cheb2evaljac(xd,cdat);

Test: Function representation and manipulation

f=@(u)1./(1+16*u.^2);              % example
df=@(u)-(32*u)./(16*u.^2 + 1).^2;  % its derivative
yi=f(xi);                          % interpolation values
fxd=f(xd);         % true function values for comparison
dfxd=df(xd);       % true derivatives for comparison

Plot interpolation points

and compare visually to true values and derivatives

figure(1);clf;
lw={'linewidth',2};
axprops=[lw,{'fontsize',16}];
plot(xd,Ex*yi',xd,fxd,'k--',xi,yi,'o',...
    xd,Ex*D*yi',xd,dfxd','k:',xi,D*yi','s',...
    lw{:},'markersize',8);
xlabel('x');
ylabel('f');
legend({'interpolation Ex*yi''','true values','Chebyshev-2 interpolation points',...
    'interp. derivative Ex*D*yi''',' true derivative','Derivative at interpolation points'},...
    'location','best');
set(gca,axprops{:});

accuracy of integral

quadrature_error=intw*yi'-atan(4)/4 %#ok<*NOPTS>

quadrature_error =
   7.0073e-08

Investigate effect of order systematically

ordrg=2:100;
figure(2);clf;
for i=1:length(ordrg)
    [xi,D,intw,cdat]=cheb2mats(ordrg(i),interval,'ynodes',1:ordrg(i)+1);
    Ex=cheb2evaljac(xd,cdat);
    pe=semilogy(ordrg(i),norm(Ex*f(xi')-f(xd'),'inf'),'ko',lw{:});
    pd=semilogy(ordrg(i),norm(Ex*D*f(xi')-df(xd'),'inf'),'ro',lw{:});
    pq=semilogy(ordrg(i),abs(intw*f(xi')-atan(4)/4),'bs',lw{:});
    hold on
end
hold off
xlabel('order');
ylabel('error');
legend([pe,pd,pq],{'interpolation error','interpolation error for derivative',...
    'error of quadrature'},...
    'location','best');
grid on
set(gca,axprops{:});

Use polynomial to solve a linear ODE x''+cx'+kx=exp(-t^2/sigma^2)*a, x(-2)=x(2)=0

cheb2mats returns data for multi-dimensional problems with the optinoal 'dim' argument.

k=1;
c=0.1;
A=[0,1;-k,-c];
a=1;
sigma=0.5;
rhs=@(t)a*exp(-t.^2/sigma^2);
interval=[-2,2];
order=100;
[xi,D,~,cdat]=cheb2mats(order,interval,'ynodes',1:order+1,'dim',2);

The ODE restricted to the polynomial

is D*y=A*y+[0;rhs(t)]. A needs to be multiplied in every point xi such that the linear part has to be blown up with kron(eye(size(xi)),A). Note that the data (D,xi,cdat) are the same for every ODE with dimension 2 on the given interval. We replace the ODE on the first node xi(1) by the boundary conditions.

Eode=D-kron(eye(order+1),A);
Ebc=[1,0,0,0;0,0,1,0]*cheb2evaljac(interval,cdat);
rhsvec=reshape([0*xi;rhs(xi)],[],1);
odesolvec=[Ebc;Eode(3:end,:)]\[0;0;rhsvec(3:end)];
odesol=reshape(odesolvec,2,order+1);
plot(xi,odesol,'.-',xi,rhs(xi),'.-',lw{:})
xlabel('t');
legend({'x(t)','x''(t)','rhs(t)'});
set(gca,axprops{:});
title('Solution of the ODE x''''+cx''+kx=a*exp(-t^2/sigma^2)');

Published with MATLAB® R2018a