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.

[xx,solnum]=ibvp_ode_split_ex(n,dom,kind)
```function [xx,solnum]=ibvp_ode_split_ex(n,dom,kind)
% Example: 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-5
% From: Matlab documentation for bvp4c
% Exact solution: solex=cos(pi*x)+erf(x/sqrt(2*ep))/erf(1/sqrt(2*ep))
% call: [x,solnum]=ibvp_ode_split_ex(64,[-1 -0.05 0 0.05 1],2);
%
tic;
m=length(dom);ep=1.e-5;
A=spalloc(n*(m-1),n*(m-1),100*n*(m-1));b=zeros(n*(m-1),1);xx=[];
xs=pd(n,[-1 1],kind);Js=prim(n,[-1,1]);Xs=mult(n,[-1 1]);Ds=deriv(n,[-1,1]);
for k=1:m-1
pdom=[dom(k),dom(k+1)];l=(pdom(2)-pdom(1))/2;med=(pdom(2)+pdom(1))/2;
x=l*xs+med;J=l*Js;X=l*Xs+med*speye(n);xx=[xx,x];
myDE;
E=zeros(m-1);E(k,k)=1;AA(1:2,:)=0;A=A+kron(E,AA);
EE=zeros(m-1,1);EE(k)=1;bb(1:2)=0;b=b+kron(EE,bb);
end;
myBC; sol=A\b;% the solution in spectral form
toc;solnum=[];
for i=1:m-1
solnum = [solnum,t2x(sol((i-1)*n+1:i*n),kind)];
end
% the solution in physical form
x=reshape(xx,(m-1)*n,1);solnumv=reshape(solnum,(m-1)*n,1);
myOUT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myDE
f=-ep*pi^2*cos(pi*x)-pi*x.*sin(pi*x);
AA=ep*speye(n)+J*X-J^2;bb=J^2*x2t(f,kind);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myBC
T=cpv(n,[dom(1),dom(m)],[dom(1),dom(m)]);% T(1,:) for xc=d1, T(2,:) for xc=dm
Tl=T(1,:);Tr=T(2,:);
A(1,1:n)=T(1,:);A(2,(m-2)*n+1:(m-1)*n)=T(2,:);
b(1)=-2;b(2)=0;
for j=2:m-1
fl=(dom(j)-dom(j-1))/2;fr=(dom(j+1)-dom(j))/2;
A((j-1)*n+1,(j-2)*n+1:(j-1)*n)=Tr;
A((j-1)*n+1,(j-1)*n+1:j*n)=-Tl;b((j-1)*n+1)=0;
A((j-1)*n+2,(j-2)*n+1:(j-1)*n)=Tr*Ds/fl;
A((j-1)*n+2,(j-1)*n+1:j*n)=-Tl*Ds/fr;b((j-1)*n+2)=0;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myOUT
solex=cos(pi*x)+erf(x/sqrt(2*ep))/erf(1/sqrt(2*ep));
figure(1);
subplot(2,1,1);semilogy(x,abs(sol),'.');grid;
title('The absolute value of the coefficients of the solution');
subplot(2,1,2);plot(x,solnumv,dom,zeros(length(dom),1),'.r');
grid;xlabel('x');ylabel('y(x)');axis([-1 1 -2.1 2.1]);
title(['The solution (error = ',num2str(norm(solnumv-solex)),')']);
figure(2);
semilogy(x,abs(solnumv-solex));
title('The error');xlabel('x');ylabel('err');grid;
figure(3);spy(A);title('spy(A)');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end```