Getting different results from Equivalent Fourier Integrals in FEniCS

mathematics
1

Hello everyone,

I’m using FEniCS (2019 version) to compute a Fourier-type integral of a function q. My first approach is to directly integrate against a complex exponential:

def exp_Fourier(x, xi):
return np.exp(-2j * (xi[0] * x[0] + xi[1] * x[1]))

q_fourier_1 = assemble(exp_proj * q_proj * dx)

I then tried an alternative formulation, which should be mathematically equivalent:

Function V_0(x,y):

def V_0(x, xi, zeta):
return np.exp(-(zeta[0] + xi[0]*1j)*x[0] - (zeta[1] + xi[1]*1j)*x[1])

Function V(x,y):

def V(x, xi, zeta):
return np.exp(((zeta[0] - xi[0]*1j)*x[0] + (zeta[1] - xi[1]*1j)*x[1]) / 2)

q_fourier_2 = assemble(V0_proj * V_proj * V_proj * q_proj * dx)

Even though both expressions should yield the same result in theory, q_fourier_1 and q_fourier_2 produce different outcomes.

Has anyone experienced this kind of issue before? Could it be related to how FEniCS (2019) handles complex expressions or projection? Any insight would be greatly appreciated!

Thanks in advance!

2

As your code is far from complete, it is really hard to give you any advice (No way anyone can reproduce the results, and various variables are undefined).

One thing that I can see already is that you are using complex numbers in the second formulation. Complex number support is not part of legacy FEniCS, but DOLFINx (see for instance:The Poisson problem with complex numbers — FEniCSx tutorial)

3

Hello,
Sorry for that, here is a completed code (for large values of xi the results are really different):

import numpy as np
from dolfin import *

Parameters

xi_val = [20, 15]
zeta_val = [2, 1]

x0_val = 0.0
y0_val = 0.0
d_val = 0.1

Define functions

def V_0(x, xi, zeta):
return np.exp((-(zeta[0] + xi[0]*1j)*x[0] - (zeta[1] + xi[1]*1j)*x[1]) )

def V_1(x, xi, zeta):
return np.exp(((zeta[0] - xi[0]*1j)*x[0] + (zeta[1] - xi[1]*1j)*x[1])/2)

def exp_Fourier(x, xi):
return np.exp(-21j(xi[0] * x[0] + xi[1] * x[1]))

Function Xi:

def phi(x, d, x0, y0):
norm = np.sqrt((x[0]-x0)**2 + (x[1]-y0)**2)
if norm < d:
return np.exp(1+1/((norm/d)**2-1))
else:
return 0.0

Function q_epsilon:

def q_function_epsilon(x, d, x0, y0):
return phi(x, d, x0, y0)

def exp_Fourier(x, xi):
return np.exp(-21j(xi[0] * x[0] + xi[1] * x[1]))

mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 32, 32)
V = FunctionSpace(mesh, “CG”, 3)

— Process q_function—

class PythonQExpression(UserExpression):
def init(self, d, x0, y0, **kwargs):
self.d = d
self.x0 = x0
self.y0 = y0
super().init(**kwargs)

def eval(self, values, x):
    values[0] = q_function_epsilon(x, self.d, self.x0, self.y0)

def value_shape(self):
    return ()

q_expr = PythonQExpression(d=d_val, x0=x0_val, y0=y0_val, degree=2)
q_proj = interpolate(q_expr, V)

— Process V0 —

class PythonV0Expression(UserExpression):
def init(self, xi, zeta, **kwargs):
self.xi = xi
self.zeta = zeta
super().init(**kwargs)

def eval(self, values, x):
    values[0] = V_0(x, self.xi, self.zeta)
def value_shape(self):
    return ()

V0_expr = PythonV0Expression(xi=xi_val, zeta=zeta_val)
V0_proj = interpolate(V0_expr, V)

— Process V —

class PythonVExpression(UserExpression):
def init(self, xi, zeta, **kwargs):
self.xi = xi
self.zeta = zeta
super().init(**kwargs)
def eval(self, values, x):
values[0] = V_1(x, self.xi, self.zeta)
def value_shape(self):
return ()
V_expr = PythonVExpression(xi=xi_val, zeta=zeta_val)
V_proj = interpolate(V_expr, V)

