clear all; global e f; % e = van der pol parameter % f = forcing parameter e=.1; f=0; % initial conditions x0=4; y0=0; % time step dt=.1; Tf=200; % time vector T = [0:dt:Tf]; X=0*T; Y=0*T; % inserting initial conditions in beginning % of X and Y vectors X(1)=x0; Y(1)=y0; % plot the limit cycle th = 0:pi/50:2*pi; xunit = 2 * cos(th); yunit = 2 * sin(th); r0 = sqrt(x0^2+y0^2); %the analytic solution from the two-time perturbation method Xtwotimes = ( 2./(sqrt(1 - (1-4/r0^2)*exp(-e*T)) )).* cos(T); figure(1),clf; % do the Runge-Kutta integration for it=2:length(T); x=X(it-1); y=Y(it-1); t=T(it-1); [k1,m1] = VP(x,y,t); [k2,m2] = VP(x+k1*dt/2,y+m1*dt/2,t+dt/2); [k3,m3] = VP(x+k2*dt/2,y+m2*dt/2,t+dt/2); [k4,m4] = VP(x+k3*dt,y+m3*dt,t+dt); X(it)=x+(k1+2*k2+2*k3+k4)*dt/6; Y(it)=y+(m1+2*m2+2*m3+m4)*dt/6; if rem(it,round(2/dt))==1 % plot X(t) figure(1); subplot(2,4,[1:2]) plot(T(1:it),X(1:it),'r','LineWidth',2); hold on; plot(T(1:it),Xtwotimes(1:it),'--b','LineWidth',2); hold off; xlabel('\$t\$','FontSize',18,'Interpreter','latex'); ylabel('\$x\$','FontSize',18,'Interpreter','latex'); title('red: RK4, blue: two-time analytic solution','FontSize',18,'Interpreter','latex'); ylim([-3 3]) ylim([-10 10]) if T(it)<80 xlim([0 100]) else xlim([0*round(Tf/4) Tf]) end % plot Y(t) subplot(2,4,[5:6]) plot(T(1:it),Y(1:it),'r','LineWidth',2) xlabel('\$t\$','FontSize',18,'Interpreter','latex'); ylabel('\$dx/dt\$','FontSize',18,'Interpreter','latex'); ylim([-5 5]) xlim([0 Tf]) % plot phase space subplot(2,4,[3,4,7,8]) plot(X(1:it),Y(1:it),'r','LineWidth',2); hold on; plot(xunit,yunit,'--m','LineWidth',2); plot(X(it),Y(it),'*r','MarkerSize',10,'LineWidth',2); hold off; xlabel('\$x\$','FontSize',18,'Interpreter','latex'); ylabel('\$dx/dt\$','FontSize',18,'Interpreter','latex'); r=1.1*max([max(X),max(Y),x0,y0,4]); axis([-r r -r r]); axis square % drawnow; end end