Approximate Solutions of Delay Differential Equations with Constant and Variable Coefficients by the Enhanced Multistage Homotopy Perturbation Method

Algorithm 1

MATLAB algorithm.

function dde23_ddeEM_mathieu_paper

disp(‘Mathieu equation solution’)

% Solution by EMHPM with zeroth and first order solution

eps=1; kappa=0.2; T=2pi; tau=T; b=−1; delta=3; % Mathieu Parameters

pntDelay=1; N=50; m_ord=5; dt=tau/(N−1); ktau=2; tspan=0,ktauT; % EMHPM Parameters

%% Solution

% Solution by dde23

tdde=linspace(0,ktautau,ktau(N−1)+1);

dde=@(t,y,z) mathieu_dde(t,y,z,kappa,delta,eps,b,tau);

te_dde=tic;sol = dde23(dde,tau,@history,tspan); toc(te_dde); xdde=deval(sol,tdde);

% Solution by zeroth EMHPM

dde_emhpm_fun=@(t,c0,tt,zm,zn,tau,N) mathieu_zeroth(t,c0,tt,zm,zn,tau,N,m_ord,

kappa,delta,eps,b,T);

tic;t0,z0=ddeEMHPM(dde_emhpm_fun,tspan,@history,tau,dt,pntDelay); toc

% Solution by First EMHPM

dde_emhpm_fun=@(t,c0,tt,zm,zn,tau,N) mathieu_first(t,c0,tt,zm,zn,tau,N,m_ord,kappa,delta,eps,b,T);

tic;t1,z1=ddeEMHPM(dde_emhpm_fun,tspan,@history,tau,dt,pntDelay); toc

% Plot results

ind0=1:2:ktau(N−1); ind1=2:2:ktau*(N−1);

Parent1=figure (1);

axes1 = axes(‘Parent’,Parent1,‘FontSize’,12,‘FontName’,‘Times New Roman’);

box(axes1,‘on’); hold(axes1,‘all’);

plot(tdde,xdde(1,:),‘Parent’,axes1,‘LineWidth’,2,‘Color’,0.502 0.502 0.502,‘DisplayName’,‘Numerical

dd23’);

plot(t0(ind0),z0(ind0,1),‘MarkerSize’,5,‘Marker’,‘o’,‘LineStyle’, ‘none’,‘DisplayName’,‘Zeroth

EMHPM’,‘Color’,0 0 0);

plot(t1(ind1),z1(ind1,1),‘MarkerSize’,7,‘Marker’,‘x’,‘LineStyle’,‘none’,‘DisplayName’,‘First

EMHPM’,‘Color’,0 0 0);

xlabel(‘∖itt’,‘FontSize’,12,‘FontName’,‘Times New Roman’);

ylabel(‘∖itx’,‘FontSize’,12,‘FontName’,‘Times New Roman’);

end

%% Mathieu definitions

function dydt = mathieu_dde(t,y,z,kapa,dlt,eps,b,T)

dydt = y(2)

    -kapay(2)−(dlt+epscos(2pit/T))y(1)+bz(1);

end

function Z = mathieu_zeroth(t,c0,tt,zm,zn,tau,N,m,kpa,dlt,eps,b,T)

Z=c0(1),c0(2); alf=dlt+epscos(2pi/Ttt);

for ik=1:m

    Z(ik+1,1)=Z(ik,2)t/ik;

    Z(ik+1,2)=−kpaZ(ik,2)−alfZ(ik,1);

    if ik==1, Z(ik+1,2)=Z(ik+1,2)+bzm(1); end

    Z(ik+1,2)=Z(ik+1,2)t/ik;

end

Z=sum(Z);

end

function Z = mathieu_first(t,c0,tt,zm,zn,tau,N,m,kpa,dlt,eps,b,T)

alf=dlt+epscos(2pi/Ttt); Z=c0(1),c0(2); Z_=0,0;

