Draw the vector field and eigenvectors in the phase portrait for Van der Pol ODE
24 ビュー (過去 30 日間)
古いコメントを表示
I have the follwing system which represent the Van der Pol oscillator with the inital condition and parameters are given. I draw the phase porrait using plot and ode45 but dont know how to draw the vector field and the eigenvectors with direction on them.
%function to solve the system with the time dependent term zero
function [dxdt] = vdp1(t,x,lambda,gamma,omega)
dxdt=zeros(2,1);
dxdt(1)=x(2);
dxdt(2)=lambda.*(1-x(1)^2)*x(2)-x(1)+gamma.*sin(omega*t);
end
%function to solve the system with the time dependent not zero
function [dxdt] = myode(t,x,gt,g,lambda,gamma,omega)
g=interp1(gt,g,t);
dxdt=zeros(2,1);
dxdt(1)=x(2);
dxdt(2)=lambda.*(1-x(1)^2)*x(2)-x(1)+g;
end
%script
lambda=[0.01 0.1 1 10 100] ;
gamma=[0 0.25];
omega=[0 1.04 1.1];
x0=[1 0];
x01=[3 0];
tspan=[0 500];
tspan1=[0 100];
%Numerical solution for the first initial value
[t,x]=ode45(@(t,x) vdp1(t,x,lambda(1),gamma(1),omega(1)),tspan,x0);
%Numerical solution for the second initial value
[t1,x1]=ode45(@(t,x) vdp1(t,x,lambda(1),gamma(1),omega(1)),tspan,x01);
%plotting x1,x2 aginst t
figure(3)
plot(x(:,1),x(:,2),'g-.')
hold on;
plot(x1(:,1),x1(:,2),'r-.')
xlabel('x1');
ylabel('x2')
legend('Solution first initial condition','Solution with the second initial condition')
title('phase portrait with t=[0 500] ,gamma=0,omega=0,lambda=0.01')
[tt1,xx1]=ode45(@(t,x) vdp1(t,x,lambda(2),gamma(1),omega(1)),tspan1,x0);
[tt2,xx2]=ode45(@(t,x) vdp1(t,x,lambda(3),gamma(1),omega(1)),tspan1,x0);
[tt3,xx3]=ode45(@(t,x) vdp1(t,x,lambda(4),gamma(1),omega(1)),tspan1,x0);
[tt4,xx4]=ode45(@(t,x) vdp1(t,x,lambda(5),gamma(1),omega(1)),tspan1,x0);
[tt11,xx11]=ode45(@(t,x) vdp1(t,x,lambda(2),gamma(1),omega(1)),tspan1,x01);
[tt22,xx22]=ode45(@(t,x) vdp1(t,x,lambda(3),gamma(1),omega(1)),tspan1,x01);
[tt33,xx33]=ode45(@(t,x) vdp1(t,x,lambda(4),gamma(1),omega(1)),tspan1,x01);
[tt44,xx44]=ode45(@(t,x) vdp1(t,x,lambda(5),gamma(1),omega(1)),tspan1,x01);
figure(6)
plot(xx1(:,1),xx1(:,2),'g-.')
hold on;
plot(xx11(:,1),xx11(:,2),'r-.')
xlabel('x1');
ylabel('x2')
legend('Solution first initial condition','Solution with the second initial condition')
title('phase portrait with t=[0 100] ,gamma=0,omega=0,lambda=0.1')
figure(8)
plot(xx2(:,1),xx2(:,2),'g-.')
hold on;
plot(xx22(:,1),xx22(:,2),'r-.')
xlabel('x1');
ylabel('x2')
legend('Solution first initial condition','Solution with the second initial condition')
title('phase portrait with t=[0 100] ,gamma=0,omega=0,lambda=1')
figure(10)
plot(xx3(:,1),xx3(:,2),'g-.')
hold on;
plot(xx33(:,1),xx33(:,2),'r-.')
xlabel('x1');
ylabel('x2')
legend('Solution first initial condition','Solution with the second initial condition')
title('phase portrait with t=[0 100] ,gamma=0,omega=0,lambda=10')
figure(12)
plot(xx4(:,1),xx4(:,2),'g-.')
hold on;
plot(xx44(:,1),xx44(:,2),'r-.')
xlabel('x1');
ylabel('x2')
legend('Solution first initial condition','Solution with the second initial condition')
title('phase portrait with t=[0 100] ,gamma=0,omega=0,lambda=100')
gt=[0 500];
g=gamma(2).*sin(omega(2).*gt);
g1=gamma(2).*sin(omega(3).*gt);
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t2,x2]=ode45(@(t,x) myode(t,x,gt,g,lambda(1),gamma(2),omega(2)),tspan,x0,opts);
[t22,x22]=ode45(@(t,x) myode(t,x,gt,g,lambda(1),gamma(2),omega(2)),tspan,x01,opts);
[t3,x3]=ode45(@(t,x) myode(t,x,gt,g,lambda(1),gamma(2),omega(3)),tspan,x0,opts);
[t33,x33]=ode45(@ (t,x) myode(t,x,gt,g1,lambda(1),gamma(2),omega(3)),tspan,x01,opts);
figure(14)
plot(x2(:,1),x2(:,2),'g-.')
hold on;
plot(x22(:,1),x22(:,2),'r-.')
xlabel('x1');
ylabel('x2')
legend('Solution first initial condition','Solution with the second initial condition')
title('phase portrait with t=[0 500] ,gamma=0.25,omega=1.04,lambda=0.01')
figure(16)
plot(x3(:,1),x3(:,2),'g-.')
hold on;
plot(x33(:,1),x33(:,2),'r-.')
xlabel('x1');
ylabel('x2')
legend('Solution first initial condition','Solution with the second initial condition')
title('phase portrait with t=[0 500] ,gamma=0.25,omega=1.1,lambda=0.01')
採用された回答
Agnish Dutta
2019 年 4 月 8 日
If you can calculate the vector field values at every point, then the resulting data can be plotted using the “quiver” function, the details of which are in the following document:
Quiver function - https://www.mathworks.com/help/matlab/ref/quiver.html
I also found a few useful resources on the internet pertinent to what you are trying to do. Refer to the “computing the vector field” section of the following website:
I believe that the following MATLAB answers page has an accepted answer relevant to your question:
0 件のコメント
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Pulsed Waveforms についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!