TP algorithm

Definition of the order parameter. The first step in the TPS algorithm is to define an order parameter, i.e. a variable that describes whether the system is in the reactants, products or in an intermediate region, as shown schematically in Fig. 21. The pyruvate and lactate regions were defined by the values of the appropriate bond... 

Fig. 21 Schematic representation of the TPS algorithm. The shaded regions are identified by the order parameter as reactants and products. The solid line is a reactive trajectory. A shooting move is shown a time slice was chosen along the reactive trajectory, momenta were perturbed, and then the system was propagated forward and backward in time, resulting in a non-reactive trajectory shown with the dashed line.

It is expected that the TP algorithm will give better accuracy than the FD and BD algorithms. However, in addition to accuracy, stability is also very important. A stable algorithm means that the solution value will remain finite as the number of solution steps increases indefinitely. Stability can be usefully studied by looking at the simple linear differential equation ... 

From an accuracy viewpoint, the TP algorithm appears very attractive. However, for this example the TP algorithm has some problems when the step size becomes too large as in fact occurs for the time step of l.e-2 shown in Figure 10.7. The curves show that the accuracy for small values of time is essentially equal to the solution value. Even more problematic is the fact that the TP solution tends to oscillate about the true solution for the 1. e-2 case as seen near a time of about 6. 

The algorithm uses the TP algorithm except for the first time interval which is subdivided into 4 sub intervals and the BD rule is used for each of the four sub intervals. These four sub intervals are handled by the repeat-until loop from lines 32 to 45. The first time this loop is encountered, it is repeated four times as jfirst is initially set to zero on line 7. Subsequent encounters of this loop cause only one execution of the loop. Before entering the main loop and the repeat-until loop for the first time, three parameters h, h2 and hx are set on line 29 to values appropriate to the BD algorithm. After the repeat-until loop is executed four times, these parameters are set to the appropriate TP parameters on line 46 for all subsequent time points. Even fiiough flic first time interval is subdivided into 4 intervals, only the solution values for the last of flic four sub intervals is saved in flic output files on lines 47 and 48. The calculation of flic solution values for flie differential equations requires only four basic steps for each desired time interval. On line 36 the yn parameter needed in the derivative approximation is evaluated. On line 37 a predicted value of flic solution variables is calculated. On line 39, time is incremented to the next time point. On line 41 nsolvQ is called to solve the set of implicit equations for flic new solution values. On line 44 new values of the derivatives are evaluated so they can be used in flic next iterative step on line 36 in the evaluation of the yn parameters. 

A More Detailed Look at Accuracy Issues with the TP Algorithm... 

Equation (10.33) is suggestive that the error in the TP algorithm is proportional to the product of the third derivative of the function being calculated and h and that by decreasing the step size by a factor of 10 the solution accuracy would increase by three orders of magnitude. However care must be exercised in using... 

Figure 10.17 Absolute errors in the numerical solution of second order equation with sinusoidal solution. Solid curve is numerical error from numerical integration with TP algorithm. Open and solid points are theoretical errors.

From the theory presented above it is expected that at any given time for a solution, the error in the TP algorithm varies with the square of the step size. These concepts can now be tested in more depth for the TP integration code and the ac-... 

Figure 10.21. Errors in a sinusoidal solution with TP algorithm obtained in Listing 10.12. Both the uncorrected and corrected solution errors are shown.

Listing 10.16. Example eode for replaeing TP algorithm with fourth order RK algorithm. 

In the above ERRj is some desired maximum relative error in the solution (sueh as l.e-5), is the old step size used in evaluation die relative error and Kew estimated step size needed to give a relative error equal to the desired maximum error of ERR. In the last equation it can be seen that if die estimated relative error is larger dian die desired maximum error the new step size will be appropriately smaller than the old step size. In this maimer the local step size can be adjusted as the solution progresses to achieve a desired relative error criterion. It should be noted that the power of 2 and 0.5 in the above equations is based upon the TP algorithm. For use with the RK algorithm die appropriate values would be 4 and 0.25 respectively since the error varies as the fourth power of the step size. For multiple variables widi several differential equations, a relative error can be associated widi each variable in die same form as die above equations. In such a case the criterion should use the largest calculated relative error among all the variables at each time step so that the relative error of all variables will satisfy some desired criterion. 

The errors aehieved with the remaining two MATLAB functions are shown in Figure 10.33. These two funetions are recommended for use with non-stiff differential equations and the results indieate that they do not provide a good solution for this example. However, the aetual error achieved is not very different from the fimetions shown in Figure 10.32 whieh are supposed to be better able to handle stiff differential equations. Shown in each figure is the error aehieved by the adaptive step size eoupled with the TP algorithm funetion odeivs() developed in this work. It is seen that for the default parameters, the aehieved error of odeivs() is better than that aehieved by any of the MATLAB funetions. 

For the TP algorithm one needs to replace the first and second derivative terms by quantities already known at solution point n and by the solution value at point +1. This can be doe if these equations are written in the form ... 

Listing 10.19. Code segment for direct solution of systems of second order differential equations by TP algorithm... 

---