for ik=1:m

    Z(ik+1,:)=Z(ik,2)t/ik, −kpaZ(ik,2)−alfZ(ik,1);

    Z_(ik+1,:)=Z_(ik,2)t/ik, −kpaZ_(ik,2)−alfZ_(ik,1);

    if ik==1,

        Z(ik+1,2)=Z(ik+1,2)+bzm(1);

        Z_(ik+1,2)=Z_(ik+1,2)+b(N−1)/tau(zn(1)−zm(1))t;

    end

    Z(ik+1,2)=Z(ik+1,2)t/ik;

    Z_(ik+1,2)=Z_(ik+1,2)t/(ik+1);

end

Z=sum(Z)+sum(Z_);

end

function out=history(t)

out=1E−3+0t 0+0t;

end

%% EMHPM algorithm for ODE solutions

function t,z= odeEMHPM(nde,tspan,z0,Deltat,pnts)

% ode solver by Enhanced Multistage Homotopy Perturbation Method

tini=tspan(1); tfin=tspan(end); tstart=tini;

tini=tini−tstart; tfin=tfin−tstart; % shifted time set to zero

% Handle errors

if tini == tfin

    error(‘The ending and starting time values must be different.’);

elseif abs(tini)> abs(tfin)

    tspan=flipud(fliplr(tspan)); tini=tspan(1); tfin=tspan(end);

end

tdir=sign(tfin−tini);

if any(tdir*diff(tspan) <= 0)

    error(‘tspan entries must be strictly sorted.’);

end

incT=Deltat/pnts;

if incT<=0

    error (‘Increasing must be greater than zero.’)

end

z(1,:)=z0′; t(1)=tini; iteT=2; % set initial values

t(iteT)=tini+incTtdir;

while tdirt(iteT)<tdir*(tfin+tdirincT)

    P_act=ceil(.99999(t(iteT)−tini)tdir/Deltat); % count sub-intervals

    c=z((Deltat/incT)(P_act−1)+1,:); % set the corresponding initial condition

    tsub(iteT−1)=t(iteT)-(P_act−1)Deltattdir; % evaluate the solution at the shifted time

    temp=tsub(iteT−1);

    z(iteT,:)=nde(temp,c,tstart+t(iteT));

    if tdirt(iteT)>=tdirtfin0.99999 % repeat for the next sub-interval

        break

    else

        iteT=iteT+1; t(iteT)=tini+(iteT−1)incTtdir;

    end

end

t=t′+tstart;

end

%% EMHPM algorithm for DDE solutions

function t z=ddeEMHPM(dde,tspan,history,tau,Deltat,pntDelta)

% dde solver by Enhanced Multistage Homotopy Perturbation Method

tini=tspan(1); tfin=tspan(end); % span where the solution is founded

if tau<Deltat

    error(‘Subtinterval must not be greater than tau’);

end

if mod(tau,Deltat)~=0

    pastDlt=Deltat;

    Deltat=tau/round(tau/Deltat);

    warning(‘Subtinterval was modified from %0.5g to %0.5g.’,pastDlt,Deltat);

end

pnt=round(tau/Deltat); % samples−1 −tau 0

incT=Deltat/pntDelta; % step for set resolution

t=−tau+tini+(0:pntpntDelta)’incT; % evaluation of the initial solution −tau 0

z=history(t); % initial behavior in the inverval −tau 0

c=z(end,:); % initial condition

m=pntpntDelta; % samples between tau

itDlt=0;

while (t(end)+eps)<tfin;

    tspan=tini+itDltDeltat;tini+(itDlt+1)Deltat; % preparing the span for next Deltat

    zm=z(length(z)−m,:)’; zn=z(1+length(z)−m,:)’; % previous tau solution

    emhpm_fun=@(time,c,T) dde(time,c,T,zm,zn,tau,m); % application of the odeEMHPM

    t_aux z_aux=odeEMHPM(emhpm_fun,tspan,z(end,:),Deltat,pntDelta);

    z=z(1:end−1,:);z_aux; t=t(1:end−1,:);t_aux; % joining the solutions

    itDlt=itDlt+1;

end

z=z(pntpntDelta+1:end,:);t=t(pntpntDelta+1:end,:);

end