Math 542: Numerical Solutions to Differential Equations

vdp_9stage
%
%
%
clear all

%
% Right-Hand-Side of the Van der Pol Equation
%
f = inline('[y(2) + mu*y(1) - mu*y(1)^3/3; -y(1)]','mu','t','y');

%mu      = input('mu = ');
mu      = 1e-4;
h       = 1/20;  % Step-size for RK-method
zN      = 0.1;     % The starting value on the z-axis
z0      = 0;
cyc     = 0;        % The counter for # of cycles
tol     = 0.01;     % Tolerance for return-map


% This looks for a limit cycle!
y_all = [];
while( abs(z0-zN)>tol )
  z0  = zN;
  cyc = cyc + 1;
  t   = 0;        % reset the clock for this iteration
  y   = [0; z0];  % set the initial point
  ctr = 0;
  go  = 1;
  % Step forward with explicit 9-stage, 7-order RK-method
  while( go == 1 )
    yn    = y(:,ctr+1);
    k1    = f(mu,t,        yn);
    k2    = f(mu,t+h/6,    yn+h*k1/6);
    k3    = f(mu,t+h/3,    yn+h*k2/3);
    k4    = f(mu,t+h/2,    yn+h*(k1/8+3*k3/8));
    k5    = f(mu,t+2*h/11, yn+h*(148*k1/1331+150*k3/1331-56*k4/1331));
    k6    = f(mu,t+2*h/3,  yn+h*(-404*k1/243-170*k3/27+4024*k4/1701+ ...
			      10648*k5/1701));
    k7    = f(mu,t+6*h/7,  yn+h*(2466*k1/2401+1242*k3/343-19176*k4/16807- ...
			      51909*k5/16807+1053*k6/2401));
    k8    = f(mu,t,        yn+h*(5*k1/154+96*k4/539-1815*k5/20384- ...
			      405*k6/2464+49*k7/1144));
    k9    = f(mu,t+h,      yn+h*(-113*k1/32-195*k3/22+32*k4/7+29403*k5/3584 ...
			      -729*k6/512+1029*k7/1408+21*k8/16));
    ctr   = ctr+1;
    ynext = yn + h*(32*k4/105 + 1771561*k5/6289920 + 243*k6/2560 + ...
		    16807*k7/74880 + 77*k8/1440 + 11*k9/270);
    if( isnan(ynext) )
      error('huh?')
    end
    y     = [y ynext];
    % Stop when we cross the z-axis from the 2nd to 1st quadrant
    if( norm(ynext)<1e-8 )
      go = 0;
    end
    if( yn(1)*ynext(1)<0 & yn(2)>0 & ynext(2)> 0 )
 	go = 0;
    end
  end
  % Save the data to file(s)...
  y = y';
  %  save('-ascii',sprintf('vdp02_%03d_z0_%.3f.dat',cyc,z0),'y');
  zN = y(length(y),2);
  %% Plot...
  plot(y(:,1),y(:,2))
  y_all = [y_all ; y];
end

plot(y_all(:,1),y_all(:,2))
grid on
title(sprintf('Limit Cycle, mu=%f',mu))
print('-depsc','-f1',sprintf('vdp_lc_mu_%f.eps',mu))
Main Content
