Code covered by the BSD License  

Highlights from
Chebpack

image thumbnail

Chebpack

by

 

15 Jul 2011 (Updated )

The MATLAB package Chebpack solves specific problems for differential or integral equations.

test2
function test2
% ep*y''+ x*y'= -ep*pi^2*cos(pi*x)-pi*x*sin(pi*x), x in [-1,1]
% y(-1)=-2, y(1)=0, ep=1.e-4
% exact solution: cos(pi*x)+erf(x/sqrt(2*ep))/erf(1/sqrt(2*ep))
%

display('Chebpack');
n=1024;dom=[-1,1];kind=2;
tic;x=pd(n,dom,kind);X=mult(n,dom);J=prim(n,dom);ep=1.e-4;
myDE;myBC;sol=A\b;toc;solnum = t2x(sol,kind);
solex=cos(pi*x)+erf(x/sqrt(2*ep))/erf(1/sqrt(2*ep));norm(solnum-solex)

display('bvp4c');

tic;sol =shockbvp;toc;
norm(sol.y(1,:)-cos(pi*sol.x)-erf(sol.x/sqrt(2*e))/erf(1/sqrt(2*e)))

display('bvp5c');
tic;sol =shockbvp(@bvp5c);toc;
norm(sol.y(1,:)-cos(pi*sol.x)-erf(sol.x/sqrt(2*e))/erf(1/sqrt(2*e)))

display('dms');
tic;N=1024;ep=1.e-4;g=[1 0 0;1 0 -2];[x,D2t,D1t,phip,phim]=cheb2bc(N,g);
f=-ep*pi^2*cos(pi*x)-pi*x.*sin(pi*x);D=ep*D2t+diag(x)*D1t;
p=ep*phip(:,2)+x.*phip(:,1);
m=ep*phim(:,2)+x.*phim(:,1);
u=D\(f-p-m);toc;solex=cos(pi*x)+erf(x/sqrt(2*ep))/erf(1/sqrt(2*ep));
norm(u-solex)

display('Chebfun');
tic;splitting on;dom = [-1, 1];
N = chebop(@(x,y) 0.0001.*diff(y,2)+x.*diff(y)+...
    0.0001.*pi.^2.*cos(pi.*x)+pi.*x.*sin(pi.*x),dom);
rhs = 0;N.lbc = @(y) y+2;N.rbc = @(y) y;
u = solvebvp(N,rhs);toc;
solex=chebfun('cos(pi*x)+erf(x/sqrt(2*1.e-4))/erf(1/sqrt(2*1.e-4))','splitting','on');
norm(u-solex)


function sol=shockbvp(solver)
if nargin < 1
  solver = 'bvp4c';
end
bvpsolver = fcnchk(solver);
options = bvpset('FJacobian',@shockJac,'BCJacobian',@shockBCJac,'Vectorized','on');
sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);
e = 0.1;
for i=2:4
  e = e/10;   
  sol = bvpsolver(@shockODE,@shockBC,sol,options);
end 
  function dydx = shockODE(x,y)
    pix = pi*x;
    dydx = [                 y(2,:)
             -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix) ];   
  end
  function res = shockBC(ya,yb)
    res = [ ya(1)+2
            yb(1)  ];
  end
  function jac = shockJac(x,y)
    jac = [ 0   1
            0 -x/e ];
  end
  function [dBCdya,dBCdyb] = shockBCJac(ya,yb)
    dBCdya = [ 1 0
               0 0 ];

    dBCdyb = [ 0 0
               1 0 ];
  end
end 

    function myDE
        f=-ep*pi^2*cos(pi*x)-pi*x.*sin(pi*x);
        A=ep*speye(n)+J*X-J^2;b=J^2*x2t(f,kind);
    end
    function myBC
        T=cpv(n,dom,dom);
        A(1,:)=T(1,:);b(1)=-2;A(2,:)=T(2,:);b(2)=0;
    end
end

Contact us