home

PDF (letter size)

Matlab source code nma_project_511.m.txt

1 Introduction

Report for course EGME 511 (Advanced Mechanical Vibration). California state university, Fullerton, CA. Van der Pol differential equation is given by \[ x^{\prime \prime }\left ( t\right ) -c\left ( 1-x^{2}\left ( t\right ) \right ) x^{\prime }\left ( t\right ) +kx\left ( t\right ) =0 \] In this analysis, we will consider the case only for positive \(c,k\).  The above equation will be solved numerically using Matlab’s ODE45 for different initial conditions, and the phase portrait (velocity vs. displacement) is plotted to show the limit cycle for different initial conditions.

2 Matlab implementation

To use ODE45, one must first convert the above second order ODE to two ODE’s, each of which is first order. Letting \(x_{1}=x, x_{2}=x^{\prime }(t)\) results in

\[ \left . \begin {array} [c]{c}x_{1}=x\\ x_{2}=x^{\prime }\left ( t\right ) \end {array} \right \} \left . \begin {array} [c]{c}x_{1}^{\prime }=x^{\prime }\left ( t\right ) =x_{2}\\ x_{2}^{\prime }=x^{\prime \prime }\left ( t\right ) =c\left ( 1-x^{2}\right ) x^{\prime }\left ( t\right ) -kx\left ( t\right ) \end {array} \right \} \left . \begin {array} [c]{c}x_{1}^{\prime }=x_{2}\\ x_{2}^{\prime }=c\left ( 1-x_{1}^{2}\right ) x_{2}-kx_{1}\end {array} \right \} \]

The system of equations to be solved by ODE45 is the following

\[\begin {pmatrix} x_{1}^{^{\prime }}\\ x_{2}^{^{\prime }}\end {pmatrix} =\begin {pmatrix} x_{2}\\ c\left ( 1-x_{1}^{2}\right ) x_{2}-kx_{1}\end {pmatrix} \]

Subject to initial conditions \(x_{1}(0) =x(0) \) and \(x_{2}(0) = x^{\prime }(0)\). In the Matlab implementation below, the values of \(c,k\) and the initial conditions are defined at the top of the code. This needs to be modified to change the initial conditions before running the program again.

3 Results and output

The program was run for a number of different initial conditions (different \(x\left ( 0\right ) \) and different \(v\left ( 0\right ) \)). In all cases, \(c=1,k=1\) was used. It is noticed that one limit cycle is reached in all cases. Both the \(x\left ( t\right ) \) and the phase plot where shown. The following are output from 3 different runs made. The title on the plots show the initial conditions used.

3.1 Test case 1

3.2 Test case 2

3.3 Test case 3

4 Conclusion

A non-linear second order ODE was solved numerically using Matlab’s ode45. The solution to the Van Der Pol was found to contain a limit cycle in the phase portrait when starting from any initial conditions.

5 Source code

function nma_project_511() 
% This function solves the Van Der Pol nonlinear ode 
% numerically using Matlab's ode45. 
% 
% Final project for 511, CSUF, spring 2009. Instructor: Dr Oh, Sang June 
% 
% Written by : Nasser M. Abbasi 
% Van Der Pol ode is  x''(t) - c (1-x^2(t)) x'(t) + k x(t) = 0 
 
t  = 0:0.001:100;   % time scale 
c  = 1; 
k  = 1; 
x0 = 10; 
v0 = 3; 
 
[t,x] = ode45( @rhs, t, [x0,v0]); 
 
subplot(2,1,1); 
plot(t,x(:,1)); 
xlabel('t'); ylabel('x(t)'); 
title(sprintf('Van Der Pol, position vs, time, c=%3.2f, k=%3.2f',c,k)); 
 
subplot(2,1,2); 
plot(x(:,1),x(:,2)); 
xlabel('x(t)'); ylabel('v(t)'); 
hold on; 
plot(x0,v0,'*r','MarkerSize',10); 
title(sprintf('Phase portait showing limit cycle. x(0)=%3.2f, v(0)=%3.2f',x0,v0)); 
axis equal 
hold off; 
 
    function dxdt=rhs(t,x) 
        dxdt_1 = x(2); 
        dxdt_2 = c*(1-x(1)^2)*x(2)-k*x(1); 
 
        dxdt = [dxdt_1 ; dxdt_2]; 
    end 
end