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.

B=Ex4_delay_main
function B=Ex4_delay_main
% use: B=Ex4_delay_main;
n=96;
K=1;display(K);dom=[0 1];kind=2;t=pd(n,dom,kind);
Pinit=3.5;zinit=x2t(3*sin(2*pi*t),kind);Uinit=[zinit;Pinit;3.4];
[~,~,~,~,P,U] = Ex4_delay_bif( n , Uinit);pause(1);
B(1:n,K)=U(1:n);B(n+1,K)=P;B(n+2,K)=lambda;
figure(2);plot(lambda,norm(B(1:n,K)),'.');grid;hold on;pause(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K=2;display(K);
Uinit=[B(1:n+1,K-1);3.41];
[~,~,~,~,P,U] = Ex4_delay_bif( n , Uinit);pause(1);
B(1:n,K)=U(1:n);B(n+1,K)=P;B(n+2,K)=lambda;
figure(2);plot(lambda,norm(B(1:n,K)),'.');grid;hold on;pause(1);alpha=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for K=3:225
    Uinit=(1+alpha)*B(:,K-1)-alpha*B(:,K-2);
    [~,~,~,~,P,U,nriter] = Ex4_delay_bif( n , Uinit);pause(1);
    B(1:n,K)=U(1:n);B(n+1,K)=P;B(n+2,K)=lambda;display([K,lambda]);
    save('BIF','B');
    figure(2);plot(B(n+2,K),norm(B(1:n,K)),'.');grid;hold on;pause(1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tt,tth,solnum,solnumh,P,U,nriter] = Ex4_delay_bif( n , Uinit)
dom=[0 1];domh=[-1,0];kind=2;%lambda=Uinit(n+2);
t=pd(n,dom,kind);th=pd(n,domh,kind);D=deriv(n,dom);
T=cpv(n,t,dom);%
options=optimset('Display','iter','TolFun',1.e-9,'TolX',1.e-9,...
    'MaxFunEvals',200000,'MaxIter',1000,'Algorithm','Levenberg-Marquardt');
[U,~,~,output]=fsolve(@fun_ex5,Uinit,options);nriter=output.iterations;
solnum(:,1)=t2x(U(1:n),kind);solnumh=solnum;P=U(n+1);lambda=U(n+2);
tt=P*t;tth=P*th;I3=find(tth>=-1.65);
figure(1);plot(tth(I3),solnumh(I3),'--',tt,solnum,'-');
xlabel('t');ylabel('y(t)');
grid;title(['\lambda= ',num2str(lambda)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=fun_ex5(U)
    F=zeros(2*n+1,1);z=U(1:n);P=U(n+1);lambda=U(n+2);
    tw1=t-1.65/P+(1-sign(t-1.65/P))/2;F1=cpv(n,tw1,dom)*z;
    tw2=t-0.35/P+(1-sign(t-0.35/P))/2;F2=cpv(n,tw2,dom)*z;
    F(1:n)=D*z+lambda*P*x2t((F1+F2).*(1+T*z),kind);
    if K==1, F(n+2)=U(n+2)-3.4;elseif K==2, F(n+2)=U(n+2)-3.41;
    else F(n+2)=(U-Uinit)'*(Uinit-B(:,K-1));end
    F(n+1)=T(1,:)*z-1;% phase condition
    F(n)=(T(1,:)-T(n,:))*z;% continuity condition
end
end
end

Contact us