I'm appealing to the experience of the people in this newsgroup in
the hope someone here could give a hint which would take us to the
solution of this tricky problem.

I know it's a long message, but I didn't want to omit any information
that maybe could be important.

Here is a concise description of the problem:

Using FEM to solve a set of coupled differential equations we found
out that authentic solutions get polluted by "spurious" (strong
oscillating) ones if their eigenvalues are close within about 10%.

And now the details:

We are using (or rather trying to) the Finite Element Method to solve
a set of four coupled diferential equations

     -d/dr FU(r) + (k/r) FU(r) + VPS(r) GU(r) + DELTA(r) GV(r) = e GU(r)
      d/dr GU(r) + (k/r) GU(r) + VMS(r) FU(r) + DELTA(r) FV(r) = e FU(r)
      d/dr FV(r) - (k/r) FV(r) - VPS(r) GV(r) + DELTA(r) GU(r) = e GV(r)
     -d/dr GV(r) - (k/r) GV(r) - VMS(r) FV(r) + DELTA(r) FU(r) = e FV(r)

in a box (0 <= r <= R). GU, FU, GV, FV are the sought solutions and
VPS, VMS and DELTA are known functions, k is an integer parameter.

Boundary conditions are:

        GU(0) = GV(0) = FU(0) = FV(0) = 0
        GU(R) = GV(R) = 0

We were using the shooting method to solve this problem, but for the
case where DELTA is an integral operator, the above equations are
integro-differential equations and the shooting method fails. So we
thought we could use FEM to solve this.

We started of course solving (with the FEM) a simpler case which can
also be solved with the shooting method, for comparison.

After discretization using linear base functions (but we also used cubic
*Hermite* functions) we got the following generalized eigenproblem:

        A X = e B X

First trouble: we started to get what people seem to call "spurious"
solutions, strong oscillating functions. I read in a FEM book that
these oscillating solutions are a natural consequence of the fact that
the FEM is exact at the nodes, but just there. Ok, no problem per se,
since these "spurious" solutions are easy to identify by the program.

But, for certain choices of the functions VMS, VPS and DELTA it does
happen that a spurious solution has an eigenvalue close to the one of
an authentic solution. In this case the authentic solution seems to get
polluted by the spurious one and it gets some undesired oscillation.

That's a well known problem, we thought. In this case we just have to
look for the "smoothest" linear combination of the two vectors, the
spurious one and the "polluted" one, and that would be our "purified"
solution. BUT then we found out that the polluted solution already
represents the smoothest linear combination!

Still worse: if a spurious eigenvalue lies within 10% (!) of a authentic
one (in the test case we can always check with the shooting method if
the solution is ok), then the authentic solution gets polluted.

We have used three different diagonalization routines to solve the
problem: one from a FEM book, ARPACK (both sparse) and a Householder
algorithm with full matrices and all gave us the same results.

The first one, which we use at most, uses a bisection method the get the
eigenvalues, followed by an inverse iteration to get the eigenvectors. I
read recently an article by Peter Spellucci where he warns about this
procedure in the case of clustered eigenvalues. Could that be the source
of our problems? Can "equal within 10%" be considered clustered?

[Just to give you a more concrete idea of the problem: we ask for
eigenvalues in the window [0, 100] and the eigenvector correspoding to
an eigenvalue 27.26 seems to get polluted by a spurious eigenvector
with eigenvalue 24.75!]

What can be done to remove the influence of the spurious solutions in
the authentic ones? Better yet, is there a way to eliminate this strong
oscillating solutions? (After using cubic Hermite base functions I'm
afraid one has to live with them!)

Thank you for your patience reading this. Your comments are very
welcome.

Emilio.

|>
|> in a box (0 <= r <= R). GU, FU, GV, FV are the sought solutions and
|> VPS, VMS and DELTA are known functions, k is an integer parameter.

also e is sought !

|>
|> Boundary conditions are:
|>
|>         GU(0) = GV(0) = FU(0) = FV(0) = 0

|>         GU(R) = GV(R) = 0
I am not well acquainted with the theory of these systems, but would think that
here only one condition is possible (substitute HU(r)==e, with the equation
 (d/dr ) HU(r) == 0 ,
then you have a system of 5 ODE's and five boundary conditions,
this should suffice to define a locally unique solution. )
otherwise there might be locally nonunique solutions, but then you can  
forget your solution method

|>
|> We were using the shooting method to solve this problem, but for the
|> case where DELTA is an integral operator, the above equations are
|> integro-differential equations and the shooting method fails. So we
|> thought we could use FEM to solve this.
I assume that in the integral equation case you have something like
    DELTA(r)GV(r) = \int _0^r GV(rho)fct(rho) d (rho)  with fct given?
    hence a linear integro differential equation ??

|>
|> We started of course solving (with the FEM) a simpler case which can
|> also be solved with the shooting method, for comparison.
|>
|> After discretization using linear base functions (but we also used cubic
|> *Hermite* functions) we got the following generalized eigenproblem:
|>
|>         A X = e B X
|>
|> First trouble: we started to get what people seem to call "spurious"
|> solutions, strong oscillating functions. I read in a FEM book that
|> these oscillating solutions are a natural consequence of the fact that
|> the FEM is exact at the nodes, but just there. Ok, no problem per se,
|> since these "spurious" solutions are easy to identify by the program.
^^^^^^^^^^^ this is surely nosense. the oscillating effect comes from
            an inadequate grid, you can be sure. try the simplest case
             \epsilon y'' + y' = 0 , y(0)=1, y(\infty)=1
            with an equidistant grid !

