Embedded Files
Solve y'(t)=z(t), z'(t)=1000(1-y(t)^2)z(t)-y(t) with y(0)=2, z(0)=0 for t between 0 and 3000. This problem has internal boundary layers, and the initial condition is nearly on the limit cycle, which has a perior of about 1615.5. Use EPSODE to solve this problem, and examine how the numerical solution depends on the choice of error tolerance. You can find some discussion of this problem in Shampine, "Evaluation of a Test Set for Stiff ODE Solvers", ACM Transactions on Mathematica Software v. 7 #4 (1981) pp 409-420.
Solution:
Backgroud:
van der Pol oscillator
y'' - mu(1-y^2)y' + y = 0
and let y' = z, then z' = mu(1-y^2)z - y. van der Pol equation has a unique stable limit cycle for each mu > 0.
Fig.(wiki) Phase portrait of the unforced Van der Pol oscillator, showing a limit cycle and the direction field.
Source Code:
This old version epsode.f (NOT used) is a stiff ordinary differential equation initial-value problem solver by Bryne and Hindmarsh. "This program uses Adams-Bashforth/Adams-Moultonpredictor/corrector methods of order at most 12 for non-stiff problems, and backward differentiation formula of order at most 5 for stiff problems." (Trangenstein, Scientific Computing)
This new version epsode.f (taken from Prof. Trangenstein's server) together with d1mach.f and VanderPol.c can be compiled by "./compile" using this script compile. The error tolerance "eps" needs to be modified in the code VanderPol.c. Run "./VanderPol", then the program outputs data (t,y,z) in the file "trajectory.dat". The results using eps = 0.000000001 and 10000 time intervals is here. One can plot y-z, y-t or z-t by these gnuplot scripts plotyz, plotyt and plotzt.
Results:
y-z phase portrait with different error tolerances. Trajectory starts to converge for eps <= 10^(-5)
The y-t and z-t plots at eps = 10^(-9) are
The stiffness of the problem is shown by the separation of time scales. While y(t) and z(t) vary slowly within some time intervals, y(t) exhibits rapid "jumps" and z(t) exhibits rapid "spikes" within very short time. On the limit cycle of the y-z phase portrait, this corresponds to two slow curved edges (top left and bottom right) and two fast straight edges (top right and bottom left).
Google Sites
Report abuse