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.

ex4.m
function Falpha=ex4(alpha)
% from: Tadeusz Jankowski, Positive solutions for second order impulsive
% differential equations involving Stieltjes integral conditions, 
% Nonlinear Analysis 74 (2011) 37753785
% use: for alpha=0:0.1:40,Falpha=ex4(alpha);
%      plot(alpha,Falpha,'.k','MarkerSize',8);pause(0.1);
%      hold on;end;grid;hold off;
% or   Falpha=ex4(3.8674);% for example
% 
n=64;kind=2;Q1=-0.5;
doml=[0,1/2];tl=pd(n,doml,kind);cosl=cos(2*pi*tl);
Dl=deriv(n,doml);Jl=prim(n,doml);Tl=cpv(n,doml,doml);CL=fact(x2t(cosl,kind));
domr=[1/2,1];tr=pd(n,domr,kind);cosr=cos(2*pi*tr);
Dr=deriv(n,domr);Jr=prim(n,domr);Tr=cpv(n,domr,domr);CR=fact(x2t(cosr,kind));
A11=Dl^2;A12=zeros(n);A21=zeros(n);A22=Dr^2;
A11(n-1,:)=Tl(1,:);A11(n,:)=Tl(2,:);A12(n,:)=-Tr(1,:);
A21(n-1,:)=Tl(1,:)*Dl;A21(n,:)=-Tl(2,:)*Dl;
A22(n,:)=Tr(1,:)*Dr;
A=[A11,A12;A21,A22];b=zeros(2*n,1);
startl=0.5;startr=0.5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xoldl=zeros(n,1);xoldl(1)=2*startl;
xoldr=zeros(n,1);xoldr(1)=2*startr;
Xold=[xoldl;xoldr];% starting iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for iter=1:300
    xlold=Xold(1:n);
    b(1:n)=x2t(f(t2x(xlold,kind)),kind);b(n-1:n)=0;
    xrold=Xold(n+1:2*n);
    b(n+1:2*n)=x2t(f(t2x(xrold,kind)),kind);b(2*n-1)=alpha;b(2*n)=Q1;
    Xnew=A\b;err=norm(Xnew-Xold);%display(err);
    if err<1.e-9, break;end
    Xold=Xnew;
end  
Falpha=((Tl(2,:)-Tl(1,:))*Jl*CL)*Xnew(1:n)+...
    (Tr(2,:)+(Tr(2,:)-Tr(1,:))*Jr*CR)*Xnew(n+1:2*n);
% Uncomment for Falpha=ex4(...);
% xl=t2x(Xnew(1:n),kind);xr=t2x(Xnew(n+1:2*n),kind);t=[tl;tr];x=[xl;xr];
% figure(1);plot(t,x);xlabel('t');ylabel('x');grid;
% title(['x''(0)= ',num2str(alpha),', F= ',num2str(Falpha)]);
% test4(alpha);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test4(alpha)
options=odeset('RelTol',1.e-9,'AbsTol',1.e-9);
[~,yl]=ode45(@fun4,tl,[0,alpha],options);
[~,yr]=ode45(@fun4,tr,[yl(end,1),yl(end,2)-0.5],options);
figure(2);
semilogy([tl;tr],abs([xl-yl(:,1);xr-yr(:,1)]));grid;
xlabel('t');ylabel('Difference between Chebpack and ode45');
title(['x''(0)= ',num2str(alpha),', F= ',num2str(Falpha)]);
function dydt=fun4(t,y)
    dydt=zeros(2,1);
    dydt(1)=y(2);
    if y(1)>=0&&y(1)<=1;F=(y(1)-1).^2/10;
    elseif y(1)>1&&y(1)<=2;F=2*(y(1)-1).^2;
    elseif y(1)>2&&y(1)<=9/2;F=2*(y(1)+3)/5;
    else F=3;
    end
    dydt(2)=-30*pi^2/(1+pi^2)*F;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=f(y)
    F=zeros(size(y));
    I1= y>=0&y<=1;F(I1)=(y(I1)-1).^2/10;
    I2=y>1&y<=2;F(I2)=2*(y(I2)-1).^2;
    I3=y>2&y<=9/2;F(I3)=2*(y(I3)+3)/5;
    I4=y>9/2;F(I4)=3;
    F=-30*pi^2/(1+pi^2)*F;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Contact us