Generalized Forced Van Der Pol Oscillator Phase Plot

29 次查看（过去 30 天）
Thomas Ward
Thomas Ward 2020-4-16
评论： george korris 2021-4-15
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

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2020-4-17
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.
.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by