Change the period/frequency of a Van der Pol Oscillator

f = @(t,x) [x(2); -lambda*(x(2)^2 + x(1)^2 - a)*x(2) - w^2*x(1)]; (1)

where the y(t) term was taken away since it was set to 0 anyway. Rather than anonymous functions, there are two functions defined at the bottom of the script code below. The relevant line for the time derivatives is

dxy = [x(2); (-lam1*x(2)^2 -lam2*x(1)^2 +lam3*a)*x(2) - w^2*x(1)]

which is the same except there are three adjustable lambda values instead of just one. Suppose you start with (1) and want to speed up the waveform by a factor of b and change its size by a factor of A, all the while keeping the same shape. This can be done with

lam1 = lambda/(b*A^2);
lam2 = lambda*b/A^2;
lam3 = lambda*b;
and w --> w*b

If you speed up a waveform by a factor of b, then the acceleration goes up by a factor of b^2. The reason for A is that if you want to keep the acceleration the same as before, you can correct by a factor of A = 1/b^2. Or of course you can set the acceleration to any size, within reason.

The initial conditions should really be adjusted to get exact agreement for short times, but I left that alone because the solution settles down pretty quickly.

% define parameters
lambda = 3;
w = 0.8;                                                                    
a = 1;                                                                      
x0 = [.3; 0];
tspan = [0 50];

% reproduce original waveform
b = 1;  A = 1;
lam1 = lambda/(b*A^2);
lam2 = lambda*b/A^2;
lam3 = lambda*b;
[t1,x1] = ode45(@(t,x) fun(t,x,lam1,lam2,lam3,a,w*b),tspan,x0);
a1 = accel(x1,lam1,lam2,lam3,a,w*b);
figure(1)
plot(t1,a1)
grid on

% speed up waveform by a factor of 2, also increases acceleration by a factor of 4 
b = 2;  A = 1;
lam1 = lambda/(b*A^2);
lam2 = lambda*b/A^2;
lam3 = lambda*b;
[t2,x2] = ode45(@(t,x) fun(t,x,lam1,lam2,lam3,a,w*b),tspan,x0);
a2 = accel(x2,lam1,lam2,lam3,a,w*b);
figure(2)
plot(t2,a2)
grid on

% demo to show size change, no speedup
b = 1;  A = 3;
lam1 = lambda/(b*A^2);
lam2 = lambda*b/A^2;
lam3 = lambda*b;
[t3,x3] = ode45(@(t,x) fun(t,x,lam1,lam2,lam3,a,w*b),tspan,x0);
a3 = accel(x3,lam1,lam2,lam3,a,w*b);
figure(3)
plot(t3,a3)
grid on

% speed up waveform by factor of 2, maintain the original size of the acceleration
b = 2;  A = 1/b^2;
lam1 = lambda/(b*A^2);
lam2 = lambda*b/A^2;
lam3 = lambda*b;
[t4,x4] = ode45(@(t,x) fun(t,x,lam1,lam2,lam3,a,w*b),tspan,x0);
a4 = accel(x4,lam1,lam2,lam3,a,w*b);
figure(4)
plot(t4,a4)
grid on

return

function dxy = fun(t,x,lam1,lam2,lam3,a,w)  
  dxy = [x(2); (-lam1*x(2)^2 -lam2*x(1)^2 +lam3*a)*x(2) - w^2*x(1)];
end

function a = accel(x,lam1,lam2,lam3,a,w)  
  a = ((-lam1*x(:,2).^2 -lam2*x(:,1).^2 +lam3*a).*x(:,2) - w^2*x(:,1));
end