Question 23.6: The van der Pol equation is a model of an electronic circuit...
The van der Pol equation is a model of an electronic circuit that arose back in the days of vacuum tubes,
\frac{d^2 y_1}{d t^2}-\mu\left(1-y_1^2\right) \frac{d y_1}{d t}+y_1=0 (E23.6.1)
The solution to this equation becomes progressively stiffer as μ gets large. Given the initial conditions, y_1(0) = dy_1/dt = 1, use MATLAB to solve the following two cases: (a) for μ = 1, use ode45 to solve from t = 0 to 20; and (b) for μ = 1000, use ode23s to solve from t = 0 to 6000.
(a) The first step is to convert the second-order ODE into a pair of first-order ODEs by defining
\frac{d y_1}{d t}=y_2
Using this equation, Eq. (E23.6.1) can be written as
\frac{d y_2}{d t}=\mu\left(1-y_1^2\right) y_2-y_1=0
An M-file can now be created to hold this pair of differential equations:
function yp = vanderpol(t,y,mu)
yp = [y(2);mu*(1−y(1)^2)*y(2)−y(1)];
Notice how the value of μ is passed as a parameter. As in Example 23.1, ode45 can be invoked and the results plotted:
>> [t,y] = ode45(@ vanderpol,[0 20],[1 1],[],1);
>> plot(t,y(:,1),'-',t,y(:,2),'--')
>> legend('y1','y2');
Observe that because we are not specifying any options, we must use open brackets [] as a place holder. The smooth nature of the plot (Fig. 23.10a) suggests that the van der Pol equation with μ = 1 is not a stiff system.
(b) If a standard solver like ode45 is used for the stiff case ( μ = 1000), it will fail miserably (try it, if you like). However, ode23s does an efficient job:
>> [t,y] = ode23s(@vanderpol,[0 6000],[1 1],[],1000);
>> plot(t,y(:,1))
We have only displayed the y_1 component because the result for y_2 has a much larger scale. Notice how this solution (Fig. 23.10b) has much sharper edges than is the case in Fig. 23.10a. This is a visual manifestation of the “stiffness” of the solution.
