clear all; global e f; % e = van der pol parameter % f = forcing parameter e=0.1; f=0; % initial conditions x0=2; y0=0; % time step dt=.1; Tf=100; % if you set this to =1 then the % program will also plot the solution % from the normal perturbation method (which resonates) plot_pert=0; % 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 circle with radius the initial conditions th = 0:pi/50:2*pi; xunit = x0 * cos(th); yunit = x0 * sin(th); r0 = sqrt(x0^2+y0^2); Xperturb = x0.* cos(T) + (3/8*x0^2*e).*(-T.*sin(T)); Xtwotimes = x0.* cos((1+3/8*x0^2*e)*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] = DUF(x,y,t); [k2,m2] = DUF(x+k1*dt/2,y+m1*dt/2,t+dt/2); [k3,m3] = DUF(x+k2*dt/2,y+m2*dt/2,t+dt/2); [k4,m4] = DUF(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,10)==1 % plot X(t) figure(1),clf; 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); if plot_pert==1 plot(T(1:it),Xperturb(1:it),'--k','LineWidth',2); end hold off; xlabel('\$t\$','FontSize',26,'Interpreter','latex'); ylabel('\$x\$','FontSize',26,'Interpreter','latex'); if plot_pert==1 ylim([-10 10]) else ylim([-5 5]) end if T(it)<80 xlim([0 100]) if plot_pert==1 xlim([0 Tf]) end else xlim([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',26,'Interpreter','latex'); ylabel('\$dx/dt\$','FontSize',26,'Interpreter','latex'); ylim([-5 5]) if T(it)<80 xlim([0 100]) if plot_pert==1 xlim([0 Tf]) end else xlim([round(Tf/4) Tf]) end % 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',26,'Interpreter','latex'); ylabel('\$dx/dt\$','FontSize',26,'Interpreter','latex'); r=max([max(X),max(Y),x0,y0,4]); axis([-r r -r r]); axis square drawnow; end end