| function revision_Hosein_paper |
| %% Inputs |
| tspan=0,10; % time span |
| x0=1; % initial condition |
| Delta_t=10; % Delta time (equal to tspan for no multistage) |
| points=100; % number of points to evaluate each Delta time |
| %% Application of the EMHPM algorithm to the proposed solutions |
| t1, z1=Semhpm(@he_stable,tspan,x0,Delta_t,points); % EMHPM by He's solution |
| t2, z2=Semhpm(@hn_unstable,tspan,x0,Delta_t,points); % EMHPM by Hosein Nia's solution |
| t3, z3=Semhpm(@od_taylor,tspan,x0,Delta_t,points); % EMHPM by Odibat's solution |
| t4, z4=ode45(@numeric,tspan,x0); % numerical ode45 solution |
| %% Display computer plots |
| plot(t1,z1,'- -',t2,z2,':',t3,z3,'-.',t4,z4,'-'); |
| legend('He - ∖ity_2','Hosein Nia et al. - ∖ity_2','Odibat - ∖ity_4', 'Numerical - ode45'); |
| ylabel('Displacement, ∖itx');xlabel(Time, t'); axis(0 10 0 2) |
| %% Solution Definitions |
| function dydt=numeric(t,y) % numerical solution |
| dydt(1,1)=y(1)-exp(t).*y(1)∧2; |
| function y=he_stable(t,x0,T) % homotopy/multistage unstable, He's solution |
| y0=x0*exp(t); |
| y1=.5*x0∧2*(exp(t)-exp(3*t)); |
| y2=.25*x0∧3*(exp(t)-2*exp(3*t)+exp(5*t)); |
| y=y0+y1+y2; |
| function y=hn_unstable(t,x0,T) % homotopy/multistage stable, Hosein Nia's solution |
| y0=x0*exp(-t); |
| y1=(2*x0-x0∧2)*t.*exp(-t); |
| y2=(2*x0-3*x0∧2+x0∧3).*t.∧2.*exp(-t); |
| y=y0+y1+y2; |
| function y=od_taylor(t,c,T) % homotopy/multistage Taylor Expansion, Odibat's solution |
| y0=c; |
| y1=(c-c.∧2)*t; |
| y2=(c.∧3-2*c.∧2+c/2)*t.∧2; |
| y3=(-c.∧4+3*c.∧3-13*c.∧2/6+c/6)*t.∧3; |
| y4=(c.∧5-4*c.∧4+29*c.∧3/6-5*c.∧2/3+c/24)*t.∧4; |
| y=y0+y1+y2+y3+y4; |
| %% Enhanced Multistage Homotopy Perturbation Method |
| function t,z=Semhpm(nde,tspan,z0,Deltat,pnts) |
| tini=tspan(1); tfin=tspan(end); tstart=tini; |
| tini=tini-tstart; tfin=tfin-tstart; |
| if tini == tfin |
| error('The last entry in tspan must be different from the first entry.'); |
| 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('The entries in tspan must strictly increase or decrease.'); |
| end |
| incT=Deltat/pnts; |
| z(1,:)=z0'; t(1)=tini; |
| iteT=2; |
| t(iteT)=tini+incT*tdir; % real time |
| while tdir*t(iteT)<tdir*(tfin+tdir*incT) |
| P_act=ceil(.99999*(t(iteT)-tini)*tdir/Deltat); % next subinterval |
| c=z((Deltat/incT)*(P_act-1)+1,:); % initial conditions |
| tsub(iteT-1)=t(iteT)-(P_act-1)*Deltat*tdir; % shifted time |
| temp=tsub(iteT-1); |
| z(iteT,:)=nde(temp,c,tstart+t(iteT)); % evaluate the solution |
| if tdir*t(iteT)>=tdir*tfin*0.99999 |
| break % last subinterval |
| else |
| iteT=iteT+1; |
| t(iteT)=tini+(iteT-1)*incT*tdir; % next subinterval |
| end |
| end |
| t=t'+tstart; |
