4.12 Solve 2nd order ODE (Van Der Pol) and generate phase plot

Problem: Solve

\[ y^{\prime \prime }\left ( t\right ) -\mu \left ( 1-y^{2}\right ) y^{\prime }\left ( t\right ) +y\left ( t\right ) =0 \]

for different values of \(\mu \) to compare effect of changing \(\mu \) on the solution trajectories. The initial conditions are

\begin{align*} y\left ( 0\right ) & =0.1\\ y^{\prime }\left ( t\right ) & =0 \end{align*}

In both Mathematica and Matlab, numerical ODE solver was used.

For Matlab, The 2nd order ODE is first converted to 2 first order ODE’s, and then solve the coupled ODE’s using ode45. In Mathematica NDSolve was used. This does not require the conversion.

Starting by writing the 2nd order ODE as 2 first order ODE’s

\[\left . {\begin {array}[c]{l}x_{1}=y\\ x_{2}=y^{\prime }\end {array}} \right \} \text {derivatives} \Rightarrow \left . {\begin {array} [c]{l}x_{1}^{^{\prime }}=y^{\prime }\\ x_{2}^{^{\prime }}=y^{^{\prime \prime }}\ \end {array}} \right \} \Rightarrow \left . {\begin {array} [c]{l}x_{1}^{^{\prime }}=x_{2}\\ x_{2}^{^{\prime }}=\mu \left ( 1-x_{1}^{2}\right ) x_{2}+x_{1}\end {array}} \right \} \]

Below is the solution of the differential equation for different values of \(\mu =1\)

Mathematica

process[mu_]:=Module[ 
      {sol,p1,p2,p3,data,system,z,x1, 
       x2,t,ode1,ode2,timeScale}, 
ode1 =x1'[t]==x2[t]; 
ode2 =x2'[t]==z(1-x1[t]^2)x2[t]-x1[t]; 
timeScale ={t,0,80}; 
sol=First@NDSolve[ 
  {ode1,ode2/.z->mu,x1[0]==0.01,x2[0]==0}, 
  {x1[t],x2[t]}, 
   timeScale,Method->{"Extrapolation", 
   Method->"LinearlyImplicitEuler"}]; 
sol = {x1[t],x2[t]}/.sol; 
p1 = Plot[Evaluate[sol[[1]]], 
      Evaluate[timeScale], 
      Frame->True, 
      FrameLabel->{"time","y(t)"}, 
      PlotLabel->"\[Mu]="<>ToString[mu]]; 
p2 = Plot[Evaluate[sol[[2]]], 
          Evaluate[timeScale], 
      Frame->True, 
      FrameLabel->{"time","y'(t)"}, 
      PlotLabel->"\[Mu]="<>ToString[mu], 
      PlotStyle->RGBColor[1,0,0]]; 
data = Table[{evaluate[sol[[1]]], 
        Evaluate[sol[[2]]]}, 
        Evaluate[timeScale]]; 
p3=ListPlot[data,Joined->True,Frame->True, 
     PlotLabel->"\[Mu]="<>ToString[mu], 
     FrameLabel->{"y(t)","y'(t)"}, 
     PlotRange->All]; 
{p1,p2,p3} 
]; 
mu={0.1,.2,1,3,5,10}; 
r = process/@mu; Grid[r]
 

Matlab

function e71 
mu=[0.1,0.2,1,3,5,10]; 
 
for n=1:length(mu) 
    process(mu(n),n-1); 
end 
 
function process(mu,n) 
timeScale=[0 80]; 
ic=[0.1;0]; 
 
[t,x]=ode45(@ode,timeScale,ic,[],mu); 
subplot(6,3,(3*n)+1); 
plot(t,x(:,1)); 
title(sprintf('mu=%1.2f',mu)); 
xlabel('time'); ylabel('y(t)'); 
 
subplot(6,3,(3*n)+2); 
plot(t,x(:,2),'r'); 
title(sprintf('mu=%1.2f',mu)); 
xlabel('time'); ylabel('y''(t)'); 
 
subplot(6,3,(3*n)+3); 
plot(x(:,2),x(:,1)); 
title(sprintf('mu=%1.2f',mu)); 
xlabel('y(t)'); ylabel('y''(t)'); 
 
 
function xdot=ode(t,x,mu) 
xdot=[x(2) ; mu*(1-x(1)^2)*x(2)-x(1)];