**NumericalMath****`****OrderStar****`**

Analysis of the stability of numerical methods for solving differential equations is a considerably more difficult task than finding the order of approximation. This difficulty is reflected in the fact that there are many different ways of defining stability. This package assists in the determination of stability regions for numerical methods.

Stability regions are important because they reflect the rate at which errors are propagated in our approximate solution. Just as there are absolute and relative errors in numerical analysis, there are also absolute and relative measures of stability. This package renders order stars, which are useful in examining the relative stability of a method. Furthermore, by specifying the comparison function to be identically 1, you can also draw regions of absolute stability.

A given numerical method for a problem can be recast into the framework of approximation theory. The goal is then to study how well this approximant behaves when compared with the solution. There is a kind of paradox here, for if the solution were known then we would have no need to resort to a numerical approximation. However, we want to establish a framework which applies to any problem in a given class. Since generically analytic solutions to problems cannot be found, it is common to study how a numerical method behaves when it is applied to a linearized system. In the area of ordinary differential equations, for example, we are interested in solutions to the system of equations:

for a generally nonlinear

. It is common to replace this system by a scalar linear problem that we can solve, namely,

where is considered to be a complex constant and we have fixed an initial condition so that the equation is uniquely determined. Stability analysis is now a study of how well a numerical solution behaves when applied to the simplified differential equation (1 ). Equation (1

) is often referred to as the scalar linear test problem, or Dahlquist's equation.

The discussion which follows concentrates on how to use the package. You should keep in mind that although the focus is on the behavior of an approximant, our underlying interest is in the numerical method from which the approximant arose when applied to some problem. For more information on this correspondence, stability analysis, and the theory of order stars, see the references at the end of this section.

Drawing order stars.

Padé approximations are rational polynomial approximants where all parameters are chosen to maximize order at some local expansion point. Certain numerical methods such as Runge-Kutta methods are related to Padé approximants to the exponential.

This loads the package.
In[1]:= **<<NumericalMath`OrderStar`**

This constructs a Padé approximant to

. The package for doing this is loaded automatically. This approximant corresponds to the forward Euler method.
In[2]:= **approx = Pade[ Exp[z], {z, 0, 1, 0} ]**

Out[2]=

This is the relative stability region, or order star of the first kind, for the forward Euler method. The pole of the approximant is highlighted.
In[3]:= **OrderStar[ approx, E^z ]**

OrderStar::sols: Warning: No finite zeros of function found using NSolve. Either inverse functions or transcendental dependencies were involved. Try specifying omitted points using options.

This is the absolute stability region for the forward Euler method, obtained as a relative comparison with 1.
In[4]:= **OrderStar[ approx, 1 ]**

Options unique to OrderStar.

When you ask for certain features to be displayed and Mathematica is unable to find these features, you will obtain the message OrderStar::sols containing more specific information relating to the problem. Solve may also issue messages such as when inverse functions are being used.

OrderStar uses heuristics in order to determine what the independent variable is. You can save time in a very complicated expression by specifying the variable to use explicitly. If there is any ambiguity in the variable choice, then input returns unevaluated and an appropriate warning message is issued, since the function will not evaluate numerically.

This indicates the variable to use and highlights points where

. This may not be possible in general if the relationship is nonalgebraic.
In[5]:= **OrderStar[approx, Exp[z], z,**

OrderStarInterpolation->True]

OrderStar::sols: Warning: No finite zeros of function found using NSolve. Either inverse functions or transcendental dependencies were involved. Try specifying omitted points using options.

In addition to True and False, the options OrderStarInterpolation, OrderStarLegend, OrderStarPoles, and OrderStarZeros can take on lists of coordinate pairs to specify points that cannot be found automatically. As well as resizing the plot legend by specifying scaled coordinates, you can specify information to the legend such as the style and size of the font to use.

The position of the legend is given in scaled coordinates using the same syntax as that of Rectangle. Font style and size information is also specified and the symbols used to represent zeros and poles are increased in size.
In[6]:= **OrderStar[ Pade[ Sinh[z-1], {z, 0, 3, 3}], Sinh[z-1],**

OrderStarLegend -> {{.6, .6}, {.98, .98}},

DefaultFont -> { "Courier", 6.5},

OrderStarSymbolSize -> 0.02 ]