|> We have used three different diagonalization routines to solve the
snip
|> The first one, which we use at most, uses a bisection method the get the
|> eigenvalues, followed by an inverse iteration to get the eigenvectors. I
|> read recently an article by Peter Spellucci where he warns about this
|> procedure in the case of clustered eigenvalues. Could that be the source
|> of our problems? Can "equal within 10%" be considered clustered?
no,   this effect occurs for very near eigenvalues only
I hope you have used a good method for the general eigenvalue problem
not transformation into the "standard" form, which in turn may introduce
large errors if B is illconditioned. In this case it is the mass matrix
and should be reasonable if anything (including integration) is done
correctly.
|> [Just to give you a more concrete idea of the problem: we ask for
|> eigenvalues in the window [0, 100] and the eigenvector correspoding to
|> an eigenvalue 27.26 seems to get polluted by a spurious eigenvector
|> with eigenvalue 24.75!]
I don't believe that the problem comes from this. If you use inverse iteration
for individual eigenvectors and have neighboring eigenvalues , then this
degrades speed and also leads to nonorthogonal (in your case B-orthogonal)
eigenvectors, to some degree. you can simply check that.

The equations shown above are what you get by substituting g(r)=G(r)/r
and f(r)=F(r)/r in the original physical equation. Because you want
the solutions to be well defined at r=0, you impose that G(0)=F(0)=0.
But it does *not* mean that the original solutions g(r)=G(r)/r and
f(r)=F(r)/r are zero at r=0, since G(r) und F(r) can have a linear
behaviour near r=0 (and that's indeed the case!).

So the aditional boundary conditions don't influence the physics of
the problem, they are just a mathematical artifact.

But still linear. In our test case fct involves a Dirac-Delta "function".

Initially I though I could be using wrong basis functions: I have a
first-order DEs but my first derivatives aren't continuous! Then I
tried a Hermite basis, where the first derivatives are continuous. Of
course the spurious solutions couldn't oscillate in that case, but the
"spuriousness" was still there, in a different way (can't explain with
words, but looking at the plots it was clear that some solutions were
still spurious).

So, could you elaborate on your comment about the grid? Do you think
the problem is the homogeneous grid? What kind of grid should I be
using? One with spacing proportional to r, I suppose from you comment
at the beginning? Could you help me to understand how the wrong choice
of grid can originate these spurious solutions?

The matrix B in this case is positive definite.

Thank you for your comments!

|>    DELTA GV(r) = \int_0^\infty GV(rho) fct(r,rho) d(rho)

 ^^^^^^^^^^ numerically you write "R" instead of \infty ?
|>
|> But still linear. In our test case fct involves a Dirac-Delta "function".

ooh ooh....    this introduces a step discontinuity in the differential equation?
never integrate numerically over such a thing. This necessarily introduces
severe oscillations.
the numerical results you report simply say that the eigensolver is o.k.
hence the trouble comes from the type of discretization.
I guess you should try the geometric grid, avoid any step or derivative
discontinuity _between grid points_ , linear elements are enough,
if this does not help, go back to your original equations and consult numerics
for singular perturbation problems. there are lots of papers on this subject.
maybe you have a series solution fro small r and can fit the two parts together,
this would be the finest solution, or you must go even deeper into these numerics

The original equations are these:

   -d/dr fu(r) + (k-1)/r fu(r) + VPS(r) gu(r) + DELTA(r) gv(r) = e gu(r)
    d/dr gu(r) + (k+1)/r gu(r) + VMS(r) fu(r) + DELTA(r) fv(r) = e fu(r)
    d/dr fv(r) - (k-1)/r fv(r) - VPS(r) gv(r) + DELTA(r) gu(r) = e gv(r)
   -d/dr gv(r) - (k+1)/r gv(r) - VMS(r) fv(r) + DELTA(r) fu(r) = e fv(r)

The boundary conditions:

        gu(R) = gv(R) = 0 for every k

and

        if (k == 1) then
           gu(0) = gv(0) = 0
        elseif (k == -1) then
           fu(0) = fv(0) = 0
        else
           gu(0) = gv(0) = fu(0) = fv(0) = 0
        endif

We do the transformation g(r) = G(r)/r and f(r) = F(r)/r, for both u
and v, and get

   -d/dr FU(r) + (k/r) FU(r) + VPS(r) GU(r) + DELTA(r) GV(r) = E GU(r)
    d/dr GU(r) + (k/r) GU(r) + VMS(r) FU(r) + DELTA(r) FV(r) = E FU(r)
    d/dr FV(r) - (k/r) FV(r) - VPS(r) GV(r) + DELTA(r) GU(r) = E GV(r)
   -d/dr GV(r) - (k/r) GV(r) - VMS(r) FV(r) + DELTA(r) FU(r) = E FV(r)

Now with

        GU(0) = GV(0) = FU(0) = FV(0) = GU(R) = GV(R) = 0 for every k.

            = DEL(r) GV(r).

As you see, it just makes the integral go away, so that we get the
equations I wrote above.

I'll try your suggestion for the grid, it'll surely improve things.

Thank you for the idea.