VFGEN § Taylor

Taylor Method

The idea of the Taylor method for approximating the solution to a system of differential equations is to expand the solution to
X'(t) = F(X) in a Taylor series: X(t+h) = X(t) + hX'(t) + (h²/2)X''(t) + (h³/6)X⁽³⁾(t) + ... The differential equation gives us X'(t). To find X''(t), we differentiate the differential equation with respect to t (t arguments are suppressed): X'' = DF(X)[X'] (We denote the r-linear differential D^rF(X) acting on the vectors V₁,...,V_r as D^rF(X)[V₁,...,V_r].) Differentiate again to find X⁽³⁾ and higher derivatives of X:

X⁽³⁾	=	D²F(X)[X',X']+DF(X)[X'']
X⁽⁴⁾	=	D³F(X)[X',X',X'] + 3D²F(X)[X'',X'] + DF(X)[X⁽³⁾]
X⁽⁵⁾	=	D⁴F(X)[X',X',X',X'] + 6 D³F(X)[X'',X',X'] + 3D²F(X)[X'',X''] + 4D²F(X)[X⁽³⁾,X'] + DF(X)[X⁽⁴⁾]

and so on.

The VFGEN Taylor Command

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

$ vfgen taylor:order=r vector_field_file.vf

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^kF(X)[V₁,...,V_k]. These could be useful in programs that analyze bifurcations or that compute invariant manifolds.

Options

order

This option specifies the highest order retained in the Taylor series. The order must be a positive integer.
Default: order=5

Bugs and Limitations

The length of the expression for the n-th derivative of an expression can grow rapidly as n increases. Currently, the program computes derivatives and prints the formulas to the C file without checking the length. Some C compilers may not be able to handle the long lines.
Some of the C variable names that are generated can also become very long. This may cause some compilers to fail.
The program assumes that the vector field is autonomous. (This should be fixed in a future version.) It will not warn you if you use this command on a nonautonomous vector field, but the results returned by the C function will not be correct.

Example 1

Here is a sample vector field file: pendulum.vf.
The file created by

$ vfgen taylor:order=7 pendulum.vf

are

A short C program that uses these files is pend_main.c. A Makefile for the pend_main program is Makefile_pend_main. The program pend_main uses the function pendulum_taylor7 to implement a simple solver for the differential equations. Here is a plot of the output of pendmain:

The following plot shows a comparison of the polynomial approximations for orders 3, 7, and 15 for the pendulum vector field. The Taylor expansion was found for θ(0)=1.0 and v(0)=1.5, with parameter values g=9.81, b=0.25, L=1, and m=1. (For example, the data for the order=3 case was created with pend_taylor3_demo.c.) Also shown is the solution generated with the GSL ODE solver (with abserr=10^-12 and relerr=10^-10).

Example 2

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

$ vfgen taylor:order=8 vanderpol.vf

A Makefile for the solver is Makefile-vanderpol_adaptive8; the solver is compiled and run as follows:

$ make -f Makefile-vanderpol_adaptive8
$ ./vanderpol_adaptive8 > vdp8.dat

The 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.dat

A solver that uses RADAU5 is created with

$ vfgen radau5:demo=yes vanderpol.vf

and 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.

The plots are indistinguishable. The adaptive Taylor method generates 4767 data points, the GSL solver generates 4904, and the RADAU5 solver generates just 1191 points. The final values of the x and y coordinates of the numerical solutions are

Solver	x(10)	y(10)
Adaptive Taylor (Order 8)	1.80308347	0.150123424
GSL	1.80308347	0.150123395
RADAU5	1.80308348	0.150123481