% ### VDPode45EX.m ### 09.21.14 % Numerically integrate the van der Pol oscillator % m*x''- P.mu*(1-x^2)*x' + k*x = A*sin(wt) clear % ----------------------------------------------------- % User input (Note: All paramters are stored in a structure) P.y0(1) = 0.01; % initial position [m] P.y0(2) = 0.0; % initial velocity [m/s] P.mu= 0.5; % damping coefficient [kg/s] P.k= 1.0; % stiffness [N m] P.m= 1.0; % mass [kg] % sinusoidal driving term P.A= 0.0; % amplitude [N] (set to zero to turn off) fD= 1.05*sqrt(P.k/P.m)/(2*pi); % freq. (Hz) [expressed as fraction of resonant freq.] % Integration limits P.t0 = 0.0; % Start value P.tf = 100.0; % Finish value P.dt = 0.01; % time step % ---------------------------------------------------------------------- % +++ % spit back out some basic derived quantities P.wr= 2*pi*fD; % convert to angular freq. % +++ % use built-in ode45 to solve [t y] = ode45('VDPfunction', [P.t0:P.dt:P.tf],P.y0,[],P); % ------------------------------------------------------ % visualize figure(1); clf; plot(t,y(:,1)); hold on; grid on; xlabel('t [s]'); ylabel('x(t) [m]') % Phase plane figure(2); clf; plot(y(:,1), y(:,2)); hold on; grid on; xlabel('x [m]'); ylabel('dx/dt [m/s]')