The Van der Pol Oscillator

Presented here is a numerical solution to the Van der PolOscillator. The form of the ODE’s has been chosen from Hairer et al. The equations are:

The Initial conditions y1(0)=-0.5, y2(0)=0, A = 2.008619861. The non-linearity in the equations is fairly clear from the plots below. The solution of y1 remains stable for several seconds between sudden changes in response, clearly following the behavior of y2. 

Advantages of the State Space Matrix Approach

The state-space approach requires an additional operation to perform the matrix inversion (or an LU-decomposition). Despite the complexity, the state-space method presents an elegant approach for ODE’s.

The form of the higher order ODE’s gets converted into 1^st order ODE’s, which can easily and readily be applied to any numerical solver. The state-space condensed matrix allows one to visualize the original nature of the problem (in K, C & M for the mass-spring-damper system). In other words, from an algorithmic stand point, all one needs to supply is the coefficient matrices along with the forcing functions and not worry about conversions to an equivalent format.

For a series of equations of higher order, the state-space method is thus more effective once the initial definition of the matrix is derived.

Solution of ODE’s using the State-Space Method

The state-space representation of

   provides a no-fuss approach to converting problems with 2^nd order ODE’s. The derivative equations take the form

 

The derivative definition thus requires a matrix inversion for A. 

 

Matrix inversion can be performed using techniques such as Gauss-Jordan. Subroutine gaussj from Numerical Recipes can be effortlessly used for the inversion. Rather than perform the matrix inversion, Press et al recommend a LU-decomposition followed by an LU-back substitution to solve for the differential form. They suggest that using the LU-decomposition is “a factor of 3 better than subroutine gaussj”.

Following the recommendation of Press et al, the preferred way to solve for the differential form is to assume it to be a set of linear equations of the form

 

The derivative form is then obtained by performing an LU-decomposition of Q, followed by the LU-backsubstitution of the LU version of Q with P. Subroutines ludcmp and lubksb work superbly efficient for this case.

The State Space Method (and Matrix Representations)

The 2^nd order differential equations describing the 3-DOF KCM system can be condensed into a matrix format as shown below: 

 

Once the component matrices for Mass, Stiffness, Damping and the Forcing Functions are defined, the overall system can be further condensed into a matrix format as: 

 

Letters in bold-italic are now vectors in matrix form. The response vectors and forcing functions are single-column matrices, while the KCM ones are square. Working with these matrices along the same lines as mentioned earlier, one can chop it up into two first order ODE’s, which can be “matricized” as follows:

 

The logic behind the above representation becomes fairly apparent when one expands the matrices to two separate equations. To state the obvious, what was originally a single set of matrix-equation in 2^nd order has now doubled.

Performing one last condensing, the above matrix then takes the form : 

 

  

where

           

The result is a matrix of 1^st order ODE’s that can be numerically integrated to obtain the solution of vector U. Such a representation of the system of 2^nd order equations is commonly referred to as the State-Space Method.

Algorithmic representation – An Example

Consider, for example, a 3-degree of freedom mass-spring-damper system as shown below. 

Each of the masses (m1 thru’ m3) is connected to springs (K1-K3) and dampers (C1-C3), with K3 and C3 connected to ground. Each of these masses can be subject to a time-dependent forcing function (F1-F3). The resultant displacement of each of these masses due to the forcing functions can be described by the 2^nd order ODE’s by applying Newton’s laws of motion. 

 

  

Using the conversion mentioned earlier, the above three 2^nd order ODE’s result in the following six 1^st order ODE’s:

 

 

 

 

 

By rearranging the above 1^st order ODE’s to solve for and , the numerical solutions can be readily obtained.