In addition to the many options unique to OrderStar, there are several that are simply passed to ContourPlot and used to produce the plot.

Options common to many graphics functions.

An important issue is whether an order star crosses the imaginary axis. Additionally, it may be interesting to illustrate symmetry around the real axis. In order to facilitate these comparisons OrderStar uses graphics options to render axes which pass through the origin.

The default plot region and subdivision are determined from essential features of the order star. However, this default can be overridden using standard Graphics options.

This defines a function.
In[7]:= **fun = Exp[Cos[z]+ I Cosh[z]];**

This constructs a Padé approximant to the function.
In[8]:= **approx = Pade[ fun, {z, 0, 2, 3}];**

We can change the default plot region and plot density using standard options. However, in doing so we can make the plot quite jagged. Increasing PlotPoints can resolve this issue.
In[9]:= **OrderStar[approx, fun,**

PlotRange -> {{-5, 5},{-5, 5}},

PlotPoints -> 40 ]

OrderStar::sols: Warning: No poles of function found using NSolve. Either inverse functions or transcendental dependencies were involved. Try specifying omitted points using options.

OrderStar::sols: Warning: No zeros of function found using NSolve. Either inverse functions or transcendental dependencies were involved. Try specifying omitted points using options.

Order stars of the second kind are useful in the study of stability for linear multistep methods.

We can also plot the order star of the second kind. Note the poor resolution near the point of expansion

. There is not a similar problem near the origin because extra work is done there by default.
In[10]:= **OrderStar[ Pade[ Log[1+z], {z, 1, 3, 2}], Log[1+z],**

OrderStarKind -> Second,

PlotPoints -> 50,

PlotRange -> {{-.5, 2},{-1, 1}} ]

By default OrderStar resolves fine features at the origin by overlaying a subplot. However, there may be times when you wish to disable this feature, or to resolve fine features in other specified regions. Increasing the overall PlotPoints setting is an inefficient solution to this problem. The option OrderStarSubPlots provides a mechanism for doing this.

By overlaying smaller contour plots in regions of poor resolution the figure is much improved. If no PlotPoints option is specified for a particular subplot, a fraction of the value for the main plot is used.
In[11]:= **OrderStar[ Pade[ Log[1+z], {z, 1, 3, 2}], Log[1+z],**

OrderStarKind -> Second,

PlotPoints -> 50,

PlotRange -> {{-.5, 2},{-1, 1}},

OrderStarSubPlots ->

{{PlotRange -> {{0.9, 1.1},{-0.1, 0.1}},

PlotPoints -> {20, 20}}} ]

Order stars provide a means of determining, at a glance, many important features of interest such as A-stability. Furthermore, by considering a relative comparison with a function, order stars manage to encrypt the order of accuracy of a numerical scheme into the stability region.

Here we include only the zeros and poles of the approximant, since there are no finite zeros or poles of

. The numerical method corresponding to this approximant is A-stable since the approximant has no poles in the left half-plane. Furthermore a count of the sectors adjoining the origin tells us that the order of approximation is 5 (one less than the number of sectors).
In[12]:= **OrderStar[ Pade[ Exp[z], {z, 0, 2, 3}], Exp[z],**

OrderStarPoles -> {True, False},

OrderStarZeros -> {True, False} ]

Expositions of the theory of stability of numerical methods can be found in , , and . More examples of applications of the package are provided in

.

**References**

E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems

, Springer, Berlin, 1991.

A. Iserles and S. P. Nørsett, Order Stars

, Chapman and Hall, London, 1991.

J. D. Lambert, Numerical Methods for Ordinary Differential Equations: The Initial Value Problem

, John Wiley and Sons, Chichester, 1991.

M. Sofroniou, "Order Stars and Linear Stability Theory", Journal of Symbolic Computation, 1996, 11

, pp. 1-31.