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?
Neeta Dsouza answered .
2025-11-20
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