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 |
|