X^{(3)} | = | D^{2}F(X)[X',X']+DF(X)[X''] |
X^{(4)} | = | D^{3}F(X)[X',X',X'] + 3D^{2}F(X)[X'',X'] + DF(X)[X^{(3)}] |
X^{(5)} | = | D^{4}F(X)[X',X',X',X'] + 6 D^{3}F(X)[X'',X',X'] + 3D^{2}F(X)[X'',X''] + 4D^{2}F(X)[X^{(3)},X'] + DF(X)[X^{(4)}] |
The taylor command of VFGEN creates a function that computes the Taylor expansion of X(t+h) to a given order. The order is specified with the order option. The general form of the command is
The files created are [name]_taylor[r].c and [name]_taylor[r].h, where [name] is the name of the vector field given in the vector field file, and [r] is the order given with the order option. The C file [name]_taylor[r].c will contain the functions [name]_vf and [name]_taylor[r]. The first computes the vector field, and the second computes the Taylor approximation. Also defined in the C file are the functions [name]_diff[k], for k=1, 2, ..., r-1. These functions compute the k-linear differentials D^{k}F(X)[V_{1},...,V_{k}]. These could be useful in programs that analyze bifurcations or that compute invariant manifolds.
order |
This option specifies the highest order
retained in the Taylor series.
The order must be a positive integer.
Default: order=5 |
A vector field file for the van der Pol equation is vanderpol.vf. The file vanderpol_adaptive8.c implements an adaptive solver that uses an order 8 Taylor method to solve the differential equation. The code works by using the Taylor polynomial at x(t) to estimate x(t+h). Then, the Taylor polynomial at x(t+h) is used to compute the value at x(t). If this value differs from the original x(t) by more than tol, the step size is reduced by 2/3, and the test is tried again. If the difference is less than tol, the solution is advanced. If, in addition, the solution is less then mintol, the step size is increased by a factor of 1.5.
(This method was implemented as a fairly simple test of the code generated by VFGEN. There are smarter ways to adapt the step size--see, for example, the documentation of the Taylor program by Zou and Jorba.)
The program uses the functions defined in the C files created by the command
$ make -f Makefile-vanderpol_adaptive8 $ ./vanderpol_adaptive8 > vdp8.datThe parameter values used in the file vanderpol_adaptive8.c are: ε=10^{-3}, tol=10^{-8}, mintol=10^{-10}.
For comparison, a solver that uses the GNU Scientific library is created (also with VFGEN):
$ vfgen gsl:demo=yes vanderpol.vf $ make -f Makefile-vanderpol $ ./vanderpol_solve epsilon=0.001 relerr=1e-7 abserr=1e-9 > vdpgsl.datA solver that uses RADAU5 is created with
$ vfgen radau5:demo=yes vanderpol.vfand then the driver program vanderpol_dr5.f is modified to use ε = 0.001, and the relative and absolute tolerances are set to 10^{-7} and 10^{-9}, respectively. The RADAU5 driver is compiled as shown in the RADAU5 section of the web page.
Solutions generated by the GSL solver and by the adaptive Taylor solver are shown in the following plot. Only the interval 8.5 < t < 10 is shown.
Solver | x(10) | y(10) |
Adaptive Taylor (Order 8) | 1.80308347 | 0.150123424 |
GSL | 1.80308347 | 0.150123395 |
RADAU5 | 1.80308348 | 0.150123481 |