Eigenvector "pollution", spurious FEM solutions, etc.
Emilio Lopes
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.
--
Emilio C. Lopes
Physik-Department (T36) mailto:Emilio...@Physik.TU-Muenchen.DE
Technische Universitaet Muenchen http://WWW.Physik.TU-Muenchen.DE/~ecl
Germany
Peter Spellucci
Emilio Lopes <Emilio.Lo...@Physik.TU-Muenchen.DE> writes:
snip
|> 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)
but numerically r near zero will introduce considerable trouble
use a fine grid with geometric refinement at r=0
|>
|> 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.
hope this helps
peter
Emilio Lopes
> |> -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)
> ^^^^^^^^^^^^warning: the homogeneous initial values make this well defined
> but numerically r near zero will introduce considerable trouble
> use a fine grid with geometric refinement at r=0
You mean the term k/r, right?
> also e is sought !
That's correct, I forgot to mention it.
> |> 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 [...]
That's right, we have "too much" boundary conditions. Physically one
needs only four (either G(0) or F(0) = 0, but not both), instead of
the six we have now.
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.
> |> 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 ??
DELTA GV(r) = \int_0^\infty GV(rho) fct(r,rho) d(rho)
But still linear. In our test case fct involves a Dirac-Delta "function".
> |> 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 !
It's nice to hear these spurious solutions aren't "natural"! I have
seem many FEM books and none but this one mentions this oscillating
"solutions".
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?
> 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.
The methods I tried (except the Household with full matrices) do not
transform the problem into the standard form.
The matrix B in this case is positive definite.
Thank you for your comments!
Emilio.
Bob Roach
discretized system. They just don't happen to be the solution you wanted. They
are also a natural consequence of the polynomial basis functions one uses for FEM
or FD or FV methods. May I suggest you try to use a basis function which arises
from the solution to the linearized problem. At least this way the basis function
will include real eigenvalues of the linear system which will more closely mimic
the solution you want than those which arise from using polynomial bases. I have
found that solutions from such bases (when they can be found) are monotonic,
robust, and extremely accurate.
good lucks,
Bob
Peter Spellucci
Emilio Lopes <Emilio...@Physik.TU-Muenchen.DE> writes:
|> Peter Spellucci wrote:
|>
|> > |> -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)
|> > ^^^^^^^^^^^^warning: the homogeneous initial values make this well defined
|> > but numerically r near zero will introduce considerable trouble
|> > use a fine grid with geometric refinement at r=0
|>
|> You mean the term k/r, right?
problem like
(d/dr)(g(r)*r) = ..............
and transformed this into "a better looking" system as above. Hence the
substitution you mentioned below?But for singular perturbaton problems,
there exist specail metohds. you cannot solve such equations by the standard
approach. using a geometric grid is the minimum which is necessary. hence use a
grid rule like
delta_r = delta0 if r>=r0, (with r0=1, say)
and
delta_r = delta_r/2
r=r-delta_r if ( r<r0, as long as r>=rleft, with a very small rleft.
that means a geometric progression to zero at the left end of the interval.
|> > |> 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 [...]
|>
|> That's right, we have "too much" boundary conditions. Physically one
|> needs only four (either G(0) or F(0) = 0, but not both), instead of
|> So the aditional boundary conditions don't influence the physics of
|> the problem, they are just a mathematical artifact.
|> 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 matrix B in this case is positive definite.
that means that the trouble will not come from the eigensolver, since you are
using
already state of the art software for this.
hope this helps
peter
Emilio Lopes
> what is called a singular perturbation problem like
> (d/dr)(g(r)*r) = .............. and transformed this into "a
> better looking" system as above. Hence the substitution you mentioned
> below?
No.
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.
> |> DELTA GV(r) = \int_0^\infty GV(rho) fct(r,rho) d(rho)
> ^^^^^^^^^^ numerically you write "R" instead of \infty ?
Yes. There is a limit for R so that larger values don't change the
results, because 1) the functions VPS(r), VMS(r) and DELTA(r) go to zero
for large values of r. 2) We have a cutoff on e.
> |> In our test case fct involves a Dirac-Delta "function".
> ooh ooh.... this introduces a step discontinuity in the differential
> equation?
DELTA GV(r) = \int_0^\infty GV(rho) DEL(rho) \delta(r-rho) d(rho)
^^^^^^^^^^^^^^^^^^^^^^ fct(r,rho)
= 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.
Emilio.