Change the period/frequency of a Van der Pol Oscillator

Hello everyone,

I'm having a trouble/issue changing the frequency of a van der Pol oscillator and getting an specific shape of the plot.

I'm using a pedestrian lateral forces model that walkers apply to bridges while walking and it uses a van der Pol oscillator as follows.

Setting the parameters (lambda = 3, a = 1, omega = 0.8, initial condition x0 = [3 -1]) I can get the acceleration of the system using the following code:

%% Init
clear variables
close all
clc
%% Code
% define parameters
lambda = 3;                                                                 % damping coefficient
% f = 0.8;       % [hz]                                                     % Frequency
% w = 2*pi*f;    % ~ 5[rad/sec]                                             % natural frequency
w = 0.8;                                                                    % natural frequency
a = 1;                                                                      % nonlinearity parameter
y = @(t) 0;
m = 90;

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

% initial conditions
x0 = [3;-1];

% time span
tspan = [0 100];

% solve using ode45
[t,x] = ode45(f, tspan, x0);

% compute acceleration
xpp = -lambda*(x(:,2).^2 + x(:,1).^2 - a).*x(:,2) - w^2.*x(:,1) + y(t);

% plot solution
figure;
plot(t, xpp(:,1), 'b', 'LineWidth', 1.5);
xlim([40 50])
xlabel('t');
ylabel('xpp');
title('Acceleration of a Van der Pol Oscillator');
grid on

The shape looks like to the way a pedestrian applies lateral forces while walking according to multiple researches, BUT, the waves are far from each other. So my problem is that I want the waves come closer each other because the pedestrian frequency is approx f = 0.8[hz] or w = 2*pi*f = 5 [rad/sec] and as soon as I change the frequency the shape changes. I've been trying for hours to find a combination of parameters that produces exactly the same shape but the waves come closer to each other. Here is what I get (not the same shape, but I have the frequency)

Any help or recomendation to find the parameters?

Answer 1

[David Goodmanson]

Hi Alexis,

If you want to preserve the shape you will have to have a larger set of parameters. I am taking as given your equation

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

[Alexis]

@David Goodmanson, that is a very elegant solution to my problem, and I think that is exactly what I was asking for. I'm going to recalibrate the parameters and write them in this topic before accepting the answer.

[Alexis]

I used the following:

lambda = 9;
a = 1;
x0 = [0.02,0.01];
f = 0.8;            % I changed w for f because it's actually a frequency f = 0.8[hz] 
b = 2*pi;           % So w = f*b = 2*pi*f
A = 1/b^2;

And I got this signal that is similar enough to the measured one.

Thanks you so much!

Answer 2

[Mathieu NOE]

hello

seems to me your code actually generates the correct frequency signal at 0.8 Hz , as soon as you stick with

f = 0.8; % [hz] % Frequency

w = 2*pi*f; % ~ 5[rad/sec] % natural frequency

and not

% w = 0.8; % natural frequency

minor tweaks / complements in your code

%% Init
clear variables
close all
clc
%% Code
% define parameters
lambda = 3;                                                                 % damping coefficient
f = 0.8;       % [hz]                                                     % Frequency
w = 2*pi*f;    % ~ 5[rad/sec]                                             % natural frequency
% w = 0.8;                                                                    % natural frequency
a = 1;                                                                      % nonlinearity parameter
y = @(t) 0;
% m = 90;       % what for ? 

% define equation
f = @(t,x) [x(2); -lambda*(x(2)^2 + x(1)^2 - a)*x(2) - w^2*x(1) + y(t)];
% initial conditions
x0 = [3;-1];
% time span
tspan = [0 20];
% solve using ode45
[t,x] = ode45(f, tspan, x0);
% compute acceleration
xpp = -lambda*(x(:,2).^2 + x(:,1).^2 - a).*x(:,2) - w^2.*x(:,1) + y(t);

% remove first seconds (transient)
ind = (t>10);
t = t(ind);
xpp = xpp(ind);

% find signal frequency by measuring "zero" (or any other value crossing
% time index

threshold = 0; % 

t0_pos1 = find_zc(t,xpp,threshold);
period = diff(t0_pos1); % delta time of crossing points
freq = 1./period; % signal frequency

figure(1)
subplot(2,1,1),plot(t,xpp,'b-',t0_pos1,threshold*ones(size(t0_pos1)),'.r','linewidth',2,'markersize',20);grid on
xlim([min(t) max(t)]);
legend('signal','signal positive slope crossing points');
xlabel('Time (s)');
ylabel('xpp');
title('Acceleration of a Van der Pol Oscillator');
subplot(2,1,2),plot(t0_pos1(2:end),freq,'b.-','linewidth',2,'markersize',12);grid on
xlim([min(t) max(t)]);
ylim([0 1]);
legend('signal rate (frequency)');
xlabel('Time (s)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)

    % positive slope "zero" crossing detection, using linear interpolation
    y = y - threshold;
    zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
    ix=zci(y);                      %find indices of + zero crossings of x
    ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
    Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Sam Chak]

The following is not a solution, but I think OP wants to preserve the shape of the wave and the wave should oscillates at a higher frequency. Something like the plot generated below (which I manipulated using time for illustration purposes only.)

%% Init
clear variables
close all
clc
%% Code
% define parameters
lambda = 3;                                                                 % damping coefficient
% f = 0.8;       % [hz]                                                     % Frequency
% w = 2*pi*f;    % ~ 5[rad/sec]                                             % natural frequency
w = 0.8;                                                                    % natural frequency
a = 1;                                                                      % nonlinearity parameter
y = @(t) 0;
m = 90;

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

% initial conditions
x0 = [3;-1];

% time span
tspan = [0 100];

% solve using ode45
[t,x] = ode45(f, tspan, x0);

% compute acceleration
xpp = -lambda*(x(:,2).^2 + x(:,1).^2 - a).*x(:,2) - w^2.*x(:,1) + y(t);

% plot solution
figure;
plot(t/4, xpp(:,1), 'b', 'LineWidth', 1.5);
xlim([10 20]), 
ylim([-2 2])
xlabel('t');
ylabel('xpp');
title('Acceleration of a Van der Pol Oscillator');
grid on

[Alexis]

Thank you both for trying to solve this.

@Mathieu NOE, yes when I change the frequency to w = 2*pi*f, I get the frequency that I want for my oscillator, but it's not the same shape. (the m I wrote is for the mass of the pedestrian that I will use for get the force Force = m*ddx)

@Sam Chak, that is the shape that I'm looking for, but I can't change the time because eventually I want to use this in a simulink MatlabFunction block. It works for visual purposes but would not work in my simulations.

As a reference I've this figure that I want to replicate (Bocian et al 2015)

[Mathieu NOE]

hello

why not digitize this picture and take that signal shape as "master curve"

simply define the time axis to make it the right frequency as you wish

that seems to be the simplest wat IMHO

[Alexis]

@Mathieu NOE, I didn't know that I could do that. But it is not what I am looking for because I want to add uncertainties in the model choosing random pedestrian parameters (the van der pol parameters) from distributions and this would only be a deterministic case (even if I make minor changes 200 times). Also I want to run it in simulink and I'm still not sure If fixed-step solver is time efficient for the multiple simulations (because you told me that I need to define the time axis "before" my simulation).

[Mathieu NOE]

ok I understand