Solving PDE with non-trivial null space

Main Forum Old-Forum - do not post here
1

Hi:

I am solving the following problem (solvability condition satisfied)

"a(u,v)- lam* b(u,v) = L(v),

where lam (non-zero) is the eigenvalue and the non-trivial null space is spanned by gfu, i.e. a(gfu,v) = lam* b(gfu,v)."

I assume that one cannot simply assemble the bilinear form "a(u,v)- lam* b(u,v) ", since lam is an eigenvalue. Is there a way to solve this PDE in netgen?

Any help is greatly appreciated.

Best,

2

Could this be realized through defining a FE space incorporating that the solution u is orthogonal to gfu? If yes, does netgen provides such modifications?

3

I don’t think you can do this purely in NGSolve without some C++ hacking, but if efficiency is not the key point you can use NGSolve to assemble the matrices, then convert them to scipy matrices and do the orthogonalization and solving in scipy…

4

Thanks for the reply. Efficiency is not an issue here :slight_smile:

5

Hello,

I think we need more information for this problem.

Have you computed the Eigenvalue/Eigenvector in advance, or should it be part of the same equation ?
Are your matrices symmetric and positive definite ?

If you know the Eigensystem, you can pose your equation orthogonal to the Eigenvector using a scalar Lagrange parameter (from NumberSpace).

Joachim

6

Hi Joachim:

Thank you very much for the reply. I have computed the eigenvalues/eigenfunctions in advance using Arnoldi solver; the matrices are symmetric and non-negative definite.

I am not sure yet how to use a scalar Lagrange parameter (from NumberSpace), is there a short illustrative example on this?

Best,
Shixu

7

Ok I’m sorry I was wrong… With Joachims comment I think it’s actually quite easy…
I think this should do what you want:

from netgen.geom2d import unit_square
from ngsolve import *

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

h1=H1(mesh,order=5,dirichlet="bottom|left|right|top",complex=True)

u,v = h1.TnT()

a = BilinearForm(h1)
a += SymbolicBFI(grad(u)*grad(v))
b = BilinearForm(h1)
b += SymbolicBFI(u*v)


a.Assemble()
b.Assemble()

evecs = GridFunction(h1,multidim=50)
vlam = ArnoldiSolver(a.mat,b.mat,h1.FreeDofs(),evecs.vecs,1)
Draw(evecs)

# get a non multidim GF out of the multidim one
evec0 = GridFunction(h1)
evec0.vec.data = evecs.vecs[0]

number = NumberSpace(mesh,complex=True)
fes = FESpace([h1,number])

(u,unum), (v,vnum) = fes.TnT()

# The orthogonality condition is satisfied because vnum = 1 means that evec0 * u = 0
# unum incorporates the evec0 part of the right hand side, if the argument of f is orthogonal to
# evec0, then unum should be 0
alam = BilinearForm(fes)
alam += SymbolicBFI(grad(u) * grad(v) - vlam[0] * u * v + unum * v * evec0)
alam += SymbolicBFI(evec0 * u * vnum)

f = LinearForm(fes)
f += SymbolicLFI(1*v)

alam.Assemble()
f.Assemble()

sol = GridFunction(fes,"solution")
sol.vec.data = alam.mat.Inverse(fes.FreeDofs()) * f.vec

Draw(sol.components[0])
8

This is great! I appreciate this short example and thank you very much!