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.

[t,y]=ex2(n)
function [t,y]=ex2(n)
% From C. Marriott and C. DeLisle, Effects of discontinuities in the 
% behavior of a delay differential equation, Physica D, 36 (1989) 198-206
% use: [t,y]=ex2(64);
%
tic;a=0.16;ep=0.02;u=0.5;xb=-0.427;tau=1;tol=1.e-9;kind=2;t=[];y=[];
grila=-12:12:120;
dom=[grila(1),grila(2)];id=2;
told=pd(n,dom,kind);t=[t,told];
yvold=0.6*ones(n,1);y=[y,yvold];% complete the history until grila(2)
while dom(2)<120
    dom=[grila(id),grila(id+1)];% the next interval
    pas;% calculate the next step and modify the history and the grid
end
toc;plot(t(:),y(:),'.k');xlabel('t');ylabel('y(t)');grid;
%y(end,end)
function pas
% calculate the solution on the current dom
tnew=pd(n,dom,kind);D=deriv(n,dom);T=cpv(n,dom(1),dom);
A=D+eye(n)/tau;A(n,:)=T;
hid=find(abs(grila-(dom(1)-12))<tol);% find the history
del=y(:,hid)-xb;
b=x2t(pi/tau*(a-u*sin(del).^2+ep*sign(sum(sign(del)))),kind);b(n)=y(end,end);
ycnew=A\b;yvnew=t2x(ycnew,kind);
% find the event y=xb on dom
S=sign((yvnew(1:n-1)-xb)).*sign((yvnew(2:n)-xb));ind=find(S==-1,1);% 
if isempty(ind),t=[t,tnew];y=[y,yvnew];id=id+1;
else z=chebroots(x2t(funex2(tnew),kind),dom,[tnew(ind)-tol,tnew(ind+1)+tol],tol);
     tpartl=pd(n,[dom(1),z],kind);tpartr=pd(n,[z,dom(2)],kind);
     t=[t,tpartl];y=[y,barycheb(tpartl,yvnew,tnew,kind)];
     t=[t,tpartr];y=[y,barycheb(tpartr,yvnew,tnew,kind)];
     grila=sort(union(grila,tpartl(end):12:120));id=id+2;
end
function f=funex2(x)
    f=barycheb(x,yvnew-xb,tnew,kind);
end
end
end

Contact us