# Change the period/frequency of a Van der Pol Oscillator

10 visualizzazioni (ultimi 30 giorni)
Alexis il 22 Mar 2023
Modificato: Alexis il 24 Mar 2023
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?
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Risposta accettata

David Goodmanson il 23 Mar 2023
Modificato: David Goodmanson il 23 Mar 2023
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
##### 2 CommentiMostra NessunoNascondi Nessuno
Alexis il 23 Mar 2023
@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 il 24 Mar 2023
Modificato: Alexis il 24 Mar 2023
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!

Accedi per commentare.

### Più risposte (1)

Mathieu NOE il 22 Mar 2023
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
##### 5 CommentiMostra 3 commenti meno recentiNascondi 3 commenti meno recenti
Alexis il 22 Mar 2023
Modificato: Alexis il 22 Mar 2023
@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 il 22 Mar 2023
ok I understand

Accedi per commentare.

### Categorie

Scopri di più su Sources in Help Center e File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by