Assignment 6: Initial Value Problem

Solve the ordinary differential equation with given initial values;
This is a van der Pol's equation with the initial values that the solution approximately starts on the limit cycle [Shampine, "Evaluation of a Test Set for Stiff ODE Solvers", ACM Transactions on Mathematica Software v. 7 #4 (1981) pp 409-420]. The time interval t choosen to be between [0, 3000] by Dr. Trangenstein was reasonable since the period of a cycle for this problem was about 1615.5 which means the interval was sufficiently long enough to be displaying several internal boundary layers.

This program outputs several data files that compare the convergence of the solution for different tolerance. The data in the files including the following values:
• Values of each time step;
• Values of y(t) at each time step;
• Values of z(t) at each time step

1. Backward Differentiation Formulas
This formula approximate the differential equation

by differentiating the polynomial that interpolates y at tn+1, tn, …, tn-k-+1.
In this way, the scheme has the form if

Let the coefficient of yn+1 to be 1, we have the following alternation of the expression

This method suggest the scheme has order k with error constant

2. Predictor-Corrector Method
Implicit linear multistep methods may be used for the initial value problems with rapid decay which could involve solving the nonlinear equation for yn+1. Fixed-point iteration is a good way to solve the nonlinear equations.
The predictor-corrector method uses Adams-Bashforth method of the same order k or order k-1 to estimate the initial guess for yn+1.

This backward differentiation formula is implemented by eoside.f due to Byrne and Hindmarsh and later modified in 1975. The iterative correction method is set to be “chord or semi-stationary newton method with a close form jacobian, which is computed in the user supplied subroutine pederv”.
This program is coded in VC 2005 with Intel Fortran Complier. The basic idea of calling Fortran in C++ is through dynamic link library (dll). The procedure of programing for this problem can be stated as:
• Compile epsode.f in Intel Fortran Compiler integrated in VC 2005. The compiler will generate two files, epsode.dll and epsode.obj, which are needed in the calling epsode subroutine;
• Add those two files in the epsode project. Define epsode function in C++ using “extern C”. Function “epsode” can generate the right result after proper variable input. Note that, the time step in epsode could be varied during iteration. In order to implement fixed time step algorithm, variable tout need to be reset when calling the epsode subroutine each time.;
• Code download

• Solutions with varies tolerance level
  
  
  
• Convergence analysis
  

Comments:
• The solution start to converge when to tolerance is smaller than 0.0001;
• The converged result shows a period of cycle of about 1615 which is consistence with the description of the result.