Code covered by the BSD License

# Chebpack

### Damian Trif (view profile)

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```