Generalized Forced Van Der Pol Oscillator Phase Plot

Hello,

I am trying to write a program to solve a forced van der pol oscillator and give its phase plot. I found a code on these forums, posted by @KSSV that solved an unforced van der pol oscillator and I just added a forcing term, and it worked. Here is the code

function VanderPol()

[t,y] = ode23(@vdp1,[0 20],[2; 0]);

plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE23');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

figure
plot(y(:,1),y(:,2))
title('Phase plane plot')
end
function dydt = vdp1(t,y)


dydt = [y(2); 3*(1-y(1)^2)*y(2)-y(1)+8*sin(4*t)];
end

I was able to go in and play with the variables to obtain different results, which was great, but ideally I would want to generalize the code such that the second line read

[t,y] = ode23(@vdp1,[0 20],[idis; ivel]);

And the line prior to end read

dydt=[y(2); mu*(1-y(1)^2)*y(2)-y(1)+A*sin(omega*t)];

And then be able to call VanderPol(idis, idel, mu, A, omega) and be able to solve that without having to go in and manually change it each time

Answer 1

[Star Strider]

Try this:

function [t,y] = VanderPol(idis, ivel, mu, A, omega)
[t,y] = ode23(@vdp1,[0 20],[2; 0]);
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE23');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')
figure
plot(y(:,1),y(:,2))
title('Phase plane plot')

function dydt = vdp1(t,y)
dydt = [y(2); mu*(1-y(1)^2)*y(2)-y(1)+A*sin(omega*t)];
end
end

Then to call it:

idis = 2;
ivel = 0;
A = 8;
omega = 4;
mu = 3;
[t,y] = VanderPol(idis, ivel, mu, A, omega);

figure
plot(t,y)
grid

I added a tweak so that you can get the integrated results from your funciton. It is otherwise what you wrote.

.

[Thomas Ward]

Thank you

[Star Strider]

As always, my pleasure!

[george korris]

hey guys is this code solving this equation: and what does: Solution of van der Pol Equation (\mu = 1) with ODE23' that the first figure has as title mean?