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 Ujn with spatial dependence Ujn = Un exp(iξj).
The scheme is then described as Un+1=GUn 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 Un=GnU0 and stability comes to the boundedness of
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:
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
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.
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/δx2 (the stability condition often reads like this); or are positive, e.g. physical constants.