NAUtil aims at automating the **stability analysis of linear finite difference schemes.**
It has been developed to treat electromagnetic problems but can be practically extended to any application provided it is discretized by a linear finite difference scheme.
__Collaboration in any application field is seeked.__

The different routines and samples have to be used in a Maple environment. Maple 10 has been used for the development. Routines have not been tested on earlier Maple versions, but it is not totally impossible that it works also.

Since we only deal with linear models, we can analyze them in the frequency domain.
Thus we assume that the scheme handles a single (vector valued) variable U_{j}^{n} with spatial dependence U_{j}^{n} = U^{n} exp(iξj).
The scheme is then described as U^{n+1}=GU^{n} and in our case the amplification matrix G does not depend on time or on the space steps δx and time step δt separately but only on the ratio δx/δt.
This ensures that U^{n}=G^{n}U^{0} and stability comes to the boundedness of ^{n}.

A necessary stability condition is that the eigenvalues of G lie in the unit disk. Only eigenvalues on the unit circle can cause instabilities and two cases can occur:

- (a)
- the eigenvalues of modulus 1 are simple and the scheme is stable. Stability analysis is performed on the characteristic polynomial of G using von Neumann analysis;
- (b)
- some eigenvalues of modulus 1 are multiple and stability is obtained if and only if the associated minimal subspaces are of dimension 1. The analysis of the characteristic polynomial does not yield any information on the associated subspaces. Matrix G has to be studied. In our study, this corresponds to degenerate cases for which matrix G has always a very specific form and minimal subspaces are easy to determine.

The localization of roots of polynomials is a difficult task if the polynomial has a large degree and/or coefficients depending on many parameters. This the case in our context. We split this difficult problem into many simpler ones.

Given a polynomial φ_{0}, we construct a sequence

φ_{m+1}(z) =
( φ_{m}^{*}(0)φ_{m}(z)
- φ_{m}(0)φ_{m}^{*}(z) ) / z.

This sequence has a decreasing order,
**Definition:**
A __Schur polynomial__ has all its roots r inside the unit circle: |r|<1.

**Theorem:**
A polynomial φ_{m} is a Schur polynomial of exact degree d if and only if φ_{m+1} is a Schur polynomial of exact degree d-1 and |φ_{m}(0)|≤|φ_{m}^{*}(0)|.

**Definition:**
A __simple von Neuman polynomial__ has all its roots r in the unit disk: |r|≤1 and eigenvalues of unit modulus are simple.

bB>Theorem:
A polynomial φ_{m} is a simple von Neumann polynomial if and only if

(i) φ_{m+1} is a simple von Neumann polynomial and |φ_{m}(0)|≤|φ_{m}^{*}(0)|;

or (ii) φ_{m+1} is identically zero and φ_{m}' is a Schur polynomial.

At each step, we have to check condition |φ_{m}(0)|≤|φ_{m}^{*}(0)| and degree.
This is comparatively simpler than locating the roots.
But these are not simple problems either.
We now have to determine the sign of a polynomial with coefficents depending on the parameters.
These parameters often lie in an interval, e.g. the ratio δt/δx or δt/δx^{2} (the stability condition often reads like this); or are positive, e.g. physical constants.

- July 13th 2006: First β-version online
- September 2006: Version 1.0 online