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.

Ex1_delay_cta( n , alpha, K )
function [tt,YY,T,Y] = Ex1_delay_cta( n , alpha, K )
% Wright equation x'(t)=-alpha x(t-1) [1+x(t)],
%                     x(t)=x0 for t in [-1,0]
% use: [tt,YY,T,Y]=Ex1_delay_cta(64,3,50);
%
tic;kind=2;dom=[-1,0];te=pd(n,dom,kind);Te=cpv(n,te,dom);
Dcal=Te*deriv(n,dom)/Te;Dcal(n,:)=0;
yinit=0.1*te.^3+te-0.5;
options = odeset('RelTol',1e-3,'AbsTol',1.e-6);
[tt,YY] = ode45(@cta,[0 K],yinit,options);toc;
%%%%%%%%%%%%%%%%%%%%%%%%%
tic;t(:,1)=pd(n,dom,kind);
y_val(:,1)=0.1*t(:,1).^3+t(:,1)-0.5;y_coef(:,1)=x2t(y_val(:,1),kind);
for k=2:K+1
    dom=[dom(1)+1,dom(2)+1];
    D=deriv(n,dom);
    F=fact(alpha*y_coef(:,k-1),n);
    A=D+F;b=-alpha*y_coef(:,k-1);
    A(n,:)=cpv(n,dom(1),dom);b(n)=y_val(n,k-1);
    y_coef(:,k)=A\b;t(:,k)=pd(n,dom,kind);
    y_val(:,k)=t2x(y_coef(:,k),kind);
end
T=reshape(t,n*(K+1),1);Y=reshape(y_val,n*(K+1),1);toc;
%%%%%%%%%%%%%%%%%%%%%%%%%
plot(T,Y,'-b');xlabel('t');ylabel('y(t)');grid;hold on;
plot(tt,YY(:,n),'.r');hold off;tt(end),T(end)
solcta=YY(end,n),solchp=Y(end)
function dy = cta(~,y)
    dy = Dcal*y;    % a column vector
    dy(n) = -alpha*y(1)*(1+y(n));
end
end

Contact us