### 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 diﬀerent values of $$\mu$$ to compare eﬀect 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 ﬁrst converted to 2 ﬁrst 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 ﬁrst 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 diﬀerential equation for diﬀerent 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)];