— Process exp_Fourier —

class PythonQFExpression(UserExpression):
def init(self, xi, **kwargs):
self.xi = xi
super().init(**kwargs)
def eval(self, values, x):
values[0] = exp_Fourier(x, self.xi)
def value_shape(self):
return ()
F_q_expr = PythonQFExpression(xi=xi_val)
F_q_proj = interpolate(F_q_expr, V)

q_fourier_1 = assemble(V0_proj * V_proj * V_proj * q_proj * dx)
q_fourier_2 = assemble(F_q_proj * q_proj * dx)

print(“q_fourier_1:”, q_fourier_1)
print(“q_fourier_2:”, q_fourier_2)

thank you for your help.

4

First of all, please format your code with 3x`, ie

```python

#add code here
```

So that the code can be copy-pasted.

5 
import numpy as np
from dolfin import *


xi_val = [2, 1]
zeta_val = [2, 1]

x0_val = 0.0    
y0_val = 0.0
d_val = 0.1


def V_0(x, xi, zeta):
    return np.exp((-(zeta[0] + xi[0]*1j)*x[0] - (zeta[1] + xi[1]*1j)*x[1]) )

def V_1(x, xi, zeta):
    return np.exp(((zeta[0] - xi[0]*1j)*x[0] + (zeta[1] - xi[1]*1j)*x[1])/2)

def exp_Fourier(x, xi):
    return np.exp(-2*1j*(xi[0] * x[0] + xi[1] * x[1]))


def phi(x, d, x0, y0):
    norm = np.sqrt((x[0]-x0)**2 + (x[1]-y0)**2)
    if norm < d:
        return np.exp(1 + 1 / ((norm/d)**2 - 1))
    else:
        return 0.0


def q_function_epsilon(x, d, x0, y0):
    return phi(x, d, x0, y0)

mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 32, 32)
V = FunctionSpace(mesh, "CG", 3)


class PythonQExpression(UserExpression):
    def __init__(self, d, x0, y0, **kwargs):
        self.d = d
        self.x0 = x0
        self.y0 = y0
        super().__init__(**kwargs)

    def eval(self, values, x):
        values[0] = q_function_epsilon(x, self.d, self.x0, self.y0)

    def value_shape(self):
        return ()

q_expr = PythonQExpression(d=d_val, x0=x0_val, y0=y0_val, degree=2)
q_proj = interpolate(q_expr, V)


class PythonV0Expression(UserExpression):
    def __init__(self, xi, zeta, **kwargs):
        self.xi = xi
        self.zeta = zeta
        super().__init__(**kwargs)

    def eval(self, values, x):
        values[0] = V_0(x, self.xi, self.zeta)

    def value_shape(self):
        return ()

V0_expr = PythonV0Expression(xi=xi_val, zeta=zeta_val, degree=3)
V0_proj = interpolate(V0_expr, V)


class PythonVExpression(UserExpression):
    def __init__(self, xi, zeta, **kwargs):
        self.xi = xi
        self.zeta = zeta
        super().__init__(**kwargs)

    def eval(self, values, x):
        values[0] = V_1(x, self.xi, self.zeta)

    def value_shape(self):
        return ()

V_expr = PythonVExpression(xi=xi_val, zeta=zeta_val, degree=3)
V_proj = interpolate(V_expr, V)


class PythonQFExpression(UserExpression):
    def __init__(self, xi, **kwargs):
        self.xi = xi
        super().__init__(**kwargs)

    def eval(self, values, x):
        values[0] = exp_Fourier(x, self.xi)

    def value_shape(self):
        return ()

F_q_expr = PythonQFExpression(xi=xi_val, degree=3)
F_q_proj = interpolate(F_q_expr, V)


q_fourier_1 = assemble(V0_proj * V_proj * V_proj * q_proj * dx)
q_fourier_2 = assemble(F_q_proj * q_proj * dx)

print("q_fourier_1:", q_fourier_1)
print("q_fourier_2:", q_fourier_2)