Change the period/frequency of a Van der Pol Oscillator

10 visualizzazioni (ultimi 30 giorni)
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?

Risposta accettata

David Goodmanson
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 Commenti
Alexis
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
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
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 Commenti
Alexis
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).

Accedi per commentare.

Community Treasure Hunt

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

Start Hunting!

Translated by